An understanding of climate change related effects on populations is essential for assessing vertebrate life history and conservation. However, these potential relationships are difficult to measure because of the relatively long generation time of many vertebrates and the sometimes subtle and long-term changes in climate. In addition, long-term data streams are not available or are not of adequate quality for many vertebrate taxa. Therefore, our efforts to examine effects of climate change on populations should focus on systems with the best ‘design’. First, climate change induced effects on vertebrates are most likely to be manifest through range shifts, particularly for migratory species. Recently observed shifts in species ranges have been documented on a global scale and are most likely to occur in breeding ranges because evidence indicates that spring has been advancing by 2-5 days per decade (Parmesan & Yohe 2003, Root et al. 2003, Parmesan 2007). Second, climate change effects are spatially heterogeneous with relatively accelerated rates of change in northern environments (Hinzman et al. 2005, Wendler & Shulski 2008) resulting in a variety of landscape-level changes such as reduction in glacial extent and an increase in permafrost temperatures in many areas (Serreze et al. 2000). Subarctic Alaska is a ‘hotspot’ of climate warming with the highest increases in seasonal temperatures over the past 50 years occurring in winter (2.2°C increase) and spring (Stafford et al. 2000).

Investigating shifts in the range of migratory species can be particularly challenging because few monitoring schemes are conducted at spatial scales that cover a majority of a species' range, particularly for birds (Parmesan 2006). However, continued changes in climate will likely have many impacts on Arctic and subarctic bird populations (Seavy et al. 2008), and there is broad evidence that avian species' ranges are shifting northwards (La Sorte & Thompson 2007), most likely due to a longer growing season and longer ice-free periods for birds that use wetlands. However, for wetland-breeding birds some of the potential benefits of climate change may be offset through the loss or shrinking of wetlands across the landscape (e.g. Klein et al. 2005, Riordan et al. 2006). Reductions in size and numbers of wetlands would reduce overall available habitat, and changes in plant and invertebrate communities (Smol et al. 2005) could alter food availability.

We used a long-term data set on a northern breeding bird that relies on wetlands for nesting, the trumpeter swan Cygnus buccinator, to complete one of the first detailed examinations of the potential effects of climate change on the breeding range of an avian species. Because the trumpeter swan primarily uses boreal forest wetlands for reproduction and because this species requires a relatively long breeding period to successfully fledge young (145-150 days; Hansen et al. 1971), changes in the availability of wetlands induced by climate change could have a substantial impact on the entire population. This species has recently recovered from near extinction in the early 1900s with populations in Alaska increasing at a rate of nearly 6% annually since 1968 (Caithamer 2001, Conant et al. 2002, Schmidt et al. 2009a). Recent research has also suggested that the population may be expanding at a higher rate in more recent years in more northern areas (Schmidt et al. 2009a), where the breeding range now extends beyond the historic northern limit (Hansen et al. 1971). This expansion into new areas has resulted in recent observations of an overlap of the breeding distribution with that of the smaller tundra swan Cygnus columbianus in some areas. Historically, the length of season required for successful breeding (ca 30-45 days less for tundra swans) may have functionally separated these two species in space, but a warming climate may induce new forms of competition if range overlap increases. Recolonization of historic habitats by similar species, the whooper swan Cygnus cygnus and mute swan Cygnus olor, have not been found to displace other waterfowl species (Poysa & Sorjonen 2000, Gayet et al. 2011). However, the introduction of mute swans to areas where they have not occurred historically can have profound impacts on other waterfowl populations and wetland habitats (Petrie & Francis 2003). Because trumpeter swans have not historically occupied tundra habitats in large numbers, range increases may negatively affect populations of tundra swans and other species using tundra wetlands.

Hypotheses and predictions

