I conducted all analyses in R (R Development Core Team 2021) and scaled and centered all variables. Code, data summaries and full results are supplied in a supplementary R Markdown script.

Aim 1 (evaluate the ubiquity of captive wing shape phenotypes) – I used the lengths of L_W, L_S, L_P, and ΔQ (P10-P3) the response variable in Bayesian logistic regression (family: ‘gaussian’) implemented in the package MCMCglmm (Hadfield 2010). I included a three-way interaction between the measured feather ID, species ID and provenance (captive/wild) as the fixed effect, and the specimen ID as the random effect. I used inverse-Wishart priors for the random effect and residual variance (V = 1, ν = 0.002). I ran the model for 100,000 iterations and set burn-in to 1,000 and used 100 for thinning for a total posterior sample of 990. All chains were checked for proper mixing, and I checked for auto correlation using the command ‘autocorr’ with 0.1 as a target threshold. I used emmeans (Lenth 2018) to estimate Z values of pairwise captive-wild contrasts for visualization with ggplot2 (Wickham 2016).

Aim 2 – To evaluate the relationship between survival of the first migration and individual phenotype, I used the binomial juvenile survival outcome as the response variable in Bayesian logistic regression (family: ‘categorical’), implemented in the package MCMCglmm (Hadfield 2010). I fitted ΔQ (P10-P5), L_T, and the index of body condition as fixed effects and release year as a random effect. I specified priors for the fixed effects using the ‘gelman.prior’ command (v=1, nu=0.02) and fixed residual variance at 1. I ran the model for 100,000 iterations and set burn-in to 1,000 and used 100 for thinning for a total posterior sample of 990. All chains were checked for proper mixing, and I checked for auto correlation using the command ‘autocorr’ with 0.1 as a target threshold.