Changes in waterfowl migration phenologies in central North America: Implications for future waterfowl conservation.

Globally, migration phenologies of numerous avian species have shifted over the past half-century. Despite North American waterfowl being well researched, published data on shifts in waterfowl migration phenologies remain scarce. Understanding shifts in waterfowl migration phenologies along with potential drivers is critical for guiding future conservation efforts. Therefore, we utilized historical (1955-2008) nonbreeding waterfowl survey data collected at 21 National Wildlife Refuges in the mid- to lower portion of the Central Flyway to summarize changes in spring and autumn migration phenology. We examined changes in the timing of peak abundance from survey data at monthly intervals for each refuge and species (or species group; n = 22) by year and site-specific temperature for spring (Jan-Mar) and autumn (Oct-Dec) migration periods. For spring (n = 187) and autumn (n = 194) data sets, 13% and 9% exhibited statistically significant changes in the timing of peak migration across years, respectively, while the corresponding numbers for increasing temperatures were 4% and 9%. During spring migration, ≥80% of significant changes in the timing of spring peak indicated advancements, while 67% of significant changes in autumn peak timing indicated delays both across years and with increasing temperatures. Four refuges showed a consistent pattern across species of advancing spring migration peaks over time. Advancements in spring peak across years became proportionally less common among species with increasing latitude, while delays in autumn peak with increasing temperature became proportionally more common. Our study represents the first comprehensive summary of changes in spring and autumn migration phenology for Central Flyway waterfowl and demonstrates significant phenological changes during the latter part of the twentieth century.

Introduction

There has been a considerable amount of literature published in recent decades documenting changes in phenological events of plants and animals in response to climate change [1-4]. Much of this attention has been focused on shifts in avian migration phenologies [5-9], with most of the studies examining changes in the migration phenology of spring migrants [10, 11]. Based on numerous studies, it is clear that climate change is influencing migration phenologies of bird species across the globe. Moreover, it appears that the migration timing of many bird species is more affected by climate change during spring migration than autumn migration [7]. For spring migration, the dominating pattern has been that migration dates have advanced over time as global temperatures have increased, while the patterns for autumn migration have been more variable, though delays appear more common than advancements [5-9].

Each year, considerable resources are invested by federal, state, and provincial agencies as well as private organizations in the conservation and management of waterfowl in North America. Yet, to date, comparatively little data have been published on phenological shifts in waterfowl migration in general [12-14], and North American waterfowl in particular [15-19]. In fact, we found only 3 North American studies that presented data for more than 3 species of waterfowl [15, 18, 19] and just 2 studies that offered any data at all on autumn migration [15, 19]. What published data exist indicate that spring migration generally occurred earlier over time, while for autumn migration, delays in migration timing were more common than advancements. Moreover, model projections of reduced snow cover and rising air temperatures under a warming climate are expected to lead to northward shifts in wintering ranges and delayed autumn migration for several dabbling duck species in the eastern third of North America [20]. Further, the rate of change in migratory passage by year or temperature may vary with latitude [8, 21].

The current lack of published data on basic changes in migration phenology for North American waterfowl is problematic considering that successful waterfowl management and conservation planning is dependent on accurate estimates of the spatio-temporal distribution of waterfowl during all stages of their life-cycles [22, 23]. Without such knowledge, it becomes nearly impossible to ensure that specific demands for resources (e.g., energy, necessary nutrients, roosting habitat) by migrating waterfowl can be met at appropriate points in time and space. This is of particular importance as mismatches between the need for resources and availability of those resources can have detrimental effects on survival and reproduction [22, 24–27]. Moreover, detailed knowledge of changes in migration phenology is necessary for guiding waterfowl harvest policies, particularly setting timing for seasons such that they correspond with occurrence of waterfowl. Therefore, it is important to understand how migration phenologies have changed in the past and how they may change in the future.

Waterfowl management in North America is administered within 4 major flyways: the Atlantic, Mississippi, Central, and Pacific flyways. Our study focused on the mid- to lower portion of the Central Flyway and encompassed Texas, Oklahoma, Kansas, Nebraska, and the eastern portion of New Mexico (east of the continental divide). Within this region, the U.S. Fish and Wildlife Service manages numerous conservation units (hereafter, refuges) many of which were established or are managed primarily for waterfowl conservation [28, 29]. As global warming continues, information on waterfowl migration phenologies within this flyway and how they may relate to changes in the climate will be critical for guiding future waterfowl harvest policies, management, and conservation strategies throughout central North America. As a significant number of the refuges within this region either currently monitor their waterfowl numbers or have done so in the past, they provide an opportunity to examine long-term changes in waterfowl migration phenology in this flyway. In this study, we utilized historical (1955–2008) nonbreeding waterfowl survey data from the Central Flyway to produce the first comprehensive summary of phenological changes in spring and autumn waterfowl migration for this region. Changes in spring and autumn migration phenology were evaluated with respect to time and site-specific temperatures.

