Motivation for study

Glucocorticoids are named for their role in several components of glucose metabolism. However, this role is especially apparent in fasted state of lab rodents. The accute stress response, mediated by corticosterone, is often interpreted as mobilizing glucose to allow animals to cope with short term challenges. But, how general is this to different taxa (birds) and under different conditions (i.e., not after controlled fasting in the lab)? Relatively few studies that directly test these relationships under natural varying conditions.

Understanding whether cort is directly causally responsible for mobilizing glucose is important for thinking about evolution of the stress response. Are the strongest cort responders also the strongest glucose responders and what does that mean for other components of stress response?

Recent papers have emphasized the fact that the response to stress is multi-faceted and that individual components are often not closely correlated within individuals (e.g. Lattin’s receptor work and citations in Romero’s SICB paper). This is really important for thinking about adaptive value of differences in cort response and how stress response evolves. Need a better understanding of causal relationships between different components of stress response under natural conditions.

Questions: Population level 1. Does glucose increase with acute corticosterone response? (yes: adult & nestling) 2. Does this response differ in populations with different energy demands? (no: adult) 3. Does this response differ at breeding stages with different demands? (no: adults) 4. Does experimental downregulation of cort via dex also reduce glucose? (no: adults & nestlings) 5. Does experimental upregulation of cort via ACTH increase glucose further? (maybe: adult, no: nestling) 6. Are there differences by age and sex? (not sure if we care to include these, I haven’t plotted them) Individual level 7. Does circulating cort (b/s/d/a) correlate with circulating glucose (b/s/d/a)? (no: adults & nestlings) 8. Does magnitude of cort change correlate with magnitude of glucose change? (no: adults & nestlings)

General conclusions (I think):

We can say pretty clearly that, as expected, an accute stress response (handling) results in an increase in both corticosterone and circulating glucose for both adults and nestlings in different populations and at different breeding stages. That is generally consistent with the well known mechanisms linking cort and glucose. However, there is no evidence that at an individual level, differences in cort response are actually correlated with differences in glucose response regardless of context. So what is happening here? Integrate with recent work on multi-component nature of stress response and importance of different measures. Implications are that strong responders in one component of stress response cannot be assumed to be strong responders in another component. This is important for thinking about how selection can operate on the stress response.

Methods & data

  1. Adults were sampled during early incubation, late incubation (not much), or provisioning. Generally, the early and late captures included measures of corticosterone and glucose at baseline and stress induced time points. In 2019 NY, we also have glucose measured post-dex and–for a subset of females at third capture only–at post ACTH injections.

  2. Treatments: All years include some treatments to adults and we can deal with these in different ways. One option is to focus only on control groups and ignore treatment. That is nice and clean and sample size is still very large, but has the disadvantage of removing almost all of the ACTH adults since that was only done in 2019.

  3. Nestlings are from 2019 only. Include base, stress, dex, then 3 days later acth injection (so a bit awkward to explain that as stress to acth).

  4. The other populations don’t have post dex or post acth. I’m on the fence about whether it adds much to include the other populations or would be more streamlined to just focus on NY and have one less thing to have to describe. But they do fit fine with the general (lack of) results.

  5. We can also add in a set of glucose validation data. Or maybe I already described that in the plumage manipulation bioarxiv paper and we can just cite it here? We have about 30 birds with glucose measured twice showing that the method is highly repeatable.

Analysis plan

I think the best approach is to just have a few large models. Maybe one for nestlings. One for NY only adults (which includes dex & acth). And one for cross populations (which includes breeding stage). Three tables would have output from those models and then figures to illustrate teh main things. My leaning would be to just specify and fit the full model for each of those and then interpret everythign from there. There could be a supplement that fits the same things with or without the non-control birds depending on what we decide makes the most sense. In general, I think the analysis could and should be kept pretty simple (even though the models themselves will be pretty large). The figures can do a pretty good job of illustrating everything, I think. There will be some decisions about which figures make most sense to include adn which to maybe cut or stick in a supplement.

Results: Plots

Reading in some data and organizing to get it all ready for plotting.

######################################################################
## Script to analyze samples of glucose and corticosterone from     ##
## tree swallow adult and nestlings collected in NY, TN, WY, & AK   ##
## from 2016-2019. The main purpose is to ask if increase in cort   ##
## is directly correlated with increase in glucose. Different       ##
## subsets of data are used for different questions becaues not     ##
## all samples were always collected and some birds were part of    ##
## manipulative experiments that could have influenced measures.    ##
##                                                                  ##
## Code by Conor Taff, last updated 25 February 2020                ##
######################################################################

## Some differences in samples collected by age, year, and location
  # AK, TN, WY: no dex or acth glucose; no nestlings
  # NY: dex & acth glucose in some years; acth only in 2019 and only
      # on 3rd caputure of females. Some years second or third capture
      # did not have glucose measures.
  # Nestlings: only from NY in 2019. acth is 3 days after b/s/d series

## Load in packages to be used for various purposes.
  pacman::p_load(plyr, lme4, ggplot2, here, scales, lmerTest, sjPlot)

## Load in main data file that includes all glucose & cort measures
  d <- read.delim(here("/1_raw_data/data_glucose_cort.txt"))
  
## Substituting 0.47 as the minimum detectable amount of corticosterone.
  # In some cases this was done manually when running ELISA, but this is
  # checking to make sure all are set the same way here.
  
  d[which(d$b_cort < 0.47), "b_cort"] <- 0.47
  d[which(d$s_cort < 0.47), "s_cort"] <- 0.47
  d[which(d$d_cort < 0.47), "d_cort"] <- 0.47
  
## Calculate delta values for cort and glucose
  d$s_resp <- d$s_cort - d$b_cort
  d$n_feed <- d$d_cort - d$s_cort
  d$a_inc  <- d$a_cort - d$s_cort
  
  d$gluc_resp <- d$s_gluc - d$b_gluc
  d$gluc_feed <- d$d_gluc - d$s_gluc
  d$gluc_ainc <- d$a_gluc - d$s_gluc
  
## Make a separate data set for adults and nestlings and within adults
  # make one that excludes all post-treatment birds and one that excludes
  # birds from all but NY. These will be used for different analyses
    da <- subset(d, d$class == "adult")   
    da2 <- subset(da, da$post_trt == "no") # excludes all post treatment measures
    da2n <- subset(da2, da2$state == "NY") # New York only
    dn <- subset(d, d$class == "nestling")
    
    da2n <- subset(da, da$state == "NY")
    
## set colors to be used for all plotting
    b_col <- "#56B4E9"  # a light blue
    s_col <- "#E69F00"  # an orange
    d_col <- "#009E73"  # a dark green
    a_col <- "#F0E442"  # yellow

By time NY adult

First, lets look at corticosterone at different timepoints and glucose at different timepoints. This is just plotting raw data and not controlling for anything like population or breeding stage, etc.

