Field data and analysis
Field sampling
We captured adult juncos by mist net or potter trap at sites in Arizona (breeding season), Colorado (breeding), Illinois (non-breeding), Montana (breeding), New Mexico (non-breeding), New York (breeding), South Dakota (non-breeding), and Wyoming (breeding), spanning 16° in latitude and 37° in longitude (Fig. 1). This work was completed with approval from the U.S. Fish and Wildlife Service (MB84376B-1 to M.S.; MB01543B-0 to Z.A.C.; MB758442 to D.L.S.; MB06336A-4 to M.D.C.; MB757670-1 to David Winkler; and MB094297-0 to Christopher Witt,), the State of Arizona Game and Fish Department (SP590760 and SP707897 to D.L.S.), Colorado Parks and Wildlife (10TRb2030A15 to M.D.C.), the Illinois Department of Natural Resources (NH13.5667 to Z.A.C.), the Montana Department of Fish Wildlife and Parks (2016-013 and 2017-067-W to M.S.), the New Mexico Department of Game & Fish (3217 to Christopher Witt), the New York State Division of Fish, Wildlife, & Marine Resources (LCP 1477 to David Winkler), the State of South Dakota Department of Game, Fish, and Parks (06-03, 07-02, 08-03 to D.L.S.), Wyoming Game and Fish (754 to M.D.C.), and the Institutional Animal Care and Use Committees at Cornell University (2001-0051 to David Winkler), the University of Illinois (13385 to Z.A.C.), the University of Montana (010-16ZCDBS-020916 to Z.A.C.), the University of South Dakota (03-08-06-08B to D.L.S), and the University of Wyoming (A-3216-01 to M.D.C.). We classified individuals to taxonomic unit based on plumage (J. h. caniceps, J. h. hyemalis, J. h. mearnsi, J. h. oreganus group, and J. p. palliatus). The J. h. oreganus group encompasses seven subspecies (J. h. montanus, J. h. oreganus, J. h. pinosus, J. h. pontilis, J. h. shufeldti, J. h. thurberi, and J. h. townsendii) with similar plumages and overlapping nonbreeding ranges26,27. For this reason, we were unable to distinguish subspecies of nonbreeding individuals within this group, though all breeding individuals were collected from the J. h. montanus range.
Field metabolic assays
Birds were transported from the site of capture to a nearby laboratory (<50 mi away) for measurement where they were held indoors with ad libitum access to food and water. We assayed the Msum of each individual using open-flow respirometry. Most measurements were completed within 1 day of capture to avoid the effects of captivity on metabolic rates, though some (n = 15) were made within 2 days of capture (mean = 0.85 ± 0.47 day). Mb (in g) was quantified before each measurement began. Msum trials were conducted during daylight hours. A single individual was placed in a metabolic chamber in a dark, temperature-controlled environment. We pumped dry, cooled heliox gas (21% O2 and 79% He) through the animal’s chamber at a constant flow rate. We subsampled the outflow current, dried it (Drierite), scrubbed it of CO2 using ascarite, and dried it again before quantifying the O2 concentration using a FoxBox (Sable Systems). We recorded the output using Expedata software v.1.9.13. Trials were conducted using static cold exposure (−5 °C) for Colorado, Illinois, Montana, New Mexico, New York, and Wyoming birds and sliding cold exposure (starting at 0 °C and descending 3 °C every 20 min) for Arizona and South Dakota birds; however, both methods have been shown to produce similar estimates of Msum57. Trials lasted until O2 consumption plateaued or declined for several minutes. We also sampled a blank chamber before and after trials to account for potential fluctuations in baseline, ambient air.
We used custom R scripts (https://github.com/Mstager/batch_processing_Expedata_files) to correct for drift in baseline O2 content, calculate oxygen consumption (VO2)58, and quantify Msum (calculated as the highest instantaneous VO2 averaged over a 5-min period; ml O2 min−1). We discarded measures characterized by large drift in baseline O2 (owing to ambient temperature fluctuations affecting the FoxBox) or inconsistent flow rates resulting in a total sample size of n = 335 individuals. Following measurements, birds were subject to different fates: either released, exposed to acclimation experiments30, or immediately euthanized and deposited in museums (Source Data). Metabolic data from South Dakota have been previously published59.
Environmental data for field sampling sites
To account for an individual’s recent acclimatization history, we retrieved weather data associated with each collection site (rounded to the nearest hundredth of a degree latitude/longitude) from the DayMet dataset using the R package daymetr v.1.460. This dataset is composed of daily weather parameter estimates derived using interpolation and extrapolation from meteorological observations for a 1 km x 1 km gridded surface across North America61. We downloaded daily estimates of minimum temperature (Tmin; °C), maximum temperature (Tmax; °C), precipitation (prcp; mm/day), water vapor pressure (vp; Pa), shortwave radiation (srad; W/m2), and daylength (dayl; s/day) for the 14 days prior to each individual’s capture date. We additionally calculated daily temperature range (Td_range; °C) as Tmax – Tmin. We selected a conservative potential acclimatization window because, given their migratory nature, we do not know how long juncos were present at a site before sampling occurred, preventing us from incorporating more extensive windows (e.g., months). Prior work in J. hyemalis suggested a window of 7–14 days would be appropriate30,62. We, therefore, calculated averages across acclimatization windows varying from 7 to 14 days preceding capture (i.e., from 1–7 to 1–14) for each weather variable. We also retrieved elevation (elev; m) for each site using the package googleway v.2.7.363.
Analyses for field data
All analyses were conducted in R v.4.0.264. To determine whether junco Msum varied with environmental variation, we constructed eight linear models with Mb, taxon, and a single environmental variable (Tmin, Tmax, prcp, dayl, vp, srad, or Td_range averaged across the acclimatization window or elev) as main effects. We chose to use one variable at a time because several environmental variables were strongly correlated in our dataset (i.e., Tmin and Tmax: r = 0.96). We first standardized continuous predictor variables by centering and dividing by two standard deviations using the package arm v.1.11-265. As an indicator of differences in flexibility, we additionally included an interaction between the environmental variable and taxon to determine if taxa differed in their response to environmental cues. Thus models took the form Msum = Mb + taxon × environmental variable. We present the resulting standardized β to characterize the effect sizes of these relationships. We then used AIC values to evaluate differences in model fit among environmental variables and with that of a null model (including only Mb and taxon as predictors) in order to identify a single “best” environmental variable and acclimatization window explaining variation in Msum. To justify the inclusion of the interaction term, we compared models with and without the term using ∆AIC. We also used ∆AIC to test the inclusion of a variable for the season in which we sampled each individual (breeding or nonbreeding). The J. h. oreganus group served as the reference taxa in reported values.
Lastly, following the results of those analyses, which emphasized the importance of thermal variability, we tested whether junco thermogenic flexibility in the field corresponded with the heterogeneity of their thermal environment. As an index of thermal heterogeneity, for each individual, we downloaded the annual temperature range (˚C) for the site of capture from the WorldClim database66, composed of spatially interpolated climate data aggregated from 1970–2000, at a resolution of 2.5′ using the R package raster v.3.3-1367. Using these values, we calculated the mean annual temperature range for each taxon. We then correlated the reaction norms resulting from our best model (i.e., the taxon-specific slopes; ml O2/min/˚C) with the mean annual temperature range for each taxon using a linear regression.
Population genetic data
Sampling, sequencing, and SNP generation
For phylogeographic reconstruction, we obtained muscle tissue samples (n = 192) from museum specimens collected across the breeding distribution of all 21 Junco taxa identified by Miller26. We tried to maximize the environmental variance in our dataset (per ref. 68) by including 2–30 individuals per taxonomic unit sampled from 94 geographic localities representing the majority of US counties, Canadian provinces, and Mexican states for which tissue samples exist (Fig. 2a). We then employed restriction-site-associated DNA (RAD)-sequencing, which offers a reduced representation of the genome that can be mined for thousands of single nucleotide polymorphisms (SNPs) among individuals69,70.
We extracted whole genomic DNA from each sample using a Qiagen DNeasy Blood and Tissue Extraction Kit and prepared RAD-libraries according to ref. 71 (protocol v.2.3 downloaded from dryad: https://datadryad.org/stash/dataset/doi:10.5061/dryad.m2271pf1). Briefly, we digested whole genomic DNA with two restriction enzymes (EcoRI and Mse1), ligated adapter sequences with unique barcodes (8–10 bp) for each individual (provided in Source Data file), performed PCR amplification (see Table S2 for primers), and then performed automated size selection of 300–400 bp fragments (Sage Science Blue Pippen). We split paired geographic samples between the two libraries such that all taxa were represented in each library of 96, pooled individuals. Samples were additionally randomized on each plate during library preparation. Libraries were sequenced on separate flow-cell lanes of an Illumina HiSeq 4000 at UC Berkeley’s V.C. Genomics Sequencing Lab, yielding over 300 million, 100-nt single-end reads per lane.
We demultiplexed reads, retained reads with intact EcoRI cut sites, removed inline barcodes and adapters, and performed quality filtering (removed Phred score <10) using process_radtags in STACKS v.2.172. This resulted in final reads of μ = 92 bp in length and 2.05 million reads per individual. We removed three individuals in each lane that failed to sequence (<100,000 reads per individual, comprising five J. hyemalis and one J. insularis). We used bwa mem73 to align reads to the J. hyemalis genome37, which we downloaded from NCBI (Accession GCA_003829775.1). An average of 91% of reads mapped and mapping success did not differ among Junco species.
We executed the STACKS pipeline to call SNPs and exported one SNP per locus for loci present in all five Junco species in vcf format (populations: -p 5 –write-random-snp), resulting in 162,856 SNPs. We further filtered the dataset using vcftools v.0.1.1674. We removed sites with mean depth of coverage across all individuals <5 (–min-meanDP) and >50 (–max-meanDP), minimum minor allele count <3 (–mac), more than 50% missing data (–max-missing), and characterized by indels (–remove-indels). We then removed five J. hyemalis individuals with >60% missing data (–remove). Finally, we removed sites with >5% missing data (–max-missing) and filtered three sites that departed from Hardy–Weinberg equilibrium assessed by species (p < 0.001), resulting in 22,006 biallelic SNPs across 181 individuals. We exported the dataset in plink.raw format for downstream analysis.
Genotype-environment association analyses
To determine if environmental variation corresponded to allelic variation across Junco, we employed the SNP dataset in an RDA. RDA is a multivariate ordination technique that has been used to identify multiple candidate loci and several environmental predictors simultaneously35,36. Because RDA requires no missing data, we first imputed data for n = 68,010 missing sites (1.7% of total sites) using the most common genotype for the species at each site. We then quantified population genetic structure using a PCA of the imputed SNP data with the R package ade4 v.1.7-1575. We downloaded 19 spatially interpolated climatic variables aggregated from 1970–2000 corresponding to the geographic origin of each specimen from the WorldClim database66 at a resolution of 2.5′ using the R package raster v.3.3-1367. We centered and standardized climate variables by centering and dividing by two standard deviations65. We excluded highly correlated variables (r ≥ 0.70) resulting in the retention of eight variables: mean diurnal temperature range (BIO2, a measure of monthly temperature variation), temperature annual range (BIO7), mean temperature of the wettest quarter (BIO8), mean temperature of the driest quarter (BIO9), mean temperature of the warmest quarter (BIO10), precipitation of the driest month (BIO14), precipitation of the warmest quarter (BIO18), and precipitation of the coldest quarter (BIO19). We used these eight climatic variables as constraining variables in RDAs executed with the package vegan v.2.5-676. We performed a partial RDA conditioned on background population structure summarized as PC1 and PC2 from the PCA. We tested for, but did not find, multicollinearity among predictor variables (variance inflation factor <5 for all variables). We assessed the significance of the RDA model and of the constrained axes using an ANOVA-like permutation technique77 in vegan at p ≤ 0.05 (n = 999 and n = 299 permutations, respectively). We used variance partitioning78 to quantify the proportion of allelic variation explained by each climatic variable with the function varpart.
Acclimation experiments
Population sampling for acclimation experiments
We combined information gained from a literature search, eBird sightings, and expert opinion to identify five focal populations for phenotypic sampling that (1) were likely to be nonmigratory, (2) represented different morphological subspecies, and (3) maximized variation in annual temperature range within the United States. These populations include J. h. aikeni of the Black Hills, a coastal population of J. h. shufeldti, a highland population of J. h. dorsalis, a sky island population of J. p. palliatus, and a well-studied, urban population of J. h. thurberi79. However, it is possible that some of these populations exhibit seasonal, altitudinal migrations within their geographic area of residence, e.g., J. p. palliatus80.
We captured ≤25 adult individuals from each focal population. Capture periods differed for each population in order to increase the likelihood that targeted individuals were resident year-round, as well as due to time and permitting constraints. For instance, one partially migratory population (J. h. aikeni) with distinct morphological features was caught in the winter to increase the likelihood that the individuals used were nonmigratory. The other four populations, which bred in areas where other, morphologically similar juncos overwinter, were captured in the breeding season when other subspecies were not present. Specifically, J. h. shufeldti (n = 20) were captured 14–15 July 2018 in Coos and Douglas Counties, OR; J. p. palliatus (n = 24) were captured 27 July 2018 in Cochise County, AZ; J. h. dorsalis (n = 25) were captured 30–31 July 2018 in Coconino County, AZ; J. h. aikeni (n = 15) were captured 6–9 March 2019 in Lawrence County, SD; and J. h. thurberi (n = 20) were captured 22–26 July 2019 in San Diego County, CA. In spite of these differences, capture season did not influence acclimation ability (see below). This work was completed with approval from the US Fish and Wildlife Service (MB84376B-1 to M.S. and MB45239B-0 to T.J.G.), the State of Arizona Game and Fish Department (E19253811 to M.S.), the State of California Department of Fish and Wildlife (13971 to M.S.), the Oregon Department of Fish and Wildlife (108-18 to M.S.), the State of South Dakota Department of Game, Fish, and Parks (19-13 to M.S.), and the Institutional Animal Care and Use Committee at the University of Montana (030-18ZCDBS-052918 to Z.A.C.).
Acclimation treatments
Within days of capture, we ground-transported all birds to facilities at the University of Montana where birds were housed individually under common conditions (23 °C with 12 h dark: 12 h light) for ≥8 weeks (μ = 63 days, range = 57–71 days) to minimize the effects of prior acclimatization. We had previously determined that a period of 6 weeks in common conditions is sufficient to reduce variation in metabolic traits among J. hyemalis individuals (Fig. S2). Following this adjustment period, we assayed Msum (see below). We allowed birds ~24 h to recover and then randomly assigned individuals from each population into treatment groups and exposed them to either cold (3 °C) or control (23 °C) temperatures. Treatments lasted 21 days in duration. Constant 12 h dark: 12 h light days were maintained for the duration of the experiment, and food and water were supplied ad libitum. The diet consisted of a 2:1 ratio by weight of white millet and black oil sunflower seed, supplemented with ground dog food, live mealworms, and water containing vitamin drops (Wild Harvest D13123 Multi Drops). These experimental conditions were chosen based on previous work in J. h. hyemalis exposed to the same temperatures, which revealed substantial increases in Msum over the same duration21. At the end of treatments, we euthanized individuals using cervical dislocation.
Eight individuals died during the capture-transport and adjustment periods (one J. h. dorsalis, four J. h. shufeldti, one J. h. thurberi, and two J. p. palliatus). Additionally, one J. h. thurberi individual exhibited lethargy upon introduction to the cold treatment, died within the first 24 h of cold acclimation, and was removed from analyses. This resulted in a total sample size of n = 95 individuals.
Metabolic assays for acclimation experiments
We quantified Msum in a temperature-controlled cabinet using open-flow respirometry both before and after acclimation treatments as described above. We measured Mb immediately before each assay. During the post-acclimation Msum trial, we measured body temperature with a pit tag inserted into the cloaca (per ref. 30) to verify that individuals were hypothermic (body temperature ≤37 °C). Msum trials were conducted using static cold exposure at −5 °C for pre-acclimation measures and −17 °C for post-acclimation measures. Although this experimental temperature difference could have elicited higher post-acclimation Msum, we did not find variation in Msum between the two time points in Control birds, suggesting that the two procedures provided similar levels of cold challenge (paired t-test: n = 47, t = −1.24, df = 46, p = 0.22). Because trials occurred at various times throughout the day (08:00–19:00), we tested for, but did not find, a linear effect of the trial start time on Msum either before or after acclimation (before: R2 = 0.0, p = 0.54; after: R2 = 0.0, p = 0.94).
Climate data for acclimation populations
We reconstructed the annual thermal regime experienced by a population using interpolated monthly climate data downloaded from the WorldClim dataset66 with package raster v.3.3-13. We extracted all 19 variables at a resolution of 2.5′ for the site of capture. However, our analysis mainly focuses on BIO7, the temperature annual range variable (Trange).
Pairwise genetic distance of acclimated populations
To quantify genetic differentiation among acclimated populations, we extracted whole genomic DNA from the muscle tissue of acclimated individuals and prepared RAD-libraries as above with an additional PCR step (per protocol v.2.6)71. We pooled all 95 individuals into a single flow-cell lane of an Illumina HiSeq X for sequencing at Novogene Corporation Inc. yielding ~1 billion, 150-nt paired-end reads. We demultiplexed reads, removed inline barcodes and adapters, and performed quality filtering (removed Phred score <10) using process_radtags in STACKS. This resulted in final reads of μ = 147 bp in length and 7 million reads per individual. We removed six individuals that sequenced poorly (<200,000 reads/individual, comprising one J. h. dorsalis, one J. p. palliatus, one J. h. shufeldti, and three J. h. thurberi). We aligned reads to the J. h. carolinensis genome with bwa mem (51% of reads mapped). We then sorted, merged, and removed read duplicates with samtools v.1.381. We used the STACKS pipeline to call SNPs from loci present in all five populations and exported one SNP per locus in vcf format resulting in 282,929 SNPs (populations: -p 5 –write-random-snp).
We filtered the dataset and estimated patterns of genetic differentiation using vcftools. We removed sites with mean depth of coverage across all individuals <20 (–min-meanDP) and >300 (–max-meanDP), minor allele frequency <5% (–min-maf), and >10% missing data (–max-missing), as well as removed indels (–remove-indels), resulting in 1069 biallelic SNPs across 89 individuals. Finally, we estimated pairwise FST among the focal taxa by calculating the weighted Weir’s theta82.
Analyses for acclimation data
We first performed a partial mantel test to ascertain that environmental and genetic distances did not covary among our sampling sites. We calculated pairwise genetic distance as FST/(1 – FST)83. We estimated pairwise environmental differences among the five sites as the Euclidean distance for five WorldClim variables (after removing redundant variables at r ≥ 0.70 from the original 19 WorldClim variables). We simultaneously controlled for geographic distance, estimated as pairwise geodesic distance among sampling sites with the package geosphere v.1.5-1084. We then employed these matrices of pairwise genetic and environmental distances in a partial Mantel test (conditioned on geographic distance) with the package vegan v.2.5-676.
We quantified the effects of our acclimation treatments while simultaneously incorporating population demography using Markov Chain Monte Carlo generalized linear mixed models45,46 with the package MCMCglmm v.2.3285. This allowed us to include pairwise FST as a random effect in all analyses. Because individual measures were used in these models, the pairwise FST among two individuals corresponded to that of their respective subspecies. We standardized all continuous predictor variables by centering and dividing by two standard deviations65. In all cases, we used default priors and ran the model for 1,000,000 iterations with a burn-in of 10,000 and a thinning interval of 100. We examined the resulting trace plots to verify proper convergence.
We first verified that phenotypic differences did not exist among treatment groups before acclimations began. To do this, we constructed a model to explain variation in pre-acclimation Msum with temperature treatment as the main effect and pairwise FST as a random effect. We then repeated this procedure to test for differences among treatments in pre-acclimation Mb. We constructed similar models to test the effect of Trange on pre-acclimation measures of each Msum and Mb. We also tested for (but did not find) an effect of capture season (breeding or nonbreeding) on pre-acclimation Msum while including Mb as a covariate (season: posterior μ = 0.42, 95% CI = [−1.48, 2.30], p = 0.57; Mb: posterior μ = 0.45, 95% CI = [0.01, 0.88], p = 0.04).
To evaluate whether environmental variation corresponded with thermogenic flexibility, we constructed a model explaining variation in thermogenic flexibility (ΔMsum; post- minus pre-acclimation Msum) with pre-acclimation Mb, temperature treatment, Trange, and a treatment × Trange interaction term as main effects and pairwise FST as a random effect. We tested for (but did not find) an effect of capture season on ΔMsum while including Mb as a covariate (season: posterior μ = 0.09, 95% CI = [−0.42, 0.58], p = 0.72; Mb: posterior μ = 0.67, 95% CI = [0.33, 1.03], p = 6.0 × 10−4).
Finally, we calculated CV, (expressed as percentages)86 for ΔMsum across cold-acclimated individuals within each population using the package raster. To test whether populations from more variable environments expressed less variability in their flexible response, we then regressed the five CV values on Trange for each population using a linear regression. We also repeated this analysis for control-acclimated birds.
Source: Ecology - nature.com