Methods

Data collection and processing

We gathered all available nonbreeding waterfowl survey data from 21 refuges providing at least 5 years of data within the mid- to lower portion of the Central Flyway (see [30] for detailed descriptions of data acquisition and entry). Survey data consisted of aerial and ground counts conducted during 1955–2008; data spanned from 7 to 54 years (mean ± SD: 26.6 ± 12.7) with the number of years with data ranging from 7 to 47 (mean ± SD: 24.0 ± 9.4) among sites. Only sites at mid- to lower latitudes were included because refuges at higher latitudes harbor many breeding waterfowl that makes it difficult to identify migrating waterfowl. Because few refuges conducted regular surveys during September or April, we restricted data to October–December for the autumn-winter (hereafter, autumn) migration period and January–March for the winter-spring (hereafter, spring) migration period to obtain comparable data for as many refuges as possible, while still including most of the autumn and spring migrations for the majority of waterfowl species of the Central Flyway [31]. Further details on the waterfowl survey data can be found in [30, 32].

We included all waterfowl species for which data were available in our analyses. However, due to the relatively recent split of the Canada goose complex into 2 different species [33], very few surveys in our study distinguished between Canada Geese (Branta canadensis) and Cackling Geese (B. hutchinsii). We therefore combined these 2 species into 1 group (Canada geese) for all analyses. Similarly, only 1 refuge differentiated between Snow Geese (Anser caerulescens) and Ross’s Geese (A. rossii); therefore, we analyzed them separately for this refuge, but treated them as a single group (light geese) for all other refuges. We also combined Greater and Lesser Scaup (Aythya marila and A. affinis; hereafter, scaup) and Common and Barrow’s Goldeneyes (Bucephala clangula and B. islandica; hereafter, goldeneyes) because they are closely related and difficult to visually separate in the field. All group sums (Canada geese, light geese, scaup, and goldeneyes) were calculated on raw count data for each individual survey before any other calculations were performed. We propagated missing values through summations so that any sum based upon a missing count value resulted in a missing value.

All individual surveys that were found to be incomplete (i.e., when only a portion of the usual survey area was surveyed) were excluded from all analyses. Following [32], we considered counts for which the proportion of unidentified birds was indeterminable or exceeded 10% as too unreliable and therefore excluded those data from all analyses. A few sites experienced minor increases in survey area over time (see S1–S4 Tables), and if not correctable by subtracting the counts for the added area, we used uncorrected counts. However, in no case was there any discernible effect on the relative number of counted individuals of any species among months as a result of increased survey area (Kent Andersson, Oklahoma State University, unpublished data). The frequency of waterfowl surveys varied within and among sites, from weekly to monthly surveys. Eleven sites offered only monthly surveys and only 3 sites consistently offered weekly surveys for greater than 50% of their time series. The timing of surveys within months often varied within sites as well. Therefore, we calculated monthly averages to make data consistent among and within sites.

To assess changes in migration phenology over time, we used changes in the timing of peak abundance, as indexed by the monthly count averages, as an indicator of population-level changes in migration timing. We defined peak abundance month as the month with the greatest count average within each refuge, species, and year for spring and autumn migration periods separately (hereafter, spring peak and autumn peak, respectively). Thus, for each specific refuge, species, year, and migration season, the month with the greatest count average within the 3-month period was considered the peak abundance month. Any individual spring or autumn migration period for any refuge, species, and year that did not provide count averages for all 3 months or for which a single peak month could not be identified (i.e., the greatest count average was represented in more than 1 month) was excluded from all analyses.

Peak abundance month exhibits several traits that makes it suitable for analyzing changes in migration phenology at the population level: 1) it represents when a large fraction of the population that uses a particular site can be found there and therefore, likely reflects high resource demands at that site; 2) it is robust to changes in survey effort and population size, which date of first arrival or last departure are not [5, 34–36]; and 3) unlike the mean or median passage date, peak abundance month is also robust against truncated distributions, such as the current survey data, as long as the distribution of monthly averages is unimodal or the highest peak occurs within the defined spring or autumn migration periods. In our case, this would likely be true for most refuges and species. Peaks occurring outside the defined migration periods could potentially present a problem. However, because we treated peak month as an ordered categorical variable, as long as the highest peak within the defined migration period occurs in the month closest to the true peak outside the defined period, the validity and interpretation of the results remain unchanged. Cases like these merely change the interpretation of the first or last months (i.e., categories) to mean that the peak occurred in or before the first month or in or after the last month, respectively, rather than just within the month itself. Further, an evaluation of all the refuges in our study offering any survey data outside October–March, revealed bimodality with the highest peak outside the defined periods were rare, and cases where the internal peak did not fall at the relevant boundary occurred in only 3% of cases (see S1 Text for further details). For specific species where particular refuges serve as wintering sites, spring and autumn peaks naturally represent peaks in the number of wintering birds rather than actively migrating birds. However, as the implications are the same from the perspective of migration phenology (i.e., an earlier peak indicates advanced arrival of migrating birds and a later peak indicates a delayed arrival), we did not distinguish between them here.