Our primary hypothesis was that recolonization patterns in trumpeter swans, as well as overall range expansions, may be related to recent increases in the ice-free period observed across the northern hemisphere (Magnuson et al. 2000) and in Alaska in particular (Sagarin & Micheli 2001). Hansen et al. (1971) originally identified a theoretical northern limit to the trumpeter swan breeding range based on the annual number of ice-free days in a given season. Based on this, we hypothesized that habitat occupancy patterns were related to breeding season length (possibly in a quadratic fashion), and a positive relationship between season length, as indexed by the number of days above freezing, and the probability of habitat occupancy would provide evidence that trumpeter swans are colonizing areas in response to breeding season length. To address this season-length hypothesis, we examined habitat occupancy throughout Alaska during 1968-2005 relative to the average number of days above freezing (and days above freezing²) during the previous five years.

Based on previous research utilizing historic Alaskan trumpeter swan survey data (Schmidt et al. 2009a,b), we also hypothesized that occupancy would generally increase due to recolonization of historic breeding areas as abundance increased, and we needed to control for this natural recolonization process to isolate any effects of season length. Although the full extent of the historic range is unknown, we expected that annual trend and/or the latitude of the survey unit would act as surrogates for the recolonization process. We also expected that the recolonization process may be non-linear (latitude²) reflecting the tendency of a recovering population to expand outward from a core range, in this case, middle latitudes.

Material and methods

Statewide surveys for trumpeter swans were first conducted in Alaska in 1968, and in 1975, 1980, 1985, 1990, 1995, 2000 and 2005 surveys covering all known breeding areas were conducted once during August and/or September (Conant et al. 1996, Caithamer 2001, Schmidt et al. 2009a). Individual survey units were 1:63,360 quadrangle maps (units) where trumpeter swans were suspected to occur and were surveyed from a single engine fixed-wing aircraft with the intent of counting every individual within each unit. Survey units were considered to be occupied in a given year if ≥1 swan was sighted, and additional units were added through time (e.g. N = 176 in 1968 to N = 780 in 2005; Fig. 1) in response to apparent range expansion statewide (Schmidt et al. 2009a). We were unable to estimate detection probability from the data available, but previous work concluded that detection probability was high and likely constant for individual swans (Schmidt et al. 2009b), minimizing the potential effects on occupancy modeling for our large sampling units. Tundra swans, which could not be separately identified from trumpeter swans, were present in a limited number of units at the western edge of the observed range of trumpeter swans, but we did not expect the limited scope of this situation and the extensive overlap of the two species within individual survey units occupied by tundra swans to influence the results for trumpeter swans. We based the identification of units to be added to the survey on annual statewide waterfowl surveys during the years between swan surveys, which provide good coverage of waterfowl habitat in the state. Surveys covered the entire known breeding range in Alaska each year (see Fig. 1).

We acquired temperature information for our survey area from weather stations (N = 14) around the state that had collected data for a majority of the period from 1964 to 2005 (see Fig. 1). We then calculated the average number of days with a minimum temperature > 0°C during the potential breeding season (April-October) of the survey year and the four years prior to the survey (Fig. 2). We assumed that this 5-year average would be representative of the conditions that may have affected occupancy during each corresponding 5-year survey. We assigned units a number of days above freezing for each survey period from the weather station nearest to the center point of each unit, although in some cases, topography (i.e. mountain ranges) was taken into account when selecting the appropriate station. Most groupings of survey units contained a weather station near the center of the group to provide a representative number of days above freezing (see Fig. 1).

We used Bayesian hierarchical occupancy models (Royle & Dorazio 2008) to investigate the effect of the number of days above freezing on the probability of occupancy for each survey unit during 1968-2005. Because our data were complicated by changes in survey extent, we used Markov Chain Monte Carlo (MCMC) methods to fit models describing variation in trumpeter swan occupancy through time relative to year, latitude and days above freezing. We standardized or scaled parameters to have a mean = 0 and SD = 1 to facilitate convergence. These methods allowed us to include all survey units (N = 780) in the analysis each year regardless of the year that they were added to the survey, and they also allowed us to investigate models containing random effects (Schmidt et al. 2009a). In particular, if a unit was not surveyed in a given year, we treated it as a missing value and it was estimated along with all of the other parameters. This permitted unequal numbers of samples from many of the units to be included in the analysis.