## Make a basic plot for NY adults showing change in cort and change in glucose
    
    par(mfrow = c(1, 2))
    
    plot(1, 1, type = "n", yaxt = "n", xaxt = "n", xlab = "", 
        ylab = expression(paste("Corticosterone (ng/", mu, "l)")), bty = "n",
        xlim = c(1.5, 5.5), ylim = c(0, 103), xaxs = "i", yaxs = "i",
        main = "Adult New York Corticosterone")
    axis(1, c(-10, 2, 3, 4, 5, 10), c("", "Baseline", "Stress", "Negative", "ACTH", ""))  
    mtext("Induced", 1, at = 3, line = 2)
    mtext("Feedback", 1, at = 4, line = 2)
    mtext("Challenge", 1, at = 5, line = 2)
    axis(2, seq(-20, 160, 20), las = 2)
    points(rep(2, nrow(da2n)) + runif(nrow(da2n), -.1, .1), da2n$b_cort, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    points(rep(3, nrow(da2n)) + runif(nrow(da2n), -.1, .1), da2n$s_cort, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    points(rep(4, nrow(da2n)) + runif(nrow(da2n), -.1, .1), da2n$d_cort, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    points(rep(5, nrow(da2n)) + runif(nrow(da2n), -.1, .1), da2n$a_cort, pch = 16,
           col = alpha("gray40", 0.6), cex = .5)
    boxplot(da2n$b_cort, add = TRUE, at = 2, axes = FALSE, outline = FALSE, col = alpha(b_col, 0.6))
    boxplot(da2n$s_cort, add = TRUE, at = 3, axes = FALSE, outline = FALSE, col = alpha(s_col, 0.6))
    boxplot(da2n$d_cort, add = TRUE, at = 4, axes = FALSE, outline = FALSE, col = alpha(d_col, 0.6))
    boxplot(da2n$a_cort, add = TRUE, at = 5, axes = FALSE, outline = FALSE, col = alpha(a_col, 0.6))
    
    
    plot(1, 1, type = "n", yaxt = "n", xaxt = "n", xlab = "", 
         ylab = "Blood Glucose (mg/dl)", bty = "n",
         xlim = c(1.5, 5.5), ylim = c(80, 420), xaxs = "i", yaxs = "i",
         main = "Adult New York Glucose")
    axis(1, c(-10, 2, 3, 4, 5, 10), c("", "Baseline", "Stress", "Negative", "ACTH", ""))  
    mtext("Induced", 1, at = 3, line = 2)
    mtext("Feedback", 1, at = 4, line = 2)
    mtext("Challenge", 1, at = 5, line = 2)
    axis(2, seq(-50, 800, 50), las = 2)
    points(rep(2, nrow(da2n)) + runif(nrow(da2n), -.1, .1), da2n$b_gluc, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    points(rep(3, nrow(da2n)) + runif(nrow(da2n), -.1, .1), da2n$s_gluc, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    points(rep(4, nrow(da2n)) + runif(nrow(da2n), -.1, .1), da2n$d_gluc, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    points(rep(5, nrow(da2n)) + runif(nrow(da2n), -.1, .1), da2n$a_gluc, pch = 16,
           col = alpha("gray40", 0.6), cex = 0.5)
    boxplot(da2n$b_gluc, add = TRUE, at = 2, axes = FALSE, outline = FALSE, col = alpha(b_col, 0.6))
    boxplot(da2n$s_gluc, add = TRUE, at = 3, axes = FALSE, outline = FALSE, col = alpha(s_col, 0.6))
    boxplot(da2n$d_gluc, add = TRUE, at = 4, axes = FALSE, outline = FALSE, col = alpha(d_col, 0.6))
    boxplot(da2n$a_gluc, add = TRUE, at = 5, axes = FALSE, outline = FALSE, col = alpha(a_col, 0.6))

The population level pattern here is pretty clear and obvious for cort. One caveat for ACTH is that those are all from only 2019 and only provisioning, but it does seem to increase cort somewhat more and certainly more than dex injection. Glucose clearly goes up from base to stress, and maybe from stress to acth (small sample), but is the same at stress and post-dex.

This is easier to see if you plot the delta values for cort and glucose.

## Make a basic plot for NY adults showing change in cort and change in glucose
    
    par(mfrow = c(1, 2))
    
    plot(1, 1, type = "n", yaxt = "n", xaxt = "n", xlab = "", 
         ylab = expression(paste(Delta, "Corticosterone (ng/", mu, "l)")), bty = "n",
         xlim = c(1.5, 4.5), ylim = c(-100, 100), xaxs = "i", yaxs = "i",
         main = "Adult New York Delta Corticosterone")
    axis(1, c(-10, 2, 3, 4, 10), c("", "Base to", "Stress to", "Stress to", ""))  
    mtext("Stress", 1, at = 2, line = 2)
    mtext("Dex", 1, at = 3, line = 2)
    mtext("ACTH", 1, at = 4, line = 2)
    axis(2, seq(-500, 500, 50), las = 2)
    abline(h = 0, lty = 2)
    points(rep(2, nrow(da2n)) + runif(nrow(da2n), -.1, .1), da2n$s_resp, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    points(rep(3, nrow(da2n)) + runif(nrow(da2n), -.1, .1), da2n$n_feed, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    points(rep(4, nrow(da2n)) + runif(nrow(da2n), -.1, .1), da2n$a_inc, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    boxplot(da2n$s_resp, add = TRUE, at = 2, axes = FALSE, outline = FALSE, col = alpha(s_col, 0.6))
    boxplot(da2n$n_feed, add = TRUE, at = 3, axes = FALSE, outline = FALSE, col = alpha(d_col, 0.6))
    boxplot(da2n$a_inc, add = TRUE, at = 4, axes = FALSE, outline = FALSE, col = alpha(a_col, 0.6))
    
    
    plot(1, 1, type = "n", yaxt = "n", xaxt = "n", xlab = "", 
         ylab = expression(paste(Delta, "Blood Glucose (mg/dl)")), bty = "n",
         xlim = c(1.5, 4.5), ylim = c(-100, 150), xaxs = "i", yaxs = "i",
         main = "Adult New York Delta Glucose")
    axis(1, c(-10, 2, 3, 4, 10), c("", "Base to", "Stress to", "Stress to", ""))  
    mtext("Stress", 1, at = 2, line = 2)
    mtext("Dex", 1, at = 3, line = 2)
    mtext("ACTH", 1, at = 4, line = 2)
    axis(2, seq(-500, 500, 50), las = 2)
    abline(h = 0, lty = 2)
    points(rep(2, nrow(da2n)) + runif(nrow(da2n), -.1, .1), da2n$gluc_resp, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    points(rep(3, nrow(da2n)) + runif(nrow(da2n), -.1, .1), da2n$gluc_feed, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    points(rep(4, nrow(da2n)) + runif(nrow(da2n), -.1, .1), da2n$gluc_ainc, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    boxplot(da2n$gluc_resp, add = TRUE, at = 2, axes = FALSE, outline = FALSE, col = alpha(s_col, 0.6))
    boxplot(da2n$gluc_feed, add = TRUE, at = 3, axes = FALSE, outline = FALSE, col = alpha(d_col, 0.6))
    boxplot(da2n$gluc_ainc, add = TRUE, at = 4, axes = FALSE, outline = FALSE, col = alpha(a_col, 0.6))

By time NY nestlings

Now do the same two plots as above but for nestlings in NY.

 par(mfrow = c(1, 2))
    
    plot(1, 1, type = "n", yaxt = "n", xaxt = "n", xlab = "", 
         ylab = expression(paste("Corticosterone (ng/", mu, "l)")), bty = "n",
         xlim = c(1.5, 5.5), ylim = c(0, 85), xaxs = "i", yaxs = "i",
         main = "Nestling New York Corticosterone")
    axis(1, c(-10, 2, 3, 4, 5, 10), c("", "Baseline", "Stress", "Negative", "ACTH", ""))  
    mtext("Induced", 1, at = 3, line = 2)
    mtext("Feedback", 1, at = 4, line = 2)
    mtext("Challenge", 1, at = 5, line = 2)
    axis(2, seq(-20, 160, 20), las = 2)
    points(rep(2, nrow(dn)) + runif(nrow(dn), -.1, .1), dn$b_cort, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    points(rep(3, nrow(dn)) + runif(nrow(dn), -.1, .1), dn$s_cort, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    points(rep(4, nrow(dn)) + runif(nrow(dn), -.1, .1), dn$d_cort, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    points(rep(5, nrow(dn)) + runif(nrow(dn), -.1, .1), dn$a_cort, pch = 16,
           col = alpha("gray40", 0.6), cex = .5)
    boxplot(dn$b_cort, add = TRUE, at = 2, axes = FALSE, outline = FALSE, col = alpha(b_col, 0.6))
    boxplot(dn$s_cort, add = TRUE, at = 3, axes = FALSE, outline = FALSE, col = alpha(s_col, 0.6))
    boxplot(dn$d_cort, add = TRUE, at = 4, axes = FALSE, outline = FALSE, col = alpha(d_col, 0.6))
    boxplot(dn$a_cort, add = TRUE, at = 5, axes = FALSE, outline = FALSE, col = alpha(a_col, 0.6))
    
    
    plot(1, 1, type = "n", yaxt = "n", xaxt = "n", xlab = "", 
         ylab = "Blood Glucose (mg/dl)", bty = "n",
         xlim = c(1.5, 5.5), ylim = c(80, 420), xaxs = "i", yaxs = "i",
         main = "Nestling New York Glucose")
    axis(1, c(-10, 2, 3, 4, 5, 10), c("", "Baseline", "Stress", "Negative", "ACTH", ""))  
    mtext("Induced", 1, at = 3, line = 2)
    mtext("Feedback", 1, at = 4, line = 2)
    mtext("Challenge", 1, at = 5, line = 2)
    axis(2, seq(-50, 800, 50), las = 2)
    points(rep(2, nrow(dn)) + runif(nrow(dn), -.1, .1), dn$b_gluc, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    points(rep(3, nrow(dn)) + runif(nrow(dn), -.1, .1), dn$s_gluc, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    points(rep(4, nrow(dn)) + runif(nrow(dn), -.1, .1), dn$d_gluc, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    points(rep(5, nrow(dn)) + runif(nrow(dn), -.1, .1), dn$a_gluc, pch = 16,
           col = alpha("gray40", 0.6), cex = 0.5)
    boxplot(dn$b_gluc, add = TRUE, at = 2, axes = FALSE, outline = FALSE, col = alpha(b_col, 0.6))
    boxplot(dn$s_gluc, add = TRUE, at = 3, axes = FALSE, outline = FALSE, col = alpha(s_col, 0.6))
    boxplot(dn$d_gluc, add = TRUE, at = 4, axes = FALSE, outline = FALSE, col = alpha(d_col, 0.6))
    boxplot(dn$a_gluc, add = TRUE, at = 5, axes = FALSE, outline = FALSE, col = alpha(a_col, 0.6))
    
     par(mfrow = c(1, 2))
    
    plot(1, 1, type = "n", yaxt = "n", xaxt = "n", xlab = "", 
         ylab = expression(paste(Delta, "Corticosterone (ng/", mu, "l)")), bty = "n",
         xlim = c(1.5, 4.5), ylim = c(-100, 100), xaxs = "i", yaxs = "i",
         main = "Nestling New York Delta Corticosterone")
    axis(1, c(-10, 2, 3, 4, 10), c("", "Base to", "Stress to", "Stress to", ""))  
    mtext("Stress", 1, at = 2, line = 2)
    mtext("Dex", 1, at = 3, line = 2)
    mtext("ACTH", 1, at = 4, line = 2)
    axis(2, seq(-500, 500, 50), las = 2)
    abline(h = 0, lty = 2)
    points(rep(2, nrow(dn)) + runif(nrow(dn), -.1, .1), dn$s_resp, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    points(rep(3, nrow(dn)) + runif(nrow(dn), -.1, .1), dn$n_feed, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    points(rep(4, nrow(dn)) + runif(nrow(dn), -.1, .1), dn$a_inc, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    boxplot(dn$s_resp, add = TRUE, at = 2, axes = FALSE, outline = FALSE, col = alpha(s_col, 0.6))
    boxplot(dn$n_feed, add = TRUE, at = 3, axes = FALSE, outline = FALSE, col = alpha(d_col, 0.6))
    boxplot(dn$a_inc, add = TRUE, at = 4, axes = FALSE, outline = FALSE, col = alpha(a_col, 0.6))
    
    
    plot(1, 1, type = "n", yaxt = "n", xaxt = "n", xlab = "", 
         ylab = expression(paste(Delta, "Blood Glucose (mg/dl)")), bty = "n",
         xlim = c(1.5, 4.5), ylim = c(-160, 130), xaxs = "i", yaxs = "i",
         main = "Nestling New York Delta Glucose")
    axis(1, c(-10, 2, 3, 4, 10), c("", "Base to", "Stress to", "Stress to", ""))  
    mtext("Stress", 1, at = 2, line = 2)
    mtext("Dex", 1, at = 3, line = 2)
    mtext("ACTH", 1, at = 4, line = 2)
    axis(2, seq(-500, 500, 50), las = 2)
    abline(h = 0, lty = 2)
    points(rep(2, nrow(dn)) + runif(nrow(dn), -.1, .1), dn$gluc_resp, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    points(rep(3, nrow(dn)) + runif(nrow(dn), -.1, .1), dn$gluc_feed, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    points(rep(4, nrow(dn)) + runif(nrow(dn), -.1, .1), dn$gluc_ainc, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    boxplot(dn$gluc_resp, add = TRUE, at = 2, axes = FALSE, outline = FALSE, col = alpha(s_col, 0.6))
    boxplot(dn$gluc_feed, add = TRUE, at = 3, axes = FALSE, outline = FALSE, col = alpha(d_col, 0.6))
    boxplot(dn$gluc_ainc, add = TRUE, at = 4, axes = FALSE, outline = FALSE, col = alpha(a_col, 0.6))

The patterns for nestlings are basically exactly the same as for adults. Cort does just what we would expect. Glucose goes up from base to stress, but does not change at post-dex or post-acth. Note that the lack of change post-acth is a bit different than adults although the adult sample size is very small.

Breeding stage

Comparison of glucose values and delta at different breeding stages. This is combining the different states, though I don’t think it would look any different in each state. The idea here is that certain stages that are more like ‘fasting’ might show clearer glucose responses.

## Breeding stage
    
    par(mfrow = c(1, 3))
    
    da_e <- subset(da2n, da2n$stage == "early_inc")
    da_l <- subset(da2n, da2n$stage == "late_inc")
    da_p <- subset(da2n, da2n$stage == "provision")
    
    plot(1, 1, type = "n", yaxt = "n", xaxt = "n", xlab = "", 
         ylab = "Blood Glucose (mg/dl)", bty = "n",
         xlim = c(1.5, 4.5), ylim = c(100, 350), xaxs = "i", yaxs = "i",
         main = "Base Glucose")
    axis(1, c(-10, 2, 3, 4, 10), c("", "Early Inc", "Late Inc", "Provision", ""))  
    axis(2, seq(-500, 500, 50), las = 2)
    points(rep(2, nrow(da_e)) + runif(nrow(da_e), -.1, .1), da_e$b_gluc, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    points(rep(3, nrow(da_l)) + runif(nrow(da_l), -.1, .1), da_l$b_gluc, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    points(rep(4, nrow(da_p)) + runif(nrow(da_p), -.1, .1), da_p$b_gluc, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    boxplot(da_e$b_gluc, add = TRUE, at = 2, axes = FALSE, outline = FALSE, col = alpha("chartreuse3", 0.6))
    boxplot(da_l$b_gluc, add = TRUE, at = 3, axes = FALSE, outline = FALSE, col = alpha("chartreuse3", 0.6))
    boxplot(da_p$b_gluc, add = TRUE, at = 4, axes = FALSE, outline = FALSE, col = alpha("chartreuse3", 0.6))
    
    
    da_e <- subset(da2n, da2n$stage == "early_inc")
    da_l <- subset(da2n, da2n$stage == "late_inc")
    da_p <- subset(da2n, da2n$stage == "provision")
    
    plot(1, 1, type = "n", yaxt = "n", xaxt = "n", xlab = "", 
         ylab = "Blood Glucose (mg/dl)", bty = "n",
         xlim = c(1.5, 4.5), ylim = c(100, 350), xaxs = "i", yaxs = "i",
         main = "Stress Glucose")
    axis(1, c(-10, 2, 3, 4, 10), c("", "Early Inc", "", "Provision", ""))  
    axis(2, seq(-500, 500, 50), las = 2)
    points(rep(2, nrow(da_e)) + runif(nrow(da_e), -.1, .1), da_e$s_gluc, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    #points(rep(3, nrow(da_l)) + runif(nrow(da_l), -.1, .1), da_l$s_gluc, pch = 16, 
     #      col = alpha("gray40", 0.6), cex = .5)
    points(rep(4, nrow(da_p)) + runif(nrow(da_p), -.1, .1), da_p$s_gluc, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    boxplot(da_e$s_gluc, add = TRUE, at = 2, axes = FALSE, outline = FALSE, col = alpha("chartreuse3", 0.6))
    #boxplot(da_l$s_gluc, add = TRUE, at = 3, axes = FALSE, outline = FALSE, col = alpha("chartreuse3", 0.6))
    boxplot(da_p$s_gluc, add = TRUE, at = 4, axes = FALSE, outline = FALSE, col = alpha("chartreuse3", 0.6))
    
    
    plot(1, 1, type = "n", yaxt = "n", xaxt = "n", xlab = "", 
         ylab = "Blood Glucose (mg/dl)", bty = "n",
         xlim = c(1.5, 4.5), ylim = c(-100, 150), xaxs = "i", yaxs = "i",
         main = "Glucose Response")
    axis(1, c(-10, 2, 3, 4, 10), c("", "Early Inc", "", "Provision", ""))  
    axis(2, seq(-500, 500, 50), las = 2)
    abline(h = 0, lty = 2)
    points(rep(2, nrow(da_e)) + runif(nrow(da_e), -.1, .1), da_e$gluc_resp, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    #points(rep(3, nrow(da_l)) + runif(nrow(da_l), -.1, .1), da_l$s_gluc, pch = 16, 
    #      col = alpha("gray40", 0.6), cex = .5)
    points(rep(4, nrow(da_p)) + runif(nrow(da_p), -.1, .1), da_p$gluc_resp, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    boxplot(da_e$gluc_resp, add = TRUE, at = 2, axes = FALSE, outline = FALSE, col = alpha("chartreuse3", 0.6))
    #boxplot(da_l$s_gluc, add = TRUE, at = 3, axes = FALSE, outline = FALSE, col = alpha("chartreuse3", 0.6))
    boxplot(da_p$gluc_resp, add = TRUE, at = 4, axes = FALSE, outline = FALSE, col = alpha("chartreuse3", 0.6))

There really doesn’t seem to be any difference in any of these glucose measures in incubation vs. provisioning.

Population comparison

Maybe there is some difference across populations where resources or energy demands differ? This probably isn’t the most compelling part of the comparisons here and the results aren’t any different from NY alone. Not sure if it is worth including or not.

Firs plotting aboslute values:

## Between population comparisons, just base and stress (no dex glucose)
    
    par(mfrow = c(1, 2))
    dany <- subset(da2, da2$state == "NY")
    dawy <- subset(da2, da2$state == "WY")
    daak <- subset(da2, da2$state == "AK")
    datn <- subset(da2, da2$state == "TN")
    
    col_ny <- "red"
    col_tn <- "orange"
    col_wy <- "yellow"
    col_ak <- "slateblue"
    
    plot(1, 1, type = "n", yaxt = "n", xaxt = "n", xlab = "", 
         ylab = expression(paste("Corticosterone (ng/", mu, "l)")), bty = "n",
         xlim = c(1.5, 5.5), ylim = c(0, 103), xaxs = "i", yaxs = "i",
         main = "Population Comparison Corticosterone")
    axis(1, c(-10, 2.5, 4.5, 10), c("", "Baseline", "Stress", ""))  
    mtext("Induced", 1, at = 4.5, line = 2)
    axis(2, seq(-20, 160, 20), las = 2)
    points(rep(2, nrow(dany)) + runif(nrow(dany), -.1, .1), dany$b_cort, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    points(rep(2.33, nrow(datn)) + runif(nrow(datn), -.1, .1), datn$b_cort, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    points(rep(2.67, nrow(dawy)) + runif(nrow(dawy), -.1, .1), dawy$b_cort, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    points(rep(3, nrow(daak)) + runif(nrow(daak), -.1, .1), daak$b_cort, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    
    points(rep(4, nrow(dany)) + runif(nrow(dany), -.1, .1), dany$s_cort, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    points(rep(4.33, nrow(datn)) + runif(nrow(datn), -.1, .1), datn$s_cort, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    points(rep(4.67, nrow(dawy)) + runif(nrow(dawy), -.1, .1), dawy$s_cort, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    points(rep(5, nrow(daak)) + runif(nrow(daak), -.1, .1), daak$s_cort, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    
    boxplot(dany$b_cort, add = TRUE, at = 2, axes = FALSE, outline = FALSE, col = alpha(col_ny, 0.6), boxwex = .5)
    boxplot(datn$b_cort, add = TRUE, at = 2.33, axes = FALSE, outline = FALSE, col = alpha(col_tn, 0.6), boxwex = .5)
    boxplot(dawy$b_cort, add = TRUE, at = 2.67, axes = FALSE, outline = FALSE, col = alpha(col_wy, 0.6), boxwex = .5)
    boxplot(daak$b_cort, add = TRUE, at = 3, axes = FALSE, outline = FALSE, col = alpha(col_ak, 0.6), boxwex = .5)
    
    
    boxplot(dany$s_cort, add = TRUE, at = 4, axes = FALSE, outline = FALSE, col = alpha(col_ny, 0.6), boxwex = .5)
    boxplot(datn$s_cort, add = TRUE, at = 4.33, axes = FALSE, outline = FALSE, col = alpha(col_tn, 0.6), boxwex = .5)
    boxplot(dawy$s_cort, add = TRUE, at = 4.67, axes = FALSE, outline = FALSE, col = alpha(col_wy, 0.6), boxwex = .5)
    boxplot(daak$s_cort, add = TRUE, at = 5, axes = FALSE, outline = FALSE, col = alpha(col_ak, 0.6), boxwex = .5)
    
    legend("topleft", c("NY", "TN", "WY", "AK"), bty = "n", pch = 22, pt.bg = c(col_ny,
        col_tn, col_wy, col_ak), pt.cex = 2)
    
    plot(1, 1, type = "n", yaxt = "n", xaxt = "n", xlab = "", 
         ylab = "Blood Glucose (mg/dl)", bty = "n",
         xlim = c(1.5, 5.5), ylim = c(80, 390), xaxs = "i", yaxs = "i",
         main = "Population Comparison Glucose")
    axis(1, c(-10, 2.5, 4.5, 10), c("", "Baseline", "Stress", ""))  
    mtext("Induced", 1, at = 4.5, line = 2)
    axis(2, seq(0, 500, 50), las = 2)
    points(rep(2, nrow(dany)) + runif(nrow(dany), -.1, .1), dany$b_gluc, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    points(rep(2.33, nrow(datn)) + runif(nrow(datn), -.1, .1), datn$b_gluc, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    points(rep(2.67, nrow(dawy)) + runif(nrow(dawy), -.1, .1), dawy$b_gluc, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    points(rep(3, nrow(daak)) + runif(nrow(daak), -.1, .1), daak$b_gluc, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    
    points(rep(4, nrow(dany)) + runif(nrow(dany), -.1, .1), dany$s_gluc, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    points(rep(4.33, nrow(datn)) + runif(nrow(datn), -.1, .1), datn$s_gluc, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    points(rep(4.67, nrow(dawy)) + runif(nrow(dawy), -.1, .1), dawy$s_gluc, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    points(rep(5, nrow(daak)) + runif(nrow(daak), -.1, .1), daak$s_gluc, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    
    boxplot(dany$b_gluc, add = TRUE, at = 2, axes = FALSE, outline = FALSE, col = alpha(col_ny, 0.6), boxwex = .5)
    boxplot(datn$b_gluc, add = TRUE, at = 2.33, axes = FALSE, outline = FALSE, col = alpha(col_tn, 0.6), boxwex = .5)
    boxplot(dawy$b_gluc, add = TRUE, at = 2.67, axes = FALSE, outline = FALSE, col = alpha(col_wy, 0.6), boxwex = .5)
    boxplot(daak$b_gluc, add = TRUE, at = 3, axes = FALSE, outline = FALSE, col = alpha(col_ak, 0.6), boxwex = .5)
    
    
    boxplot(dany$s_gluc, add = TRUE, at = 4, axes = FALSE, outline = FALSE, col = alpha(col_ny, 0.6), boxwex = .5)
    boxplot(datn$s_gluc, add = TRUE, at = 4.33, axes = FALSE, outline = FALSE, col = alpha(col_tn, 0.6), boxwex = .5)
    boxplot(dawy$s_gluc, add = TRUE, at = 4.67, axes = FALSE, outline = FALSE, col = alpha(col_wy, 0.6), boxwex = .5)
    boxplot(daak$s_gluc, add = TRUE, at = 5, axes = FALSE, outline = FALSE, col = alpha(col_ak, 0.6), boxwex = .5)

And then plotting delta values:

## Population delta glucose and delta cort
    
    
    par(mfrow = c(1, 2))
    dany <- subset(da2, da2$state == "NY")
    dawy <- subset(da2, da2$state == "WY")
    daak <- subset(da2, da2$state == "AK")
    datn <- subset(da2, da2$state == "TN")
    
    col_ny <- "red"
    col_tn <- "orange"
    col_wy <- "yellow"
    col_ak <- "slateblue"
    
    plot(1, 1, type = "n", yaxt = "n", xaxt = "n", xlab = "", 
         ylab = expression(paste(Delta, "Corticosterone (ng/", mu, "l)")), bty = "n",
         xlim = c(1.5, 3.5), ylim = c(-20, 200), xaxs = "i", yaxs = "i",
         main = "Population Comparison Corticosterone")
    axis(1, c(-10, 2.5, 10), c("", "Base to", ""))  
    mtext("Stress", 1, at = 2.5, line = 2)
    axis(2, seq(-50, 260, 50), las = 2)
    points(rep(2, nrow(dany)) + runif(nrow(dany), -.1, .1), dany$s_resp, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    points(rep(2.33, nrow(datn)) + runif(nrow(datn), -.1, .1), datn$s_resp, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    points(rep(2.67, nrow(dawy)) + runif(nrow(dawy), -.1, .1), dawy$s_resp, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    points(rep(3, nrow(daak)) + runif(nrow(daak), -.1, .1), daak$s_resp, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    
    boxplot(dany$s_resp, add = TRUE, at = 2, axes = FALSE, outline = FALSE, col = alpha(col_ny, 0.6), boxwex = .5)
    boxplot(datn$s_resp, add = TRUE, at = 2.33, axes = FALSE, outline = FALSE, col = alpha(col_tn, 0.6), boxwex = .5)
    boxplot(dawy$s_resp, add = TRUE, at = 2.67, axes = FALSE, outline = FALSE, col = alpha(col_wy, 0.6), boxwex = .5)
    boxplot(daak$s_resp, add = TRUE, at = 3, axes = FALSE, outline = FALSE, col = alpha(col_ak, 0.6), boxwex = .5)
    
    legend("topleft", c("NY", "TN", "WY", "AK"), bty = "n", pch = 22, pt.bg = c(col_ny,
                                                                                col_tn, col_wy, col_ak), pt.cex = 2)
    
    abline(h = 0, lty = 2)
    
    
    plot(1, 1, type = "n", yaxt = "n", xaxt = "n", xlab = "", 
         ylab = expression(paste(Delta, "Blood Glucose (mg/dl)")), bty = "n",
         xlim = c(1.5, 3.5), ylim = c(-150, 225), xaxs = "i", yaxs = "i",
         main = "Population Comparison Glucose")
    axis(1, c(-10, 2.5, 10), c("", "Base to", ""))  
    mtext("Stress", 1, at = 2.5, line = 2)
    axis(2, seq(-500, 5000, 50), las = 2)
    points(rep(2, nrow(dany)) + runif(nrow(dany), -.1, .1), dany$gluc_resp, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    points(rep(2.33, nrow(datn)) + runif(nrow(datn), -.1, .1), datn$gluc_resp, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    points(rep(2.67, nrow(dawy)) + runif(nrow(dawy), -.1, .1), dawy$gluc_resp, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    points(rep(3, nrow(daak)) + runif(nrow(daak), -.1, .1), daak$gluc_resp, pch = 16, 
           col = alpha("gray40", 0.6), cex = .5)
    
    boxplot(dany$gluc_resp, add = TRUE, at = 2, axes = FALSE, outline = FALSE, col = alpha(col_ny, 0.6), boxwex = .5)
    boxplot(datn$gluc_resp, add = TRUE, at = 2.33, axes = FALSE, outline = FALSE, col = alpha(col_tn, 0.6), boxwex = .5)
    boxplot(dawy$gluc_resp, add = TRUE, at = 2.67, axes = FALSE, outline = FALSE, col = alpha(col_wy, 0.6), boxwex = .5)
    boxplot(daak$gluc_resp, add = TRUE, at = 3, axes = FALSE, outline = FALSE, col = alpha(col_ak, 0.6), boxwex = .5)
    
    abline(h = 0, lty = 2)

I think the main takeaway here is that despite some pretty clear differences in corticosterone response, there aren’t really any clear differences in glucose response. All four populations do show an increase in glucose from base to stress though.

Individual cort-glucose relationships

So all of the above shows that glucose does go up with the stress response in general, but doesn’t say anything about whether individuals with a stronger corticosterone response also have a stronger glucose response. The scatterplots below are showing all of these relationships for adults and nestlings. This is limited to NY birds, but for now these plots aren’t controlling for anything else (treatment, breeding stage, etc).

First plotting circulating cort levels vs. circulating glucose levels at each timepoint for adults.

 ## Scatter plots of cort values and glucose values for adults
    par(mfrow = c(2, 2))
    
    plot(1, 1, type = "n", yaxt = "n", xaxt = "n", xlab = expression(paste("Corticosterone (ng/", mu, "l)")),
         bty = "n", xlim = c(0, 38), ylim = c(80, 370), xaxs = "i", yaxs = "i", 
         ylab = "Blood Glucose (mg/dl)", main = "New York Adults Baseline")
    axis(1, seq(-20, 200, 5))    
    axis(2, seq(0, 500, 50), las = 2)    
    points(da2n$b_cort, da2n$b_gluc, pch = 21, bg = alpha(b_col, 0.6), cex = 0.8)
    m <- lm(b_gluc ~ b_cort, data = da2n)
    newx = seq(min(na.omit(da2n$b_cort)), max(na.omit(da2n$b_cort)), 0.05)
    conf_interval <- predict(m, newdata = data.frame(b_cort = newx), interval = "confidence",
                             level = 0.95)
    polygon(c(rev(newx), newx), c(rev(conf_interval[, 3]), 
                                  conf_interval[, 2]), col = alpha("gray40", 0.4), border = NA)
    lines(newx, conf_interval[, 1], lwd = 2)
    
    plot(1, 1, type = "n", yaxt = "n", xaxt = "n", xlab = expression(paste("Corticosterone (ng/", mu, "l)")),
         bty = "n", xlim = c(0, 110), ylim = c(80, 370), xaxs = "i", yaxs = "i", 
         ylab = "Blood Glucose (mg/dl)", main = "New York Adults Stress Induced")
    axis(1, seq(-20, 200, 20))    
    axis(2, seq(0, 500, 50), las = 2)    
    points(da2n$s_cort, da2n$s_gluc, pch = 21, bg = alpha(s_col, 0.6), cex = 0.8)
    m <- lm(s_gluc ~ s_cort, data = da2n)
    newx = seq(min(na.omit(da2n$s_cort)), max(na.omit(da2n$s_cort)), 0.05)
    conf_interval <- predict(m, newdata = data.frame(s_cort = newx), interval = "confidence",
                             level = 0.95)
    polygon(c(rev(newx), newx), c(rev(conf_interval[, 3]), 
                                  conf_interval[, 2]), col = alpha("gray40", 0.4), border = NA)
    lines(newx, conf_interval[, 1], lwd = 2)
    
    plot(1, 1, type = "n", yaxt = "n", xaxt = "n", xlab = expression(paste("Corticosterone (ng/", mu, "l)")),
         bty = "n", xlim = c(0, 48), ylim = c(80, 370), xaxs = "i", yaxs = "i", 
         ylab = "Blood Glucose (mg/dl)", main = "New York Adults Post Dex")
    axis(1, seq(-20, 200, 10))    
    axis(2, seq(0, 500, 50), las = 2)    
    points(da2n$d_cort, da2n$d_gluc, pch = 21, bg = alpha(d_col, 0.6), cex = 0.8)
    m <- lm(d_gluc ~ d_cort, data = da2n)
    newx = seq(min(na.omit(da2n$d_cort)), max(na.omit(da2n$d_cort)), 0.05)
    conf_interval <- predict(m, newdata = data.frame(d_cort = newx), interval = "confidence",
                             level = 0.95)
    polygon(c(rev(newx), newx), c(rev(conf_interval[, 3]), 
                                  conf_interval[, 2]), col = alpha("gray40", 0.4), border = NA)
    lines(newx, conf_interval[, 1], lwd = 2)
    
    plot(1, 1, type = "n", yaxt = "n", xaxt = "n", xlab = expression(paste("Corticosterone (ng/", mu, "l)")),
         bty = "n", xlim = c(0, 110), ylim = c(80, 420), xaxs = "i", yaxs = "i", 
         ylab = "Blood Glucose (mg/dl)", main = "New York Adults Post ACTH")
    axis(1, seq(-20, 200, 10))    
    axis(2, seq(0, 500, 50), las = 2)    
    points(da2n$a_cort, da2n$a_gluc, pch = 21, bg = alpha(a_col, 0.6), cex = 0.8)
    m <- lm(a_gluc ~ a_cort, data = da2n)
    newx = seq(min(na.omit(da2n$a_cort)), max(na.omit(da2n$a_cort)), 0.05)
    conf_interval <- predict(m, newdata = data.frame(a_cort = newx), interval = "confidence",
                             level = 0.95)
    polygon(c(rev(newx), newx), c(rev(conf_interval[, 3]), 
                                  conf_interval[, 2]), col = alpha("gray40", 0.4), border = NA)
    lines(newx, conf_interval[, 1], lwd = 2)

And for nestlings:

## Scatter plots of cort values and glucose values for nestlings
    par(mfrow = c(2, 2))
    
    plot(1, 1, type = "n", yaxt = "n", xaxt = "n", xlab = expression(paste("Corticosterone (ng/", mu, "l)")),
         bty = "n", xlim = c(0, 32), ylim = c(80, 380), xaxs = "i", yaxs = "i", 
         ylab = "Blood Glucose (mg/dl)", main = "New York Nestling Baseline")
    axis(1, seq(-20, 200, 5))    
    axis(2, seq(0, 500, 50), las = 2)    
    points(dn$b_cort, dn$b_gluc, pch = 21, bg = alpha(b_col, 0.6), cex = 0.8)
    m <- lm(b_gluc ~ b_cort, data = dn)
    newx = seq(min(na.omit(dn$b_cort)), max(na.omit(dn$b_cort)), 0.05)
    conf_interval <- predict(m, newdata = data.frame(b_cort = newx), interval = "confidence",
                             level = 0.95)
    polygon(c(rev(newx), newx), c(rev(conf_interval[, 3]), 
                                  conf_interval[, 2]), col = alpha("gray40", 0.4), border = NA)
    lines(newx, conf_interval[, 1], lwd = 2)
    
    plot(1, 1, type = "n", yaxt = "n", xaxt = "n", xlab = expression(paste("Corticosterone (ng/", mu, "l)")),
         bty = "n", xlim = c(0, 110), ylim = c(80, 380), xaxs = "i", yaxs = "i", 
         ylab = "Blood Glucose (mg/dl)", main = "New York Nestling Stress Induced")
    axis(1, seq(-20, 200, 20))    
    axis(2, seq(0, 500, 50), las = 2)    
    points(dn$s_cort, dn$s_gluc, pch = 21, bg = alpha(s_col, 0.6), cex = 0.8)
    m <- lm(s_gluc ~ s_cort, data = dn)
    newx = seq(min(na.omit(dn$s_cort)), max(na.omit(dn$s_cort)), 0.05)
    conf_interval <- predict(m, newdata = data.frame(s_cort = newx), interval = "confidence",
                             level = 0.95)
    polygon(c(rev(newx), newx), c(rev(conf_interval[, 3]), 
                                  conf_interval[, 2]), col = alpha("gray40", 0.4), border = NA)
    lines(newx, conf_interval[, 1], lwd = 2)
    
    plot(1, 1, type = "n", yaxt = "n", xaxt = "n", xlab = expression(paste("Corticosterone (ng/", mu, "l)")),
         bty = "n", xlim = c(0, 40), ylim = c(80, 370), xaxs = "i", yaxs = "i", 
         ylab = "Blood Glucose (mg/dl)", main = "New York Nestling Post Dex")
    axis(1, seq(-20, 200, 10))    
    axis(2, seq(0, 500, 50), las = 2)    
    points(dn$d_cort, dn$d_gluc, pch = 21, bg = alpha(d_col, 0.6), cex = 0.8)
    m <- lm(d_gluc ~ d_cort, data = dn)
    newx = seq(min(na.omit(dn$d_cort)), max(na.omit(dn$d_cort)), 0.05)
    conf_interval <- predict(m, newdata = data.frame(d_cort = newx), interval = "confidence",
                             level = 0.95)
    polygon(c(rev(newx), newx), c(rev(conf_interval[, 3]), 
                                  conf_interval[, 2]), col = alpha("gray40", 0.4), border = NA)
    lines(newx, conf_interval[, 1], lwd = 2)
    
    plot(1, 1, type = "n", yaxt = "n", xaxt = "n", xlab = expression(paste("Corticosterone (ng/", mu, "l)")),
         bty = "n", xlim = c(0, 82), ylim = c(80, 420), xaxs = "i", yaxs = "i", 
         ylab = "Blood Glucose (mg/dl)", main = "New York Nestling Post ACTH")
    axis(1, seq(-20, 200, 10))    
    axis(2, seq(0, 500, 50), las = 2)    
    points(dn$a_cort, dn$a_gluc, pch = 21, bg = alpha(a_col, 0.6), cex = 0.8)
    m <- lm(a_gluc ~ a_cort, data = dn)
    newx = seq(min(na.omit(dn$a_cort)), max(na.omit(dn$a_cort)), 0.05)
    conf_interval <- predict(m, newdata = data.frame(a_cort = newx), interval = "confidence",
                             level = 0.95)
    polygon(c(rev(newx), newx), c(rev(conf_interval[, 3]), 
                                  conf_interval[, 2]), col = alpha("gray40", 0.4), border = NA)
    lines(newx, conf_interval[, 1], lwd = 2)

Individual delta relationships

Now plot the same things as above but for delta values of cort and glucose.

## Delta cort vs. delta glucose in adults and nestlings
    par(mfrow = c(3, 2))
    
    plot(1, 1, type = "n", yaxt = "n", xaxt = "n", xlab = expression(paste(Delta, "Corticosterone (ng/", mu, "l)")),
         bty = "n", xlim = c(-50, 100), ylim = c(-100, 200), xaxs = "i", yaxs = "i", 
         ylab = expression(paste(Delta, "Blood Glucose (mg/dl)")), main = "New York Adults Stress Response")
    axis(1, seq(-200, 200, 50))    
    axis(2, seq(-500, 500, 50), las = 2)    
    points(da2n$s_resp, da2n$gluc_resp, pch = 21, bg = alpha("gray40", 0.6), cex = 0.8)
    m <- lm(gluc_resp ~ s_resp, data = da2n)
    newx = seq(min(na.omit(da2n$s_resp)), max(na.omit(da2n$s_resp)), 0.05)
    conf_interval <- predict(m, newdata = data.frame(s_resp = newx), interval = "confidence",
                             level = 0.95)
    polygon(c(rev(newx), newx), c(rev(conf_interval[, 3]), 
                                  conf_interval[, 2]), col = alpha("gray40", 0.4), border = NA)
    lines(newx, conf_interval[, 1], lwd = 2)
    
    plot(1, 1, type = "n", yaxt = "n", xaxt = "n", xlab = expression(paste(Delta, "Corticosterone (ng/", mu, "l)")),
         bty = "n", xlim = c(-30, 85), ylim = c(-100, 200), xaxs = "i", yaxs = "i", 
         ylab = expression(paste(Delta, "Blood Glucose (mg/dl)")), main = "New York Nestling Stress Response")
    axis(1, seq(-500, 500, 20))    
    axis(2, seq(-500, 500, 50), las = 2)    
    points(dn$s_resp, dn$gluc_resp, pch = 21, bg = alpha("gray40", 0.6), cex = 0.8)
    m <- lm(gluc_resp ~ s_resp, data = dn)
    newx = seq(min(na.omit(dn$s_resp)), max(na.omit(dn$s_resp)), 0.05)
    conf_interval <- predict(m, newdata = data.frame(s_resp = newx), interval = "confidence",
                             level = 0.95)
    polygon(c(rev(newx), newx), c(rev(conf_interval[, 3]), 
                                  conf_interval[, 2]), col = alpha("gray40", 0.4), border = NA)
    lines(newx, conf_interval[, 1], lwd = 2)
    
    plot(1, 1, type = "n", yaxt = "n", xaxt = "n", xlab = expression(paste(Delta, "Corticosterone (ng/", mu, "l)")),
         bty = "n", xlim = c(-80, 5), ylim = c(-100, 100), xaxs = "i", yaxs = "i", 
         ylab = expression(paste(Delta, "Blood Glucose (mg/dl)")), main = "New York Adults Neg Feedback")
    axis(1, seq(-200, 200, 25))    
    axis(2, seq(-500, 500, 50), las = 2) 
    points(da2n$n_feed, da2n$gluc_feed, pch = 21, bg = alpha("gray40", 0.6), cex = 0.8)
    m <- lm(gluc_feed ~ n_feed, data = da2n)
    newx = seq(min(na.omit(da2n$n_feed)), max(na.omit(da2n$n_feed)), 0.05)
    conf_interval <- predict(m, newdata = data.frame(n_feed = newx), interval = "confidence",
                             level = 0.95)
    polygon(c(rev(newx), newx), c(rev(conf_interval[, 3]), 
                                  conf_interval[, 2]), col = alpha("gray40", 0.4), border = NA)
    lines(newx, conf_interval[, 1], lwd = 2)
    
    plot(1, 1, type = "n", yaxt = "n", xaxt = "n", xlab = expression(paste(Delta, "Corticosterone (ng/", mu, "l)")),
         bty = "n", xlim = c(-100, 25), ylim = c(-100, 104), xaxs = "i", yaxs = "i", 
         ylab = expression(paste(Delta, "Blood Glucose (mg/dl)")), main = "New York Nestling Neg Feedback")
    axis(1, seq(-200, 200, 50))    
    axis(2, seq(-200, 200, 50), las = 2)    
    points(dn$n_feed, dn$gluc_feed, pch = 21, bg = alpha("gray40", 0.6), cex = 0.8)
    m <- lm(gluc_feed ~ n_feed, data = dn)
    newx = seq(min(na.omit(dn$n_feed)), max(na.omit(dn$n_feed)), 0.05)
    conf_interval <- predict(m, newdata = data.frame(n_feed = newx), interval = "confidence",
                             level = 0.95)
    polygon(c(rev(newx), newx), c(rev(conf_interval[, 3]), 
                                  conf_interval[, 2]), col = alpha("gray40", 0.4), border = NA)
    lines(newx, conf_interval[, 1], lwd = 2)
    
    plot(1, 1, type = "n", yaxt = "n", xaxt = "n", xlab = expression(paste(Delta, "Corticosterone (ng/", mu, "l)")),
         bty = "n", xlim = c(-25, 100), ylim = c(-90, 90), xaxs = "i", yaxs = "i", 
         ylab = expression(paste(Delta, "Blood Glucose (mg/dl)")), main = "New York Adults ACTH Response")
    axis(1, seq(-200, 200, 25))    
    axis(2, seq(-500, 500, 25), las = 2)    
    points(da2n$a_inc, da2n$gluc_ainc, pch = 21, bg = alpha("gray40", 0.6), cex = 0.8)
    m <- lm(gluc_ainc ~ a_inc, data = da2n)
    newx = seq(min(na.omit(da2n$a_inc)), max(na.omit(da2n$a_inc)), 0.05)
    conf_interval <- predict(m, newdata = data.frame(a_inc = newx), interval = "confidence",
                             level = 0.95)
    polygon(c(rev(newx), newx), c(rev(conf_interval[, 3]), 
                                  conf_interval[, 2]), col = alpha("gray40", 0.4), border = NA)
    lines(newx, conf_interval[, 1], lwd = 2)
    
    plot(1, 1, type = "n", yaxt = "n", xaxt = "n", xlab = expression(paste(Delta, "Corticosterone (ng/", mu, "l)")),
         bty = "n", xlim = c(-50, 70), ylim = c(-210, 200), xaxs = "i", yaxs = "i", 
         ylab = expression(paste(Delta, "Blood Glucose (mg/dl)")), main = "New York Nestling ACTH Response")
    axis(1, seq(-200, 200, 25))    
    axis(2, seq(-500, 500, 50), las = 2)    
    points(dn$a_inc, dn$gluc_ainc, pch = 21, bg = alpha("gray40", 0.6), cex = 0.8)
    m <- lm(gluc_ainc ~ a_inc, data = dn)
    newx = seq(min(na.omit(dn$a_inc)), max(na.omit(dn$a_inc)), 0.05)
    conf_interval <- predict(m, newdata = data.frame(a_inc = newx), interval = "confidence",
                             level = 0.95)
    polygon(c(rev(newx), newx), c(rev(conf_interval[, 3]), 
                                  conf_interval[, 2]), col = alpha("gray40", 0.4), border = NA)
    lines(newx, conf_interval[, 1], lwd = 2)

There really doesn’t seem to be anything going on here. In a couple of cases, it seems like there might be a hint of a trend, but in those cases it is actually in the opposite direction predicted! Individuals with higher cort or a bigger increase in cort have lower glucose or a smaller increase in glucose. There certainly doesn’t seem to be ANY evidence in any of these plots that individuals with higher cort are also the ones that have higher glucose.

When you look through the tables below you can see that there are some cases where there are signficant relationships, but in most cases they are either in the opposite direction as predicted or else the effects are really really tiny (check the r squared) and are only significant because we have a very large sample size. There might be a couple of cases with more real looking effects…

Results: Tables

In all of these models I’ve scaled the values so that everything is in units of standard deviations. That makes it easier to interpret all of the coefficients.

Adult NY glucose

Fitting a series of models for NY adults with either glucose or change in glucose as the response variable and corticosterone, mass, and sex as predictors. Note that other predictors are possible but some are only available for a subset of individuals or years or are complicated by lots of different values in each year (like treatment). We could discuss the best way to deal with those or what to put in these models.

First fitting with absolute glucose levels:

## Models
    
    m1 <- lmer(scale(b_gluc) ~ scale(b_cort)*scale(mass) + sex + (1|band), data = da2n)
    m2 <- lmer(scale(s_gluc) ~ scale(s_cort)*scale(mass) + sex + (1|band), data = da2n)
    m3 <- lmer(scale(d_gluc) ~ scale(d_cort)*scale(mass) + sex + (1|band), data = da2n)
    m3.5 <- lm(scale(a_gluc) ~ scale(a_cort)*scale(mass), data = da2n)
    
    tab_model(m1, m2, m3, m3.5)
  scale(b_gluc) scale(s_gluc) scale(d_gluc) scale(a_gluc)
Predictors Estimates CI p Estimates CI p Estimates CI p Estimates CI p
(Intercept) 0.01 -0.07 – 0.09 0.822 0.04 -0.05 – 0.14 0.390 -0.20 -0.46 – 0.07 0.145 0.09 -0.68 – 0.86 0.817
scale(b_cort) 0.06 -0.01 – 0.13 0.113
sex: M -0.08 -0.31 – 0.16 0.538 -0.37 -0.62 – -0.11 0.005 0.30 -0.19 – 0.79 0.228
scale(b_cort):scale(mass) 0.04 -0.04 – 0.11 0.349
scale(mass) -0.10 -0.17 – -0.03 0.005 -0.08 -0.15 – 0.00 0.058 0.15 -0.08 – 0.38 0.214 0.09 -0.43 – 0.61 0.736
scale(s_cort) -0.06 -0.15 – 0.03 0.198
scale(s_cort):scale(mass) 0.09 0.01 – 0.17 0.031
scale(d_cort) -0.54 -0.92 – -0.16 0.005
scale(d_cort):scale(mass) -0.13 -0.44 – 0.18 0.422
scale(a_cort) -0.29 -1.32 – 0.73 0.568
scale(a_cort):scale(mass) -0.11 -0.75 – 0.54 0.739
Random Effects
σ2 0.97 0.81 0.69  
τ00 0.04 band 0.16 band 0.24 band  
ICC 0.04 0.17 0.26  
N 325 band 321 band 98 band  
Observations 761 577 102 45
Marginal R2 / Conditional R2 0.015 / 0.055 0.037 / 0.199 0.113 / 0.339 0.035 / -0.036

And then with delta glucose levels:

    m4 <- lmer(scale(gluc_resp) ~ scale(s_resp)*scale(mass) + sex + (1|band), data = da2n)
    m5 <- lmer(scale(gluc_feed) ~ scale(n_feed)*scale(mass) + sex + (1|band), data = da2n)
    m6 <- lm(scale(gluc_ainc) ~ scale(a_inc)*scale(mass), data = da2n)
    
    tab_model(m4, m5, m6)
  scale(gluc_resp) scale(gluc_feed) scale(gluc_ainc)
Predictors Estimates CI p Estimates CI p Estimates CI p
(Intercept) 0.05 -0.05 – 0.14 0.318 -0.17 -0.47 – 0.12 0.252 0.30 -0.35 – 0.95 0.357
scale(s_resp) 0.02 -0.07 – 0.11 0.679
scale(mass) 0.03 -0.04 – 0.11 0.393 0.09 -0.16 – 0.34 0.482 0.23 -0.23 – 0.69 0.320
sex: M -0.35 -0.61 – -0.09 0.009 0.51 0.01 – 1.01 0.048
scale(s_resp):scale(mass) 0.10 0.02 – 0.19 0.020
scale(n_feed) 0.02 -0.22 – 0.26 0.894
scale(n_feed):scale(mass) -0.05 -0.30 – 0.19 0.670
scale(a_inc) 0.42 -0.60 – 1.44 0.414
scale(a_inc):scale(mass) 0.40 -0.23 – 1.02 0.208
Random Effects
σ2 0.85 0.58  
τ00 0.13 band 0.43 band  
ICC 0.13 0.43  
N 318 band 98 band  
Observations 568 102 45
Marginal R2 / Conditional R2 0.027 / 0.156 0.038 / 0.448 0.085 / 0.018

Population comparison models

Fitting models similar to above but there are fewer here since it is only base and stress levels. I’ve include an interaction with state and cort in each model.

## Population comparison
    
    m1 <- lmer(scale(b_gluc) ~ scale(b_cort)*scale(mass) + sex + state*scale(b_cort) + 
                (1|band), data = da2)
    m2 <- lmer(scale(s_gluc) ~ scale(s_cort)*scale(mass) + sex + state*scale(s_cort) + 
                 (1|band), data = da2)
    
    m3 <- lmer(scale(gluc_resp) ~ scale(s_resp)*scale(mass) + sex + state*scale(s_resp) + 
                 (1|band), data = da2)
    
    tab_model(m1, m2, m3)
  scale(b_gluc) scale(s_gluc) scale(gluc_resp)
Predictors Estimates CI p Estimates CI p Estimates CI p
(Intercept) -0.10 -0.28 – 0.08 0.262 0.00 -0.20 – 0.20 0.984 0.11 -0.09 – 0.30 0.286
scale(b_cort) 0.12 -0.08 – 0.33 0.247
sex: M -0.12 -0.34 – 0.10 0.290 -0.40 -0.63 – -0.17 0.001 -0.40 -0.63 – -0.17 0.001
state: NY 0.28 0.08 – 0.48 0.006 0.13 -0.09 – 0.34 0.247 -0.06 -0.27 – 0.15 0.569
state: TN -0.02 -0.38 – 0.33 0.897 -0.05 -0.40 – 0.29 0.772 -0.25 -0.58 – 0.08 0.139
state: WY -0.28 -0.54 – -0.03 0.030 0.15 -0.14 – 0.45 0.316 0.35 0.06 – 0.64 0.019
scale(b_cort):scale(mass) -0.06 -0.14 – 0.02 0.128
scale(b_cort):stateNY 0.02 -0.20 – 0.24 0.849
scale(b_cort):stateTN -0.46 -0.97 – 0.06 0.082
scale(b_cort):stateWY -0.30 -0.57 – -0.03 0.031
scale(mass) -0.11 -0.19 – -0.04 0.004 -0.03 -0.11 – 0.05 0.435 0.07 -0.01 – 0.16 0.078
scale(s_cort) 0.01 -0.17 – 0.20 0.880
scale(s_cort):scale(mass) 0.00 -0.08 – 0.09 0.921
scale(s_cort):stateNY -0.13 -0.37 – 0.11 0.284
scale(s_cort):stateTN -0.04 -0.44 – 0.37 0.855
scale(s_cort):stateWY 0.03 -0.20 – 0.26 0.807
scale(s_resp) -0.03 -0.22 – 0.15 0.722
scale(s_resp):scale(mass) 0.04 -0.04 – 0.13 0.314
scale(s_resp):stateNY 0.06 -0.17 – 0.30 0.602
scale(s_resp):stateTN -0.06 -0.46 – 0.34 0.774
scale(s_resp):stateWY 0.11 -0.12 – 0.34 0.352
Random Effects
σ2 0.87 0.80 0.88
τ00 0.11 band 0.17 band 0.08 band
ICC 0.11 0.18 0.09
N 603 band 591 band 586 band
Observations 900 762 754
Marginal R2 / Conditional R2 0.054 / 0.159 0.031 / 0.202 0.047 / 0.129

Nestling models

Similar models to NY adults above, but here including nest identity as a random effect.

First for aboslute values:

## Nestling models
    dn$ubox <- paste(dn$site, dn$box, sep = "_")
    
    m1 <- lmer(scale(b_gluc) ~ scale(b_cort)*scale(mass) + (1|ubox), data = dn)
    m2 <- lmer(scale(s_gluc) ~ scale(s_cort)*scale(mass) + (1|ubox), data = dn)
    m3 <- lmer(scale(d_gluc) ~ scale(d_cort)*scale(mass) + (1|ubox), data = dn)
    m4 <- lmer(scale(a_gluc) ~ scale(a_cort)*scale(massd15) + (1|ubox), data = dn)
    
    tab_model(m1, m2, m3, m4)
  scale(b_gluc) scale(s_gluc) scale(d_gluc) scale(a_gluc)
Predictors Estimates CI p Estimates CI p Estimates CI p Estimates CI p
(Intercept) 0.09 -0.17 – 0.34 0.508 0.01 -0.23 – 0.25 0.954 0.01 -0.24 – 0.27 0.916 0.00 -0.17 – 0.18 0.969
scale(b_cort) -0.08 -0.24 – 0.08 0.320
scale(mass) -0.06 -0.20 – 0.09 0.456 0.11 -0.04 – 0.27 0.148 0.07 -0.08 – 0.22 0.381
scale(b_cort):scale(mass) 0.14 -0.00 – 0.29 0.054
scale(s_cort) -0.15 -0.29 – -0.01 0.039
scale(s_cort):scale(mass) 0.10 -0.05 – 0.24 0.188
scale(d_cort) -0.14 -0.28 – -0.00 0.043
scale(d_cort):scale(mass) 0.07 -0.09 – 0.24 0.386
scale(a_cort) 0.05 -0.11 – 0.22 0.516
scale(massd15) 0.33 0.17 – 0.49 <0.001
scale(a_cort):scale(massd15) -0.05 -0.19 – 0.09 0.483
Random Effects
σ2 0.40 0.40 0.33 0.81
τ00 0.60 ubox 0.52 ubox 0.61 ubox 0.10 ubox
ICC 0.60 0.57 0.65 0.11
N 43 ubox 41 ubox 41 ubox 40 ubox
Observations 182 176 175 157
Marginal R2 / Conditional R2 0.032 / 0.614 0.041 / 0.585 0.020 / 0.656 0.111 / 0.206

And then for delta values:

    m5 <- lmer(scale(gluc_resp) ~ scale(s_resp)*scale(mass) + (1|ubox), data = dn)
    m6 <- lmer(scale(gluc_feed) ~ scale(n_feed)*scale(mass) + (1|ubox), data = dn)
    m7 <- lmer(scale(gluc_ainc) ~ scale(a_inc)*scale(massd15) + (1|ubox), data = dn)
    
    tab_model(m5, m6, m7)
  scale(gluc_resp) scale(gluc_feed) scale(gluc_ainc)
Predictors Estimates CI p Estimates CI p Estimates CI p
(Intercept) -0.03 -0.21 – 0.15 0.748 -0.00 -0.19 – 0.18 0.972 0.02 -0.22 – 0.26 0.856
scale(s_resp) -0.16 -0.32 – -0.00 0.044
scale(mass) 0.17 0.00 – 0.34 0.047 0.05 -0.13 – 0.23 0.606
scale(s_resp):scale(mass) 0.05 -0.13 – 0.24 0.570
scale(n_feed) 0.04 -0.12 – 0.20 0.627
scale(n_feed):scale(mass) 0.04 -0.18 – 0.25 0.730
scale(a_inc) -0.06 -0.24 – 0.11 0.471
scale(massd15) 0.12 -0.05 – 0.29 0.173
scale(a_inc):scale(massd15) 0.04 -0.11 – 0.20 0.605
Random Effects
σ2 0.81 0.87 0.61
τ00 0.15 ubox 0.14 ubox 0.41 ubox
ICC 0.15 0.14 0.40
N 41 ubox 41 ubox 39 ubox
Observations 172 170 151
Marginal R2 / Conditional R2 0.049 / 0.196 0.005 / 0.143 0.024 / 0.416

Discussion

Fill in from above. Papers to discuss in context:

Newman-Lee et al. 2020. J Exp Biol. Investigating the relationship between corticosterone and glucose in a reptile.

Romero & Gormally. 2019. ICB. How truly conserved is the ‘well-conserved’ vertebrate stress response?

Lattin et al. 2015. Endocrinology. Are receptor concentrations correlated across tissues within indviduals? A case study examining glucocorticoid and mineralocorticoid receptor binding.

  • Obviously lots of others to add here, these ones just came to mind. The Newman-Lee paper can I think serve as a really good template for writing our paper since it is a very similar data set but in snakes. They even have the same parallel structure of correlational data + dex + acth and their results are pretty similar to ours.

  • There are some relationships with sex and mass that can be discussed and interpreted here even though they aren’t really the focus of the main questions posed above.

Conclusions

  1. At a population level, acute stress response is associated with an increase in glucose as predicted by our understanding of the direct consequences of corticosterone.

  2. Despite clear population level patterns, there is little evidence that individual variation in corticosterone regulation bears any relatinship with glucose regulation.

  3. We might have suspected that the relationship between cort and glucose would be context dependent, but even after checking across breeding stages, populations, and in adults vs. nestlings, there don’t seem to be any subsets where these relationships hold.

  4. The few cases where there is some evidence for a relationship (e.g. adult post-dex) it actually seems to be in the opposite direction as expected if intra-individual cort variation was directly causally responsible for glucose.

  5. Experimental manipulation of cort trajectory by dex or acth also does not seem to shift glucose (or maybe only a very little in adults).

  6. No simple causal relationship at the individual level, suggests need for more nuanced view of stress phenotype and multiple measures.