The more sporadically a specific waterfowl species occurs at a given location and the lower the average peak abundance, the more likely it is that random events will dominate any observed patterns in peak abundance, obscuring actual changes in migration phenology. Therefore, we limited phenology analyses to species-refuge combinations with a total of at least 5 non-zero data points, less than 25% zeros, and an average peak count of 200 birds or more for the migration season considered. Of the 21 sites, 2 were located in Nebraska (North Platte and Crescent Lake), 3 in Kansas (Kirwin, Flint Hills, and Quivira), 4 in Oklahoma (Salt Plains, Washita, Deep Fork, and Tishomingo), 2 in New Mexico (Bosque del Apache and Bitter Lake), and 10 in Texas (Texas Point, Attwater Prairie Chicken, McFaddin, Anahuac, Brazoria, San Bernard, Big Boggy, Matagorda Island, Aransas, and Laguna Atascosa; Fig 1). Although Matagorda Island is technically a unit of Aransas, it was surveyed independently and therefore, we considered it separately. All 21 sites were designated as National Wildlife Refuges by the U.S. Fish and Wildlife Service.

Fig 1

Locations of the National Wildlife Refuges included in this study.

Abbreviation represents Texas Point (TX Point).

By focusing on peak abundance month (as determined by the greatest monthly count average), the only assumption needed for our analyses was that surveys were methodologically and spatially consistent within each distinct 3-month period. Most of the included data sets show a relatively high degree of methodological and spatial consistency over their entire data series [32]. This consistency becomes even greater when differences between years and migration seasons are not considered (Kent Andersson, Oklahoma State University, unpublished data). Moreover, even the absence of a written survey protocol (a fairly common problem among these surveys; [32]) is unlikely to lead to such vast differences in methodology over as short a time interval as 3 consecutive months as to alter the rank order of the months. Hence, we were confident that the identified peak month corresponded to the month with the greatest average number of waterfowl for each species and location.

