M1 – Protocol for selection experiment:
A starting population of 49 male and 49 female random-bred HS/IBG mice, representing 35 families, was obtained from the Institute of Behavioral Genetics at the University of Colorado, Boulder, CO, USA. No more than three mice were from any one family, and no more than two of either sex were from any one family. Mice were divided in a stratified fashion such that 12-13 mating pairs were randomly allocated to each of 4 blocks where no family was represented more than once (regardless of sex). We created four blocks from the starting population so that each selection treatment was replicated four times. Within each block, mice were allocated to three treatments: a randomly bred control group with no selection (CONT), directional selection for increased mass-independent MMR (high-MMR), and antagonistic selection between decreased mass-independent BMR and increased mass-independent MMR (ANTAG). Individual blocks were bred approximately 4 weeks apart to allow sufficient time to complete the physiological measurements. We measured BMR and MMR using the methods described in the main paper and detailed previously [1].
Prior to phenotypic analyses and selection on individuals for breeding the next generation, we screened the data from each generation for outliers using least squares multiple regression of BMR or MMR on body mass and other covariates (treadmill, age, observer (i.e. person who conducted the treadmill test), and which particular BMR chamber was used). Statistical analyses were performed using SAS, v. 9.1 (SAS Institute, Cary, NC, USA). If standardized residuals from these regressions were greater than absolute value 3.0, observations were omitted from the analyses. Most BMR outliers were positive suggesting that the mice had not been quiescent so that the measurements should not be considered to represent mice at their BMR. Similarly, most MMR outliers were negative, which likely resulted from mice that did not reach their MMR due to submaximal running effort. Significant covariates were included in the models as fixed effects.
Metabolic rates typically correlate strongly with mass. Hence, to minimize confounding metabolic effects with size effects, we selected on the mass-independent metabolism (i.e., on residuals from regressions of metabolism on body mass). Within each of the 4 blocks there were 3 treatments groups producing a total of 12 lines of mice per generation. Within each line, 13 pairings in general were made to produce offspring for the next generation, with the goal of obtaining 10 successful litters from each line. The breeding protocol for each selection treatment is described in the main paper.
Literature cited:
1.Wone B., Sears M.W., Labocha M.K., Donovan E.R., Hayes J.P. 2009 Genetic variances and covariances of aerobic metabolic rates in laboratory mice. Proc R Soc Lond, Ser B: Biol Sci276(1673), 3695-3704. (doi:10.1098/rspb.2009.0980).
M2 – Details of quantification of cytokines
We used Luminex technology (Luminex 200, Invitrogen Corporation, Carlsbad, CA 92008) to measure blood serum levels of four cytokines: tumor necrosis factor-alpha (TNF-), interleukin-6 (IL-6), interlukin-1 (IL-1), and granulocyte macrophage colony-stimulating factor (GM-CSF). Samples were processed using a Mouse Inflammatory Four-Plex Antibody Bead Kit (Invitrogen Corporation, Carlsbad, CA 92008, Catalog #LMC0003). We followed the procedure recommended by the kit manufacturers except that we used 25 L of blood serum and of assay standards instead of 50 L. We adjusted the results for the altered dilution. Luminex technology basically combines a traditional sandwich enzyme-linked immunsorbent assay (ELISA) with flow though cytometry, so that multiple signaling molecules can be quantified from a single sample [1, 2]. The blood samples had to be run on two trays due to space constraints. The values from the two trays were standardized for statistical analysis. Our average standard curves ranges were as follows: 8.24-5706.87 pgml-1for GM-CSF, 2.68-6402.06 pg ml-1for IL-1,5.14-11540.19 pg ml-1for IL-6, and 7.57-1878.29 pg ml-1for TNF-.
Literature cited:
1.Dasso J., Lee J., Bach H., Mage R.G. 2002 A comparison of ELISA and flow micro sphere-based assays for quantification of immunoglobulins. J Immunol Methods263(1-2), 23-33.
2.Vignali D.A.A. 2000 Multiplexed particle-based flow cytometric assays. J Immunol Methods243(1-2), 243-255
M3 - ELISA details:
A single ELISA was run on the sample from each individual following the protocol described by Baze et al. [1]. Because mice are small, we could only take a small serum sample and were thus limited to a single ELISA for each individual. Thus, we chose a secondary antibody that allowed us to quantifyIgA, IgG, and IgM in a single ELISA over the more traditional approach of quantifying IgG and IgM separately. This approach would give us the most comprehensive picture of the antibody response given the limiting amount of serum available.
Briefly, microtiter plates (96-well MaxiSorp polystyrene plates, Fisher Scientific, Fairlawn, NJ) were coated with KLH (50 L of 10 g ml-1, Sigma-Aldrich catalog # H7017), sealed, and incubated overnight at 4 C. The plates were washed 4 times with PBS containing 0.05% Tween (PBS-T) and wells were blocked with 5% nonfat dry milk in PBS. Sealed plates were incubated for 2 h at 4 C and washed 4 times with PBS-T. Starting with a thawed serum sample of 14 L, samples were diluted to 1:10, 1:100, 1:1000, and 1:1000 with PBS-T and then 50 L of each dilution was added in duplicate to the wells of the antigen coated plate. Plates were sealed and incubated at 4C for 45 min, and then washed (4 times) with PBS-T. 50 L peroxidase-conjugated, goat anti-mouse immunoglobulins (polyvalentIgA, IgG, IgM, 1:10,000 dilution, Sigma-Aldrich catalog # A0412) were added to the wells, then plates were sealed and incubated at 4 C for 45 min. Plates were washed with PBS-T and 50 L of freshly prepared TMB Microwell Peroxidase Substrate (catalog # 50-76-03,KPL, Gaithersburg, MD) was added to each well. Plates were protected from light during the enzyme reaction. After 30 min, the reactions were stopped by adding 50 L of 0.5 M phosphoric acid to each well. The optical density for each well was determined at 450 nm using a Spectra-Max micro-ELISA reader (Molecular Devices, Hercules, CA) [1]. The response variable for our analysis was the average optical density of the serum samples diluted 1:1000. We subtracted a negative control from each optical density measurement to account for background noise caused by unavoidable nonspecific binding. We also standardized across trays using a control, so all optical densities were comparable. Optical density of the serum sample diluted 1:1000 was chosen because it fell within the linear part of the dilution curve.
Literature cited:
1.Baze M.M., Hunter K., Hayes J.P. 2011 Chronic hypoxia stimulates an enhanced response to immune challenge without evidence of an energetic tradeoff. J Exp Biol214, 3255-3268. (doi:10.1242/jeb.054544).
M4 - Details of statistical analyses
General procedure: All analyses were conducted following the “protocol” suggested in Zuur et al. [1] for finding the optimal combination of random components (possibly including variance covariates and random effects) and fixed effects (see section electronic supplement M5 for annotated example). For each response variable, restricted maximum likelihood (REML) models were first fit with the “beyond optimal” combination of fixed effects (i.e., the model that contained all reasonable fixed effects), and the optimal random effect structure was determined by likelihood ratio tests (LRT) of nested fitted models and visual examination of diagnostic plots of standardised residuals. After finding the optimal random component, nested models fit with maximum likelihood(ML) were compared with likelihood ratio tests (LRT) and appropriate fixed effects selected by sequentially dropping model terms. Significance of pairwise contrasts between factor levels was determined post-hoc with Tukey all-pair comparisons of REML models.
SHAM and immune treatment (LPS or KLH) difference: The random structure for LPS mice varied by response variable; TNF- was modeled with a random effect of block and variance covariates of LPS treatment and line, IL-6 was modeled with a random effect of block and variance covariate of line, and GMC-SF was modeled with a random effect of block and variance covariates of LPS treatment and block. The random structure for KLH mice included variance covariates of KLH treatment and selection treatment. Fixed effects were selection treatment, immune treatment (LPS or KLH), and the interaction of selection treatment with immune treatment.
Selection treatment level (LPS or KLH mice only): The random component for TNF- and GMC-SF models included a variance covariate of selection treatment, and no specific random component was deemed necessary for the IL-6 models. The KLH models included a random effect of line and a variance covariate of selection treatment. Fixed effects in the full models were mass, BMR, MMR, and selection treatment.
BMR and MMR (LPS or KLH mice only): Because the mice we studied were subsets of mice from the larger selection experiment, we also tested for differences in BMR and MMR in each of these subsets of mice. That is, for mice used to study the inflammatory response and the humoral response, we tested whether selection treatment affected mass-independent BMR and mass-independent MMR. Body mass at the time of metabolic measurement was a covariate in these analyses to make the metabolic rate variables mass-corrected. For LPS mice, BMR models included a variance covariate of line, whereas MMR models included a random effect of line. For KLH mice, BMR models included a variance covariate of line, and MMR models did not require specified random components. These analyses enabled us to determine if BMR and MMR in the subsets of mice we randomly chose for the experiment matched the pattern exhibited by the whole population of mice (i.e. generation 8).
Literature Cited:
1.Zuur A.F., Ieno E.N., Walker N.J., Saveliev A.A., Smith G.M. 2009 Mixed effects models and extensions in ecology with R. New York, NY, Springer Science+Business Media, 574 p.
M5 - Example annotated R code of model selection process
Below is an annotated example of model selection procedure we used to develop our statistical models. The example details the statistics for antibody concentration from the adaptive immune experiment. This code demonstrates how we determined that the mice injected with the saline sham (SHAM) differed from the mice inject with keyhole limpet haemocyanin (KLH), that selection treatment was not an important predictor of antibody concentrations in mice injected with KLH, and that maximal metabolic rate (MMR) and mass are important predictors of antibody concentration in mice injected with KLH. A text description of the modeling procedure for the difference in antibody concentrationsbetween SHAM and KLH mice is presented in the methods of the main text. For a text summary of the model selection procedure and a description of the final “best fit” model for all response variables please see electronic supplement M4.
# all analyses performed in R 2.14.1, with packages
# nlme, multcomp, and contrasts
#
# R Development Core Team (2011). R: A language and environment
# for statistical computing. R Foundation for Statistical
# Computing, Vienna, Austria. ISBN 3-900051-07-0,
# URL
# loading in relevant libraries
library(nlme)
library(multcomp)
library(contrasts)
# loading in data from .csv file
klh.all <- read.csv("KLH.csv")
# making sure all factor variables are treated as such
# renaming SELEC as Selec for consistency with other data sets
#LINETYPE = selection treatment (coded as a number), BLK = selection block,
# Line = line (1-12), Selec = selection treatment (coded as text).
klh.all$LINETYPE <- factor(klh.all$LINETYPE)
klh.all$BLK <- factor(klh.all$BLK)
klh.all$Line <- factor(klh.all$LINE)
klh.all$Selec <- factor(klh.all$SELEC)
klh$Selec <- klh$SELEC
# removing records with missing response variable, OD4
# (OD4 is the response variable for antibody concentrations associated with KLH injections
klh <- klh.all[!is.na(klh.all[,"OD4"]),]
# visual exploration of data
hist(klh$OD4, breaks=13)
boxplot(OD4 ~ Line, data=klh)
boxplot(OD4 ~ Selec, data=klh)
boxplot(OD4 ~ KLHtreat, data=klh)
#KLHtreat = the adaptive immune response group (SHAM or KLH)
boxplot(OD4 ~ BLK, data=klh)
##################################################
# first, looking for treatment (injection) effect to determine if the KLH injections
# actually elevated antibody concentrations above those seen in the SHAM injected mice
##################################################
# fitting first model with no specified random component
# but all intended fixed effects (beyond optimal model)
# this will be our base comparison for LRT tests
klh.gls1 <- gls(OD4 ~ KLHtreat * Selec, data=klh)
# extracting standardized residuals to check for
# heterogeneity in variance, other model diagnostics
E <- resid(klh.gls1, type="normalized")
Fit <- fitted(klh.gls1)
plot(x=Fit, y=E, xlab="Fitted Values", ylab="Stand. Resid", main="OD4")
# very strong heterogeneity in variance as evidenced by fan shape
# to plot
boxplot(E ~ Line, data= klh, xlab="Lines", ylab="Stand. Resid")
# only minor heterogeneity in Line
boxplot(E ~ KLHtreat, data= klh, xlab="KLHtreat", ylab="Stand. Resid")
# lots of heterogeneity in treatment
boxplot(E ~ Selec, data= klh, xlab="Selection", ylab="Stand. Resid")
# some heterogeneity also in Selec
boxplot(E ~ Line + KLHtreat, data= klh, xlab="Lines", ylab="Stand. Resid")
boxplot(E ~ Selec + KLHtreat, data= klh, xlab="Selection", ylab="Stand. Resid")
# this combination captures a lot of heterogeneity
boxplot(E ~ BLK, data= klh, xlab="Block", ylab="Stand. Resid")
# a bit of variance here too
boxplot(E ~ BLK + Selec, data= klh, xlab="Block", ylab="Stand. Resid")
boxplot(E ~ BLK + KLHtreat, data= klh, xlab="Block", ylab="Stand. Resid")
# checking for normality of residuals
qqnorm(klh.gls1, main="KLH", abline=c(0,1))
# not at all normal
# starting to explore various random component structures
# “weights” term here allows variance to differ between specified
# strata; here, the strata are KLH treatment group
klh.gls2 <- gls(OD4 ~ KLHtreat * Selec, data= klh, weights=varIdent(form = ~1 | KLHtreat))
# anova comparison of nested models, likelihood ratio test (LRT) reported
anova(klh.gls1, klh.gls2)
# much better by AIC and LRT, dAIC of hundreds
# making and checking diagnostic plots of standardized residuals
E2 <- resid(klh.gls2, type="normalized")
Fit <- fitted(klh.gls2)
plot(x=Fit, y=E2, xlab="Fitted Values", ylab="Stand. Resid")
#other than a few outliers, this looks pretty good
boxplot(E2 ~ Line, data= klh, xlab="Lines", ylab="Stand. Resid")
boxplot(E2 ~ KLHtreat, data= klh, xlab="LPStreat", ylab="Stand. Resid")
# still minor heterogeneity in line
boxplot(E2 ~ Selec, data= klh, xlab ="Lines", ylab="Stand. Resid")
qqnorm(klh.gls2) # pretty good
# process continued for reasonable combinations
# of variance covariates and random effects
# models evaluated by LRT and diagnostic plots
# not all models shown
klh.gls4 <- gls(OD4 ~ KLHtreat * Selec, data= klh, weights=varIdent(form = ~1 | KLHtreat * Line))
anova(klh.gls4, klh.gls1)
# final model selected for random effect component
# on to assessing the fixed effects
# specifically checking whether KLH treatment and its
# interaction are significant
# first, refit full model with maximum likelihood (ML).
# Must use ML to compare models with different fixed effects rather than
# using restricted maximum likelihood (REML).
klh.gls4.ml <- gls(OD4 ~ KLHtreat * Selec, method="ML", data= klh, weights=varIdent(form = ~1 | KLHtreat * Line))
# next, fit a reduced model with same random component, in ML
klh.gls4.selec.ml <- gls(OD4 ~ Selec, method="ML", data= klh, weights=varIdent(form = ~1 | KLHtreat * Line))
# finally, compare with anova for LRT
anova(klh.gls4.ml, klh.gls4.selec.ml)
# comparison is significant, indicating dropped model terms
# were significant and substantially helped model fit
# Model df AIC BIC logLik Test L.Ratio p-value
#klh.gls4.ml 1 24 -232.7969 -178.4926 140.3984
#klh.gls4.selec.ml 2 21 -185.1099 -137.5937 113.5550 1 vs 2 53.68693 <.0001
#########################################
# move to selection level for KLH
##########################################
# removing those records of SHAM treated mice
klh.only <- klh[klh$KLHtreat=="KLH",]
# removing records without MMR measurements
klh.only <- klh.only[!is.na(klh.only[,"MMR"]),]
# exploratory visual analysis of data
boxplot(OD4 ~ Line, data=klh.only)
boxplot(OD4 ~ Selec, data=klh.only)
boxplot(OD4 ~ BLK, data=klh.only)
plot(klh.only$BMR, klh.only$OD4)
plot(klh.only$MMR, klh.only$OD4)
plot(klh.only$KLHmass4, klh.only$OD4)
# appears to be patterns in the factors, but not continuous
# (other than mass)
# again, first fitting a model with no specified random component
klh.gls1.trt <- gls(OD4 ~ Selec + KLHmass4 + MMR + BMR, data=klh.only)
# checking the residuals for heterogeneity and normality
E <- resid(klh.gls1.trt, type="normalized")
Fit <- fitted(klh.gls1.trt)
plot(x=Fit, y=E, xlab="Fitted Values", ylab="Stand. Resid", main="KLH")
# two possible outliers, not much of a pattern
boxplot(E ~ Line, data= klh.only, xlab="Lines", ylab="Stand. Resid")
boxplot(E ~ Selec, data= klh.only, xlab="Selection", ylab="Stand. Resid")
boxplot(E ~ BLK, data= klh.only, xlab="Block", ylab="Stand. Resid")
qqnorm(klh.gls1.trt, abline=c(0,1)) # really not too bad
# after exploring plausible random component structures
# settled on the structure here as optimal
# contains a random effect of Line and variance covariate
# of selection treatment
klh.gls8.trt <- lme(OD4 ~ Selec + KLHmass4 + MMR + BMR, random = ~ 1 | Line, data= klh.only, weights=varIdent(form = ~1 | Selec))
anova(klh.gls8.trt, klh.gls1.trt) # sig LRT, AIC=26.62765
E2 <- resid(klh.gls8.trt, type="normalized")
Fit <- fitted(klh.gls8.trt)
plot(x=Fit, y=E2, xlab="Fitted Values", ylab="Stand. Resid")
# better than gls1
boxplot(E2 ~ Line, data= klh.only, xlab="Lines", ylab="Stand. Resid")
boxplot(E2 ~ Selec, data= klh.only, xlab ="Lines", ylab="Stand. Resid")
boxplot(E2 ~ BLK, data= klh.only, xlab="Block", ylab="Stand. Resid")
qqnorm(klh.gls8.trt, abline=c(0,1)) # good
# refitting in ML to examine fixed effects
klh.gls8.trt.ml <- lme(OD4 ~ Selec + KLHmass4 + MMR + BMR, random = ~ 1 | Line, weights=varIdent(form = ~1 | Selec), method="ML", data= klh.only)