When constructing our model set, we assumed that there was a trend in occupancy through time because of the apparent increase in the species' range and the increasing trend in the overall population (Schmidt et al. 2009a). This assumption helped us to control for general population increases and identify any influence of changes in average temperature on occupancy. We suspected that differences in recolonization rates at various latitudes might help explain changes in occupancy probability that were not explained by the season-length covariates alone. We used the Deviance Information Criterion (DIC) to select among competing models containing combinations of these variables (Spiegelhalter et al. 2002) and then added a normally distributed random effect of survey year as an intercept adjustment to the top fixed-effects model to help account for any remaining annual heterogeneity in the data. The resulting model set contained seven models. DIC model selection procedures are analogous to those of Akaike's Information Criterion, and we considered models within 10 units of one another to be competitive.

We used R 2.4.0 (R Development Core Team 2006) and WinBUGS 1.4.1 (Spiegelhalter et al. 2004) to fit our set of candidate models. We generated three independent Markov chains up to 200,000 iterations in length, and convergence properties were assessed using the Gelman-Rubin diagnostic (Gelman & Rubin 1992, Brooks & Gelman 1998; Rhat < 1.1 = convergence) and a visual inspection of the chains. Initial iterations were discarded as burn-in, and estimates of the posterior distributions of the parameters were calculated from the remaining iterations. The number of iterations necessary for convergence varied depending on model complexity. We estimated goodness-of-fit (P) in a manner similar to that used by Link & Sauer (2002) and Schmidt et al. (2009a) which compares the fit of the model to the fit of a replicate data set generated during each update. By comparing the observed and replicate data, we were able to assess whether the data could have come from the proposed model. Good fit occurs when P = 0.5, and major failures of the model occur when P < 0.01 or P > 0.99 (Gelman et al. 2004:175). Parameter estimates indicate means and 95% Bayesian credible intervals.

Results

The 5-year average number of days above freezing observed at each weather station generally decreased as latitude increased, although variation was high (see Figs. 1 and 2). Model selection results indicated support for a model containing quadratic relationships of both latitude and the number of days above freezing. Models containing both of these variables performed much better (DIC = 4,380) than reduced models containing the latitude (DIC = 4,402) or days above freezing (DIC = 4,460) covariates alone, and the inclusion of a random year effect improved the final model (DIC = 4,367; Table 1). Model fit was determined to be satisfactory based on the Bayesian P-value (P = 0.44 for the simplest model).

Individual logit-scale parameter estimates from the best approximating model in the model set were fairly imprecise (Table 2), indicating uncertainty in their magnitude and direction, likely due to correlation between the latitude and days above freezing variables. Despite uncertainty in individual estimates and annual fluctuations (Fig. 3), predicted probability of occupancy increased from approximately 0.3 to 0.7 at both low and high latitudes and from approximately 0.7 to 0.9 at mid-latitudes over the observed range of the number of days above freezing (Fig. 4). Estimates for the general trend through time tended to be positive (0.238; 95% CI: −0.018 - 0.476). Estimates from models including only latitude or days above freezing covariates were unambiguous in their directions. Latitude had negative linear (-0.279; 95% CI: −0.361 - −0.198) and quadratic (-0.393; 95% CI: −0.473 - −0.314) effects, and days above freezing had a positive linear (11.726; 95% CI: 7.162 - 15.94) and negative quadratic (-3.767; 95% CI: −5.246 - −2.170) effects on occupancy on the logit scale in these two simpler models, respectively.

Discussion

Strong support for models with a relationship between habitat occupancy and both the average number of days above freezing and recolonization parameters supported our hypotheses that observed patterns of occupancy were related to breeding season length as well as the overall population recovery. These models provide evidence that swan habitat use is directly related to the average annual temperature conditions at a given location. Although we do not see strong trending in the temperature data used in our analysis, the relationships we observed indicate that range expansion has likely occurred in Alaska as a result of climate warming. This conclusion is based on other climate research conducted over a longer time frame indicating that the timing of spring breakup has been advancing (Keyser et al. 2000), and a phase shift in the Pacific Decadal Oscillation has resulted in generally warmer conditions in Alaska since 1976 (Hartmann & Wendler 2005). These changes could explain why more peripheral areas of the current trumpeter swan breeding range were not used until later years, after additional warming had occurred.