To assess the influence of site-specific temperature on waterfowl migration phenology, we downloaded daily maximum temperatures for each refuge centroid from the PRISM (Parameter-elevation Regressions on Independent Slopes Model) [37] data set for the spring (January–March) and autumn (October–December) migration periods for all relevant years (PRISM Climate Group, Oregon State University, http://prism.oregonstate.edu, created 19 April 2021). Because data prior to 1981 were not available from the PRISM data set, we downloaded daily maximum temperatures for the nearest weather station with available data for years prior to 1981 from the National Oceanic and Atmospheric Administration’s Global Historical Climatology Network—Daily Database [38, 39] (https://www.ncdc.noaa.gov/cdo-web/, accessed 19 April 2021). To ensure these temperature data were comparable to PRISM data, we also downloaded overlapping temperature data for 1981–2008 and performed linear regressions on monthly averages for the overlapping data; and because all slopes were near 1.0 (range: 0.98–1.01), we considered these data sources comparable. We then calculated the average daily maximum temperature for each migration season, year, and refuge. Thus, each resulting temperature average corresponded to an entire spring or autumn migration season for a specific year and refuge.

Data analyses

We analyzed each refuge and species for changes in spring (n = 187) and autumn (n = 194) peak by year and average daily maximum temperature separately using PROC LOGISTIC in SAS 9.4 (SAS Institute, Inc., Cary, NC). Data were either binomial (i.e., when all peaks fell within 2 months; spring: n = 67, autumn: n = 69) or trinomial (i.e., when all 3 months were represented; spring: n = 100, autumn: n = 114). For the purpose of statistical analysis, the cases where all peaks fell within the same month were trivial, as no temporal change in peak could be detected in those cases (spring: n = 20, autumn: n = 11).

Binomial data sets were analyzed using logistic regression with penalized maximum likelihood estimation. The main reason for using penalized maximum likelihood was to obtain finite maximum likelihood estimates for data sets with complete separation [40] (n = 4 and n = 2 with year and average maximum daily temperature as explanatory variables, respectively, for both spring and autumn), but it has the added benefit of reducing bias resulting from small sample size [41]. Trinomial data sets were analyzed using ordinal logistic regression with a proportional odds version of the cumulative logit model [42]. This was based on the assumption of an underlying latent continuous response variable (i.e., timing of abundance peak) that satisfied an ordinary regression model with a similar dispersion across the continuum of the predictor variables (i.e., year or temperature) [43-45]. Such cumulative logit models are sensitive to location effects rather than changes in variability with changing values of the predictor variable [42]. Thus, this assumption seemed reasonable for our data where the underlying latent variable, timing of abundance peak, undoubtedly was continuous and unlikely to differ dramatically in dispersion over relevant years or temperatures. Moreover, a score test of the proportional odds assumption (available in PROC LOGISTIC [46]) performs poorly when data are sparse, such as when comparatively few observations fall in 1 of the response categories or in the case of continuous predictors [45]. As our predictor variables were continuous and data were often sparse in 1 response category, the validity of any score test for the assumption of proportional odds would be highly questionable. Therefore, we assumed the proportional odds assumption to be valid for all trinomial data sets a priori.

Penalized maximum likelihood estimation for ordinal logistic regression was not supported in SAS 9.4 [46]. We therefore transformed the single trinomial data set that exhibited complete separation (i.e., spring peaks by year for Ross’s Goose at Bitter Lake) to a binary data set by combining the single March peak with the 2 February peaks into 1 category before analysis. Because the underlying data exhibited a monotonic increase in peak month over time, this transformation was conservative. For ordinal categorical data, regression parameters often have asymmetric distributions when most observations fall in the highest or lowest category, when sample size is small, or when complete separation exists [40, 42]. Because many of our data sets exhibited 1 or more of these characteristics, we used likelihood ratio tests to identify statistically significant trends and based all confidence intervals on profile likelihood estimates [40, 42]. All modeled probabilities were cumulated over higher values of peak month; thus, a positive slope indicates that abundance peaks later in the season over time or with increasing temperature and a negative slope that it peaks earlier.

The power of a statistical test is dependent on the sample size with a smaller sample size resulting in lower power. Further, in ordinal logistic regression, fewer possible outcome categories also translate into lower power [42]. As our sample sizes were generally small and the number of outcome categories for the ordinal logistic regressions merely 3, it was likely that smaller shifts in migration phenology would not manifest as statistically significant. Hence, to identify general patterns for either year or temperature that may otherwise go undetected, we followed [1] and performed exact binomial tests on the distribution of all positive and negative slope values derived from the logistic regressions, regardless of whether statistical significance was obtained during the regression. We did this for regressions with year and temperature separately. Our null hypothesis was equal probability (i.e., P = 0.5) of a non-zero slope being positive or negative (only non-zero slopes were included in these analyses as slope values of zero were neither positive nor negative [1, 47]). For each explanatory variable (year and temperature), we performed this analysis on the entire data set, as well as on individual refuges across all species and on individual species across all refuges for spring and autumn migration periods, separately.

To investigate if the prevalence of advancements or delays in peak month by year or temperature was related to latitude, we used linear regression of the proportion of non-zero slopes (obtained from the logistic regressions) that indicated advancement (i.e., negative slopes) or delays (i.e., positive slopes) for each refuge and predictor variable as a function of refuge centroid latitude. In all cases, scatterplots of standardized predicted values vs. standardized residuals for the linear regressions indicated that these data met the assumptions of homogeneity of variance and linearity and the residuals were approximately normally distributed. The significance level was set to α = 0.05. All tests were two-tailed.

Results

Shifts over time

Of 187 spring data sets, 13% exhibited statistically significant changes in the timing of peak migration across years. The dominant pattern of change during spring was towards earlier migration peaks over time, both across all non-zero spring slopes (61% negative slopes, n = 167, exact binomial test: P = 0.005; S1 Table) and among statistically significant trends (80% negative slopes, n = 25). When examining non-zero slopes for individual species separately across all refuges, only Green-winged Teal (Anas crecca) displayed a dominance of slopes that was statistically significant with 78% of slopes indicating earlier spring peaks over time (Table 1). For individual refuges separately across all species, Aransas, Attwater Prairie Chicken, Matagorda Island, and Quivira displayed significantly skewed distributions of positive and negative trends, all of which indicated a dominance of shifts towards earlier spring peaks over time (83%, 100%, 100%, and 92% negative trends, respectively; Table 1). Species trends showing shifts towards earlier spring peaks over time became proportionally less common with increasing refuge latitude (least squares linear regression: b = -0.027, r = 0.20, F = 4.58, P = 0.046; Fig 2).

Table 1

Skewness in the distribution of slopes indicating advancements or delays in spring peak across years by waterfowl species and refuge.

Species	n a	% neg b	P c	Refuge	n a	% neg b	P c
American Wigeon	14	43	0.79	Anahuac	11	55	1.00
Blue-winged Teal	11	82	0.07	Aransas	12	83	0.04
Bufflehead	1	100	1.00	Attwater Prairie Chicken	7	100	0.02
Canada geese	12	42	0.77	Big Boggy	8	38	0.73
Canvasback	4	75	0.63	Bitter Lake	6	17	0.22
Cinnamon Teal	1	0	1.00	Bosque del Apache	9	33	0.51
Common Merganser	5	40	1.00	Brazoria	9	67	0.51
Gadwall	13	77	0.09	Crescent Lake	1	0	1.00
Goldeneyes	1	0	1.00	Flint Hills	8	75	0.29
Greater White-fronted Goose	9	67	0.51	Kirwin	10	30	0.34
Green-winged Teal	18	78	0.03	Laguna Atascosa	10	70	0.34
Light geese	14	71	0.18	Matagorda Island	8	100	0.008
Mallard	10	40	0.75	McFaddin	11	73	0.23
Northern Pintail	16	50	1.00	North Platte	6	33	0.69
Northern Shoveler	14	71	0.18	Quivira	13	92	0.003
Redhead	8	88	0.07	Salt Plains	14	64	0.42
Ring-necked Duck	1	0	1.00	San Bernard	7	57	1.00
Ross’s Goose	1	0	1.00	Texas Point	5	60	1.00
Ruddy Duck	3	67	1.00	Tishomingo	5	0	0.06
Scaup	10	50	1.00	Washita	7	57	1.00
Snow Goose	1	0	1.00

a Only non-zero slopes were included as slope values of zero were neither positive nor negative [1, 47].

b Percent non-zero logistic regression slopes that indicated advancement of spring abundance peak across years (i.e., negative slopes).

c P-value from exact binomial test under the null hypothesis of equal probability (i.e., P = 0.5) of a non-zero slope being positive or negative. Tests were performed for each waterfowl species across all refuges and for each refuge across all species.

Fig 2

Proportion of species showing advancing spring peak across years by refuge latitude.

Only non-zero slopes were included. Trend line obtained by least squares linear regression.

For autumn, 9% of data sets (n = 194) exhibited statistically significant changes in the timing of peak migration across years, with 67% of significant trends (n = 18) indicating peaks occurring later over time. There was no obvious pattern in peak shifts when all non-zero slopes were considered (52% positive trends, n = 183, exact binomial test: P = 0.55; S2 Table). No individual species exhibited a statistically significant skew in its autumn distributions of positive and negative trends across all refuges (Table 2). When considering individual refuges separately across all species, only Matagorda Island showed a significantly skewed distribution of positive and negative trends with 100% of slopes indicating autumn peaks occurring later over time (Table 2). Species trends showing advancing autumn peaks over time showed a tendency to become proportionally less common with increasing refuge latitude, but the relationship was not statistically significant (least squares linear regression: b = 0.015, r = 0.10, F = 2.04, P = 0.17; Fig 3).

Table 2

Skewness in the distribution of slopes indicating advancements or delays in autumn peak across years by waterfowl species and refuge.

Species	n a	% pos b	P c	Refuge	n a	% pos b	P c
American Wigeon	18	50	1.00	Anahuac	10	80	0.11
Blue-winged Teal	12	50	1.00	Aransas	12	58	0.77
Bufflehead	1	0	1.00	Attwater Prairie Chicken	6	50	1.00
Canada geese	10	40	0.75	Big Boggy	7	14	0.13
Canvasback	4	25	0.63	Bitter Lake	8	38	0.73
Cinnamon Teal	1	0	1.00	Bosque del Apache	9	44	1.00
Common Merganser	3	33	1.00	Brazoria	11	45	1.00
Gadwall	16	44	0.80	Crescent Lake	3	67	1.00
Goldeneyes	3	33	1.00	Deep Fork	2	100	0.50
Greater White-fronted Goose	12	75	0.15	Flint Hills	9	56	1.00
Green-winged Teal	17	35	0.33	Kirwin	11	64	0.55
Light geese	13	69	0.27	Laguna Atascosa	11	27	0.23
Mallard	12	58	0.77	Matagorda Island	8	100	0.008
Northern Pintail	16	56	0.80	McFaddin	11	36	0.55
Northern Shoveler	17	47	1.00	North Platte	7	71	0.45
Redhead	10	80	0.11	Quivira	13	62	0.58
Ring-necked Duck	3	67	1.00	Salt Plains	14	64	0.42
Ross’s Goose	1	0	1.00	San Bernard	9	33	0.51
Ruddy Duck	3	67	1.00	Texas Point	6	33	0.69
Scaup	9	56	1.00	Tishomingo	10	60	0.75
Snow Goose	1	100	1.00	Washita	6	17	0.22
Wood Duck	1	100	1.00

a Only non-zero slopes were included as slope values of zero were neither positive nor negative [1, 47].

b Percent non-zero logistic regression slopes that indicated delayed autumn abundance peak across years (i.e., positive slopes).

c P-value from exact binomial test under the null hypothesis of equal probability (i.e., P = 0.5) of a non-zero slope being positive or negative. Tests were performed for each waterfowl species across all refuges and for each refuge across all species.

Fig 3

Proportion of species showing delayed autumn peak across years by refuge latitude.

Only non-zero slopes were included. Trend line obtained by least squares linear regression.

Shifts in relation to temperature

The dominant pattern during spring was for higher temperatures being associated with earlier migration peaks, both across all spring data sets (63% negative slopes, n = 167, exact binomial test: P = 0.001; S3 Table) and among statistically significant trends (88% negative slopes, n = 8). However, only 4% of spring data sets (n = 187) exhibited statistically significant changes in the timing of peak migration with increasing temperatures. When examining positive and negative trends for individual species separately across all refuges, only Northern Pintail (Anas acuta) displayed a dominance of slopes that was statistically significant with 81% of slopes indicating earlier spring peaks with increasing temperatures (Table 3). For individual refuges across all species, Attwater Prairie Chicken and Kirwin displayed significantly skewed trend distributions of positive and negative slopes, with both showing earlier spring peaks with increasing temperatures as the dominant pattern (100% and 90% negative slopes, respectively; Table 3). Refuge latitude showed no discernable relationship with the proportion of species trends showing advancement in spring peak with increasing temperature (least squares linear regression: b = -0.007, r = 0.03, F = 0.46, P = 0.5; Fig 4).

Table 3

Skewness in the distribution of slopes indicating advancements or delays in spring peak with increasing temperature by waterfowl species and refuge.

Species	n a	% neg b	P c	Refuge	n a	% neg b	P c
American Wigeon	14	57	0.79	Anahuac	11	45	1.00
Blue-winged Teal	11	45	1.00	Aransas	12	75	0.15
Bufflehead	1	100	1.00	Attwater Prairie Chicken	7	100	0.02
Canada geese	12	58	0.77	Big Boggy	8	50	1.00
Canvasback	4	50	1.00	Bitter Lake	6	67	0.69
Cinnamon Teal	1	0	1.00	Bosque del Apache	9	33	0.51
Common Merganser	5	100	0.06	Brazoria	9	67	0.51
Gadwall	13	62	0.58	Crescent Lake	1	0	1.00
Goldeneyes	1	100	1.00	Flint Hills	8	50	1.00
Greater White-fronted Goose	9	56	1.00	Kirwin	10	90	0.02
Green-winged Teal	18	72	0.10	Laguna Atascosa	10	70	0.34
Light geese	14	71	0.18	Matagorda Island	8	38	0.73
Mallard	10	40	0.75	McFaddin	11	64	0.55
Northern Pintail	16	81	0.02	North Platte	6	67	0.69
Northern Shoveler	14	64	0.42	Quivira	13	77	0.09
Redhead	8	63	0.73	Salt Plains	14	79	0.06
Ring-necked Duck	1	0	1.00	San Bernard	7	71	0.45
Ross’s Goose	1	100	1.00	Texas Point	5	40	1.00
Ruddy Duck	3	100	0.25	Tishomingo	5	60	1.00
Scaup	10	50	1.00	Washita	7	29	0.45
Snow Goose	1	0	1.00

a Only non-zero slopes were included as slope values of zero were neither positive nor negative [1, 47].

b Percent non-zero logistic regression slopes that indicated advancement of spring abundance peak with increasing site-specific average Jan–Mar daily maximum temperature (i.e., negative slopes).

c P-value from exact binomial test under the null hypothesis of equal probability (i.e., P = 0.5) of a non-zero slope being positive or negative. Tests were performed for each waterfowl species across all refuges and for each refuge across all species.

Fig 4

Proportion of species showing advancing spring peak with increasing temperature by refuge latitude.

Only non-zero slopes were included. Trend line obtained by least squares linear regression.

For autumn, 9% of data sets (n = 194), exhibited statistically significant changes in the timing of peak migration with increasing temperature, with 67% of significant trends (n = 18) indicating peaks occurring later with increasing temperature. There was no obvious pattern in peak shifts when all non-zero slopes were considered (54% positive trends, n = 183, exact binomial test: P = 0.30; S4 Table). No species exhibited a statistically significant skew in its autumn distributions of positive and negative trends across all refuges (Table 4). When considering individual refuges separately across all species, only Kirwin showed a significantly skewed distribution of positive and negative trends with 91% of slopes indicating autumn peaks occurring later with increasing temperature (Table 4). Species trends indicating delayed autumn peaks with increasing temperature became proportionally more common with increasing refuge latitude (least squares linear regression: b = 0.023, r = 0.30, F = 8.20, P = 0.010; Fig 5).

Table 4

Skewness in the distribution of slopes indicating advancements or delays in autumn peak with increasing temperature by waterfowl species and refuge.

Species	n a	% pos b	P c	Refuge	n a	% pos b	P c
American Wigeon	18	39	0.48	Anahuac	10	50	1.00
Blue-winged Teal	12	67	0.39	Aransas	12	58	0.77
Bufflehead	1	100	1.00	Attwater Prairie Chicken	6	33	0.69
Canada geese	10	60	0.75	Big Boggy	7	57	1.00
Canvasback	4	25	0.63	Bitter Lake	8	50	1.00
Cinnamon Teal	1	0	1.00	Bosque del Apache	9	44	1.00
Common Merganser	3	100	0.25	Brazoria	11	45	1.00
Gadwall	16	44	0.80	Crescent Lake	3	100	0.25
Goldeneyes	3	67	1.00	Deep Fork	2	50	1.00
Greater White-fronted Goose	12	58	0.77	Flint Hills	9	44	1.00
Green-winged Teal	17	65	0.33	Kirwin	11	91	0.01
Light geese	13	62	0.58	Laguna Atascosa	11	55	1.00
Mallard	12	75	0.15	Matagorda Island	8	75	0.29
Northern Pintail	16	50	1.00	McFaddin	11	45	1.00
Northern Shoveler	17	29	0.14	North Platte	7	86	0.13
Redhead	10	60	0.75	Quivira	13	46	1.00
Ring-necked Duck	3	33	1.00	Salt Plains	14	57	0.79
Ross’s Goose	1	0	1.00	San Bernard	9	22	0.18
Ruddy Duck	3	67	1.00	Texas Point	6	17	0.22
Scaup	9	67	0.51	Tishomingo	10	70	0.34
Snow Goose	1	100	1.00	Washita	6	50	1.00
Wood Duck	1	0	1.00

a Only non-zero slopes were included as slope values of zero were neither positive nor negative [1, 47].

b Percent non-zero logistic regression slopes that indicated delayed autumn abundance peak with increasing site-specific average Oct–Dec daily maximum temperature (i.e., positive slopes).

c P-value from exact binomial test under the null hypothesis of equal probability (i.e., P = 0.5) of a non-zero slope being positive or negative. Tests were performed for each waterfowl species across all refuges and for each refuge across all species.

Fig 5

Proportion of species showing delayed autumn peak with increasing temperature by refuge latitude.

Only non-zero slopes were included. Trend line obtained by least squares linear regression.

Discussion

Given that Earth’s climate is projected to continue changing as a result of ongoing global warming [48], it is critical that we identify and understand changes in migration phenologies of migratory birds. For migratory waterfowl, such information is particularly important because stopover habitats used during migration are critical for acquisition of nutrients necessary for continuing migration, successful reproduction, and survival [22, 24–27]. Changing migration phenologies may generate mismatches in their ability to rely on such habitats. Our results indicated that the timing of waterfowl migration in the mid- to lower portion of the Central Flyway has undergone significant changes since the 1950s. Statistically significant trends for changes in peak abundance month were relatively rare during both autumn and spring periods (4–13% of all data sets). However, the dominant pattern for changes in spring migration was that of peak abundance month occurring earlier over time and with increasing temperatures, both for statistically significant trends only and across all non-zero slopes. During autumn migration, only statistically significant trends showed any discernable pattern, with delays in peak abundance month being more common than advancements over time and with increasing temperatures.

Given the coarseness of the temporal resolution of the dependent variable (i.e., peak abundance month) and the necessity to treat it as an ordinal categorical variable with only 3 possible outcomes for each migration period, it was not unexpected that relatively few species-refuge trends manifested as statistically significant. Even with our generally relatively long survey series (mean ± SD: 26.6 ± 12.7 years), only very strong relationships would likely result in statistically significant trends. Given this, the number of detected phenological changes was surprisingly high and suggests that rather dramatic phenological shifts in waterfowl migration have occurred in North America during the latter part of the twentieth century.

With global temperatures expected to continue increasing [48], we expect migration timing to continue changing for many waterfowl species, leading to gradual changes in the timing of their occurrence at specific refuges. Thus, for individual refuges, resource management may need to be adjusted such that availability of resources match waterfowl abundances and requirements. Such optimization of refuge management would maximize conservation as well as fiscal outcomes. Further, the timing of local harvest seasons must coincide with presence of waterfowl in great enough numbers to allow rewarding hunting opportunities or it may negatively affect hunter participation. This, in turn, may directly lead to reduced funding for waterfowl conservation through loss of revenue from sales of Federal Migratory Bird Hunting Conservation Stamps (i.e., duck stamps; [49]).

Numerous studies have documented advancement in spring migration in birds including some waterfowl [5-9]. Our results, indicating a widespread advancement in the timing of spring migration for waterfowl in the Central Flyway that appears connected to increasing temperatures, agree well with the general consensus of these previous studies. The apparent northward shifts in wintering distribution in several waterfowl species over the last century [31, 50] could explain at least some of this pattern as it may greatly accelerate the apparent migratory passage [8]. We also found that the proportion of trends indicating advancements versus delays in spring peak decreased with increasing latitude. This appears congruent with the findings of [8].

Mismatches between the need for and availability of resources required for successful breeding can have detrimental effects on a species reproductive success [51-53]. Hence, advancing spring waterfowl migration may be explained as an adaptive response if the conditions required for successful breeding, such as the availability of food and open water, occurred earlier over time. Indeed, the last century has seen increasing temperatures [48] concomitant with a gradual advancement of spring [54-56] and corresponding advancements of phenological events in many plants and animals [2, 3, 57–59] at a range of latitudes. However, as migrating waterfowl cannot have direct knowledge of the conditions at breeding areas, this also requires corresponding advancements in migration cues or needed resources along the migration route, or detrimental mismatches may still result [51].

The costs and benefits associated with and ecological and environmental constraints on autumn migration timing are poorly understood and likely highly variable depending on the species’ life history strategy [7, 10]. Generally, waterfowl migrate from breeding grounds to avoid deteriorating conditions linked to declining temperatures as winter approaches [60]. Indeed, temperature and snow cover as well as associated large-scale weather patterns have been identified as important factors affecting the timing of autumn migration for several waterfowl species [19, 61–64]. While temperature and related weather phenomena may be the ultimate drivers of autumn migration, many other factors, such as wind speed and direction [63, 65, 66], precipitation [66], human disturbance [67, 68], competition [69, 70], and predation pressure [71] could influence a specific species’ autumn migration phenology, leading to complex patterns [72].

Unfortunately, published data on changes in autumn migration phenology for waterfowl are scarce. What data exist suggest that there is much variation but delays appear more common than advancements [8, 13, 18, 19]. Our results confirmed this general pattern with significant delays being twice as common as advances, and with seasonal averages of local temperatures appearing to have some explanatory value. That delays in autumn peaks appeared more affected by temperature at higher latitudes support the theory that temperatures below freezing and snow cover are important drivers behind autumn migration for many waterfowl species in central North America. As the likelihood of lasting snow cover and extended periods of temperatures below freezing decreases with decreasing latitude [38, 39] (https://www.ncdc.noaa.gov/cdo-web/, accessed 16 July 2021), the effects of warming temperatures on the timing of autumn migration may be less pronounced at lower latitudes.

Given the variation in latitude and data series length and time span among refuges, it was not surprising that few species showed significantly skewed trend distributions across refuges. It was, however, noteworthy that 4 refuges exhibited significantly skewed trend distributions across species towards advancing spring peaks over time. This intimates that advancing spring migration may be widespread among species.

By necessity, studies of migration phenologies for larger birds like waterfowl generally rely on counts of actively migrating birds at stopover sites and therefore miss any birds that do not stop at the survey location, as well as individual birds that arrive and depart between consecutive surveys. This was the case for most waterfowl counts included in our analyses as well. But, it is generally assumed that counts of birds at stopover sites provide a good index of the actual number of birds that migrate past, although this remains an untested assumption [73].

Survey data from stopover or wintering sites may also involve different populations with different migration phenologies of the same species whose population sizes have changed asynchronously over time. This can cause apparent phenological shifts even if the respective migration phenologies remain unchanged over time. Unfortunately, reliable data for historical changes in different sub-populations of waterfowl species utilizing the Central Flyway appear rare to non-existent [31]. Given the large differences in biology of different waterfowl species, we judged this scenario to be an unlikely explanation of the generally consistent patterns across many species. Similarly, large-scale changes in relative dominance by groups divided by age, sex, or breeding investment (single versus multiple broods and successful versus unsuccessful nesting attempts) are unlikely to offer a cogent explanation of the observed patterns.

Our study represents the first large-scale analysis of changes in migration phenology for waterfowl in the Central Flyway. We show that waterfowl migrating within the Central Flyway have experienced significant changes in their spring and autumn migration phenology since the 1950s and that these changes appear linked to warming temperatures. As migration phenologies change, it is critical that management and conservation efforts are realigned to match the changing spatiotemporal distribution of waterfowl; otherwise, mismatches with critical food and habitat resources needed by migrating and wintering waterfowl may occur.

Phenological changes in spring abundance peak across years for waterfowl at Central Flyway National Wildlife Refuges.

(DOCX)

Click here for additional data file.

Phenological changes in autumn abundance peak across years for waterfowl at Central Flyway National Wildlife Refuges.

(DOCX)

Click here for additional data file.

Phenological changes in spring abundance peak with increasing temperature for waterfowl at Central Flyway National Wildlife Refuges.

(DOCX)

Click here for additional data file.

Phenological changes in autumn abundance peak with increasing temperature for waterfowl at Central Flyway National Wildlife Refuges.

(DOCX)

Click here for additional data file.

Further details on the choice of phenological parameter.

(DOCX)

Click here for additional data file.