Support for a non-linear latitude effect is most likely explained by range edge effects. The population was concentrated at middle latitudes within our survey area during the first years of our survey. The negative quadratic effect of latitude could indicate that it took some time for colonization outside of these areas to occur, even though habitat was available. Waterfowl in general, and swans in particular, are highly philopatric (Anderson et al. 1992), so substantial time may be required for swans to colonize newly available habitat, which would lower the average occupancy probability. Juveniles are the most likely dispersers (Lindberg et al. 1998) and usually take 4-5 years to reach breeding age, and this lag could have slowed colonization even after habitats became suitable. Additionally, northern areas might not be as likely to be occupied because of variable weather conditions and the time needed for colonization to occur, and southern areas might also have suboptimal weather conditions during brood rearing (high rainfall levels and cool spring temperatures) making these areas less desirable for breeding. We suspect that weather patterns at range edges (e.g. cool and rainy in coastal areas) could counteract the positive effects of the number of days above freezing and help explain why occupancy is higher in the interior relative to the Arctic and coastal areas. Mountainous terrain also limits the number of suitable breeding wetlands in many areas of southern Alaska, explaining the lower occupancy rates in those areas.

The lack of precision in our estimates suggests uncertainty in the direction and magnitude of the effects. However, this is not surprising considering the partial correlation between latitude and days above freezing. We initially considered the days above freezing variable without latitude in the models, but this simpler model was not supported by the model selection results. In fact, models containing both latitude and days above freezing lowered the DIC by > 20 units over models containing either of these variables individually, indicating that both were important predictors. The correlation between these two variables made it difficult to precisely estimate the effects of each individual covariate, but together they provided a better representation of the data than either covariate alone. Inclusion of only latitude or season length may lead to erroneous or misleading conclusions because both appear to be important predictors. Regardless, together they indicate a strong positive effect of season length on occupancy probability. The support for a random year effect indicates that other sources of annual variation (e.g. precipitation and water levels) may also influence habitat occupancy, and the inclusion of this term accounted for these unmeasured variables.

Nearly 40 years of survey data, combined with the fact that 75-80% of the trumpeter swans in North America breed in Alaska, provided us with a unique opportunity to examine the effects of breeding season length over most of the breeding range of this species. For trumpeter swans, it appears that any negative effects of wetland shrinkage have so far been outweighed by the positive effects of increased available habitat due to extended ice-free periods. It is possible that some wetlands are no longer suitable due to shrinkage, while some deeper types, such as ponds and lakes, may be becoming shallower and therefore more suitable for breeding. If this is the case, we would expect range expansion to continue if climate warming continues. Schmidt et al. (2009b) found some evidence that cygnet production at lower latitudes is slowing, and this may indicate that breeding habitats are becoming saturated or that climate effects such as wetland shrinkage are affecting habitat suitability. For this population in particular, availability of sufficient wintering habitat may ultimately limit population growth (Mitchell 1994).

Based on our findings for trumpeter swans, we would expect continued climate warming to affect populations of boreal forest-nesting species in the future. Habitat requirements for swans are very similar to those of many boreal forest-nesting species, suggesting that the effects of climate warming on wetland-breeding birds could be widespread. Avian distributions are shifting northwards (La Sorte & Thompson 2007) and abundances of some waterfowl species in the Prairie Pothole Region are positively related to precipitation (Forcey et al. 2007), implying that pond loss would have negative population effects (Sorenson et al. 1998). Some species might be expected to experience changes in breeding phenology if seasons lengthen (Winkler et al. 2002, Møller et al. 2010), although these changes may not always have positive impacts. An advancement of spring conditions could cause a mismatch between breeding and food source availability (Both & Visser 2001, Inouye et al. 2000) or substantial changes in habitat (Ward et al. 2005). Either of these could negatively affect populations through reductions in breeding habitat and reproductive success. Although the future impacts of climate change on individual avian species are unknown, long-term monitoring programs like the trumpeter swan survey provide the best opportunity to identify population-level effects and should be encouraged. Because climate warming is expected to be more pronounced in northern regions, concentration of additional monitoring efforts in boreal forest habitats may help to identify species at risk of decline.