The Northern Harrier (Circus hudsonius, hereafter harrier) is regarded as a dietary generalist (Smith et al. 2011). However, the breeding productivity of harriers in Wisconsin (Hamerstrom et al. 1985) and New Brunswick (Simmons et al. 1986) has been linked to the multiyear cycle in springtime local abundance of meadow voles (Microtus pennsylvanicus, hereafter vole), one of the harriers' primary food sources. The harrier-vole association is one in which harrier reproduction is stimulated to the point that its courtship behavior changes with high vole abundance in the spring (Hamerstrom 1986). The potential for cycles in harrier productivity that are synchronous across areas larger than the Wisconsin or New Brunswick study sites seems to have had little formal examination. Berger and Taylor (2003) reported that breeding harrier numbers in southern Manitoba tended to have local peaks about 8 yr apart, but synchrony across the province was weak. The North American Breeding Bird Survey (Sauer et al. 2013) has insufficient data on harrier numbers for the higher latitudes of the harrier breeding range. In autumn, however, harriers migrate almost completely out of the northernmost part of their summer range (Bildstein 2006, Smith et al. 2011), and the annual tallies of autumn migrants that have recently become available (Hawk Migration Association of North America 2018) may now be used to interpret population trends (Farmer et al. 2007). It is expected that each major migration watch site in the north may tally autumn migrants from a large summer region, and that a composite of tallies from longitudinally distant sites would be an indicator of harrier abundance across a wide area at higher latitudes. Although long-term cycles in autumn harrier counts have not been clearly documented, positive interannual correlations among counts from different watch sites have been found (Farmer et al. 2008). Synchrony of local changes in vole abundance across large distances in North America has not been established (Krebs 2013).

We studied this record of autumn counts of migrant harriers for possible evidence of cycles synchronous across broad areas by using autocorrelation analyses (Turchin 2003) to test the null hypothesis that the counts were not cyclic. Further, we tested for correlation between spring counts and autumn counts of migrant harriers in the same region to seek inferences about breeding productivity. We assessed the year-to-year differences in regional atmospheric conditions during migration as one possible influence on the numbers of harriers counted. Variation in the counting process from year to year could also influence the count pattern, so we compared counts for other diurnal raptor species to those for harriers at the same sites.

Methods

We analyzed 35 consecutive years of autumn migration tallies from four major observation sites in northeastern North America where harrier count totals were substantial and expected to be independent of each other within the same autumn: Hawk Ridge Bird Observatory (Hawk Ridge), Holiday Beach Migration Observatory (Holiday Beach), Hawk Mountain Sanctuary (Hawk Mountain), and Cape May Bird Observatory (Cape May; Fig. 1, Table 1). We also compared autumn counts with counts of north-bound harriers at Derby Hill Bird Observatory (Derby Hill; Fig. 1, Table 1) from both the preceding spring and the following spring. For the region containing these five sites and the land north of them, the North Atlantic Oscillation (NAO) is the predominant broad-scale atmospheric pattern (Hurrell and Deser 2009).

Figure 1.

Location of five raptor count sites for the Northern Harrier migration cycle study, 1983–2017. Sites denoted by a circle had counts made during autumn migration, and the site denoted by a triangle had counts made during spring migration.

Table 1.

Locations and observation intervals for counts of migrating raptors.

Total visual counts of harriers passing in flight (  Supplemental Materials (38_rapt-54-01-04_s03.pdf)) were from Hawk Migration Association of North America (2018) for 1983–2017 (autumn) and 1983–2018 (spring), the longest uninterrupted record of autumn counts for this set of watch sites. Variables were evaluated for normality with the Shapiro-Wilk test. Total seasonal counts at each site were log₁₀-transformed, as were the three-site sum (Interior Sum) and each two-site sum of the counts from the three interior autumn sites for the whole season. The same transformation was performed on the four monthly Interior Sums and three monthly counts for Cape May. After we inspected plots of log₁₀-transformed counts against year, ordinary least-squares linear regression was used for temporal detrending. For every data set, the linear model had a lower small-sample Akaike's Information Criterion (AIC_C) score than the second-order polynomial model did. All log₁₀-transformed counts and their sums were detrended, whether the linear trend model was statistically significant or not. Temporal autocorrelation and partial autocorrelation analyses were performed on the regression residuals. Standard deviations of the regression residuals were computed, representing amplitudes of count fluctuations (Turchin 2003). Non-transformed (raw) counts for the Interior Sum, Cape May, and Derby Hill were linearly detrended by ordinary least-squares regression; in every case the AIC score for the linear model was lower than it was C for the second-order polynomial model. We used the Fisher-Pearson g₁ as the skewness of the 1-yr rate of change in the detrended raw counts for the Interior Sum and Cape May, r_t = ln (N_t/N_t–1), where N_t is the detrended raw count rounded to the nearest integer for year t. The g₁ was also calculated as the skewness for the residuals of the Interior Sum and Cape May log₁₀-transformed counts.

We used the North Atlantic Oscillation Index (NAOI) based on the linear principal components for station-based pressures (National Center for Atmospheric Research Staff 2018) for the September–November and March–May periods ( Supplemental Materials (38_rapt-54-01-04_s03.pdf)). Residuals of log₁₀-transformed harrier counts were tested for association with the NAOI by using Pearson's r. Annual September–November NAOI values were detrended by ordinary least-squares linear regression before temporal autocorrelation analysis was performed on the residuals to test for multiyear atmospheric cycles.

Differences in harrier counts among years could be influenced by uncontrolled variation in the observation process or conditions for observation. Such variations may be evidenced by their influence on the counts of other species. As we had for harriers, we tested residuals of log₁₀-transformed counts for 14 other abundant diurnal raptor species summed from the same three interior autumn sites and same calendar intervals ( Supplemental Materials (38_rapt-54-01-04_s03.pdf)) for temporal autocorrelation. These transformed sums of counts were detrended by ordinary least-squares regression models that were linear, second-order polynomial, third-order polynomial, or fourth-order polynomial, using the model that had the lowest AIC_C score, regardless of statistical significance. Pearson's r was used to test for association among the several harrier count residuals, as well as between the residuals for the harrier Interior Sum and each of the other diurnal raptor species. All statistical procedures ( Supplemental Materials (38_rapt-54-01-04_s03.pdf)) were carried out with Statistix 9 (Analytical Software, Tallahassee, FL).

Results

Autocorrelation coefficients (r_k) for residuals (detrended log₁₀) of autumn harrier counts were significantly (P < 0.05) positive at the k = 4 yr lag for Hawk Ridge, Holiday Beach, Hawk Mountain, the Interior Sum of all three sites, Hawk Ridge + Holiday Beach, Holiday Beach + Hawk Mountain, and Hawk Ridge + Hawk Mountain (Table 2). No other lags had a significant positive autocorrelation for the Interior Sum (Fig. 2). Consecutive peaks or troughs in the detrended raw Interior Sum were 4 or 5 yr apart (Fig. 3). Holiday Beach almost always had the highest fraction, 36% to 69% (mean = 54%) of the raw Interior Sum in any one year. The strongest periodicity of the Interior Sum residuals was for counts made in October and September (Table 2), with nonsignificant 4-yr autocorrelations for August and November (Fig. 4). A mean of 47% of the whole-season raw Interior Sum was counted in September and 37% in October. Amplitudes of residuals ranged between 0.113 and 0.237 across the data sets (Table 3).

Table 2.

Statistically significant (P < 0.05) values of the temporal autocorrelation coefficient for Northern Harriers at the 4-year lag (r_4-yr). Each data set represents the full observation interval unless otherwise indicated.

Table 3.

Standard deviations (amplitudes) of the residuals for all Northern Harrier data sets. Skewness of the 1-yr rate of change (r_t) of the detrended raw count for Cape May and Interior Sum. Skewness of the residuals for Cape May and Interior Sum.

Figure 2.

Temporal autocorrelograms for residuals of log₁₀ autumn-migrating Northern Harriers (Interior Sum, IS; Cape May, CM), residuals of the autumn North Atlantic Oscillation Index (NAOI), and partial autocorrelogram for log₁₀ Interior Sum (IS Partial), for years 1983–2017. Markers filled black represent r_k or partial-r_k different from 0 at P < 0.05, where k is the lag (yr). IS and IS Partial had identical r_k values of –0.480 at the k = 2-yr lag, with the IS marker left unfilled so that both can be seen.

Figure 3.

Detrended raw counts of migrating Northern Harriers for years 1983–2017: Interior Sum (IS), Cape May (CM), and Derby Hill (DH).

Figure 4.

Temporal autocorrelograms for residuals of log₁₀ monthly Interior Sum counts of migrating Northern Harriers for years 1983–2017: August (A), September (S), October (O), November (N). Solid black markers represent r_k different from 0 at P < 0.05, where k is the lag (yr).

The partial autocorrelogram of the residuals of the full-season Interior Sum had significant partial-r_k values at lags of k = 2, 3, and 4 yr (Fig. 2). This indicates a higher-order oscillation (Turchin 2003), with multiyear delayed feedback underlying the harrier population changes. The Interior Sum's autocorrelation and partial autocorrelation patterns (Fig. 2) allowed us to reject a hypothesis that the Interior Sum's trajectory represents merely noncyclic changes in the number of migrating harriers and errors in estimation. Collectively, these point to a substantial element of periodicity in the Interior Sum, consistent with a hypothesis of a regionally cyclic harrier population size. Of the 14 other diurnal raptor species we analyzed for the same interior watch sites, the r_4-yr of the residual of their Interior Sum was significantly positive for only Sharp-shinned Hawk (Accipiter striatus), Rough-legged Hawk (Buteo lagopus), and Turkey Vulture (Cathartes aura), albeit weaker than for harrier (r_4-yr = 0.401, 0.367, and 0.353, respectively;  Supplemental Materials (38_rapt-54-01-04_s03.pdf)). Harrier residuals were positively correlated with those from the same year for Sharp-shinned Hawk (r = 0.662, P < 0.001) and Turkey Vulture (r = 0.344, P = 0.043), but negatively correlated with those for Rough-legged Hawk (r = –0.414, P = 0.013).

The waves of the Interior Sum's harrier cycle were asymmetric, with decreases steeper than increases. Their skewness (Table 3) was moderately (Bulmer 1979) negative, with the mean of the 12 1-yr decreases = –0.489 and that of the 22 1-yr increases = 0.290. The skewness of the residuals of the Interior Sum (Table 3) indicated at most a rather weak tendency for the troughs to be sharper than the peaks.

In contrast, 1983–2017 autumn harrier counts at the Atlantic coastal Cape May (Fig. 3) had residuals that show no strong evidence of a 4-yr cycle, but had a significant negative r_k at the k = 5 yr lag (Fig. 2). Cape May's mean of its 17 1-yr decreases = –0.257 and that of its 17 1-yr increases = 0.286, producing a skewness (Table 3) that was not strong evidence of asymmetry. Despite the lack of a 4-yr cycle at Cape May, residuals of the Cape May and Interior Sum harrier counts in the same year were positively correlated (r = 0.484, P = 0.003). Those residuals were negatively correlated with the North Atlantic Oscillation Index (NAOI) for the full September–November period of the same year for both Cape May (r=–0.427, P=0.011) and the Interior Sum (r= –0.334, P = 0.050). For both sets of counts, the negative correlation appeared to result mainly from a lack of years in which both the seasonal NAOI and the residual of the count were strongly negative (Fig. 5). When residuals for single months were compared to the same month's NAOI, only September at Cape May had a significant correlation (Table 4). No temporal autocorrelations for the autumn NAOI were significant (Fig. 2). The correlations of the NAOI with residuals of Interior Sum Sharp-shinned Hawk (r = 0.242), Turkey Vulture (r = 0.061), and Rough-legged Hawk (r=0.056) were not significant.

Figure 5.

Residuals of log₁₀ autumn-migrating Northern Harriers for Cape May (CM) and Interior Sum (IS) plotted against the autumn North Atlantic Oscillation Index for the same year, 1983–2017.

Table 4.

Correlations (r) between the residuals of log₁₀ single-month Northern Harrier counts and the NAOI values for the same month, 1983–2017, with P the probability that an r value at least as extreme would differ from 0 by chance.

Residuals of log₁₀ counts of north-bound harriers in the previous spring at Derby Hill (Fig. 3 shows detrended raw counts) had no correlation with the autumn count residuals of the Interior Sum (r = 0.036) or Cape May (r = –0.110), supporting an inference that the Interior Sum's cycle in autumn migrants resulted from a cycle in harrier productivity during the immediately previous breeding season. Residuals of Derby Hill spring counts in 1984–2018 were positively correlated with the previous autumn's residuals at Cape May (r = 0.360, P = 0.034) and those of the Interior Sum (r = 0.350, P = 0.039). The residuals of the spring counts did not correlate with the March–May NAOI for 1983–2018 (r = –0.031).

Discussion

Each of the three interior autumn sites had autocorrelation at the 4-yr lag (r_4-yr) that is smaller for its own count data alone than it is for the two- and three-site sums in which it is included (Table 2). This may be taken as evidence that each site's count represents a spatial subsample of the same regional cycle. The highest r_4-yr is the one for the Interior Sum and the lowest r_4-yr among the pairwise sums is the one for the sites with the greatest geographic separation, Hawk Mountain and Hawk Ridge. Some of the harriers counted at each of the three interior autumn sites would be expected to have come from populations that were more synchronous in productivity, the remainder from less-synchronous or noncyclic ones. The exceptionally high autocorrelation coefficient for Interior Sum October counts (Table 2) suggests that tallies of harriers from more-synchronous sites are centered in that month. Harriers migrating from more-synchronous summer populations would be expected to make up their greatest fractions of the monthly Interior Sum in October of peak years, with harriers from poorly synchronous or noncyclic summer populations tending to make up their greatest fractions of the October Interior Sum in trough years (Fig. 4).

An alternative hypothesis is that the cycle in autumn tallies represents cyclic interannual changes in region-wide spatial patterns of migration rather than in summer productivity. In this scenario the number of regional migrants would not cycle but their proximity to the three interior watch sites would. Regional atmospheric dynamics could affect the migration routes taken, but any aspect of those dynamics that is captured by the season-long NAOI does not show strong interannual periodicity (Fig. 2). This alternative hypothesis could be tested by examining the structure of migrant harrier populations. Adult males can be distinguished from adult

female and juvenile harriers when observed in flight (Godfrey 1986), and would be expected to compose a smaller fraction of the autumn tally in years when there is higher summer production of juveniles, which would be the peak years if productivity is the basis for the cycle. This alternative hypothesis would also need to explain the positive correlation of autumn counts with north-bound harrier counts at Derby Hill in the following spring. Because few of the other 14 species studied show a 4-yr cycle in the count, we should be able to reject a 4-yr cycle in general conditions for observation or in observer effort/skill as an external cause for the harrier cycle. Sharp-shinned Hawks also had a significantly positive r_2-yr ( Supplemental Materials (38_rapt-54-01-04_s03.pdf)), which could mean that the significant r_4-yr represents an echo of a 2-yr cycle. We offer no explanation for a 4-yr cycle for Turkey Vultures. The Rough-legged Hawk cycle is strongly out of phase with the harrier cycle, ruling out an external cause common to both.

The negatively skewed frequency distribution of r_t for the Interior Sum (Table 3) is more likely to be seen when fluctuations result more from demographic processes such as a productivity cycle than from exogenously driven changes in the counts; the nonsignificant skewness of the residuals of the Interior Sum (Table 3) indicates that the residuals' trajectory does not have the especially sharp peaks that are often expected for a predator in a cycle (Turchin 2003). If anything, there is a weak tendency toward sharper troughs. Harriers nest on the ground and have increased depredation rates on their eggs, nestlings, and incubating females—by other species that prey on voles—when vole abundance is low (Smith et al. 2011). Regionally synchronous decreased availability of preferred prey (voles) and concomitant increased losses to predation may deepen the troughs and make skewness in r_t more negative.

Consumer-resource theory (Murdoch et al. 2003) suggests that the harrier's ≥ 1 yr age at first reproduction with repeated reproduction by some adults in subsequent years (Smith et al. 2011) makes it unlikely that a 4-yr cycle results from processes internal to its population. It is much more probable that interior autumn harrier counts oscillate in response to a large subset of the regional vole population that has an internal 4-yr synchronous cycle, which is consistent with theory (Murdoch et al. 2003) because voles become reproductive at just a few months of age (Banfield 1974). The Interior Sum's autocorrelogram (Fig. 2) has a dominant period of 4 yr, which matches a commonly observed period for vole population cycles (Turchin 2003). The amplitude of the harrier cycle (Table 3) is one to two orders of magnitude less than the amplitudes reported for various species of voles at high latitudes (Turchin 2003), as would be expected for a vole-responsive predator.

Where are Summer Populations More Synchronous? The breeding areas for harriers counted at Cape May must differ substantially from those counted at the Interior Sum sites, given their contrasting autocorrelation patterns. Research on Sharp-shinned Hawks found that many of the birds trapped at Cape May during autumn migration were banded at breeding sites near the Atlantic and many of the autumn migrants trapped at more-interior sites were banded at more-interior breeding sites (Goodrich and Smith 2008). Knowing the locations of the breeding grounds of the harriers counted at these autumn migration sites would enable research on population dynamics of voles and their predators to be focused there. However, other than Bildstein's (2006) assertion that many harriers seen in autumn migration at Cape May had spent the summer in Atlantic coastal wetlands, broad-scale information connecting migrant harriers to their nesting areas is not available. Most of the high-latitude summer area for harriers in eastern Canada (Bird Studies Canada 2006) is in remote wetlands difficult for banders to access.

Future research should consider analysis of trace elements (Bortolotti 2010) or stable isotopes (Haché et al. 2012) to discover geographic connections between migrant harriers and their breeding areas. Notably, trapping of some migrant raptors to obtain samples of feathers or other tissue already takes place at some of these autumn observation sites (Bildstein 2006); future success with this or other methods would depend on continued occurrence of the regional cycle. The stronger periodicity in October at the interior sites (Fig. 4) could sharpen the power to narrow the geographic origin of cycling harriers if chemical composition data were segregated by month for trapped hatch-year birds.

Comparability of Count Calendar Intervals. Counts began 15 August at Hawk Ridge and Hawk Mountain, but not until 1 September at the other two autumn sites (Table 1). Hawk Ridge's August harrier counts were 2%–38% of its August–October total in the same year, with a mean of 11%. Hawk Mountain's August harrier counts were 1%–10% of its August–November total, with a mean of 5%. Counts for all of November at Hawk Ridge began in 1993, and would have added 0%–9% to their year's harrier total, with a mean of only 2%. The mean date of autumn harrier passage for eastern North America trended later by about 2.2 d from 1985 to 2012 (Therrien et al. 2017). For Derby Hill, the mean date of spring harrier passage trended earlier by about 5.4 d from 1980 to 2011, with much less change in the dates for the first and third quartiles of its annual tally (Sullivan et al. 2016). These variations are small enough that they do not lead us to question the general patterns that we found.

Further Implications. Functionally, the interior harrier cycle could impart cyclic interannual changes in predation activity (Bauer and Hoye 2014) during autumn migration stops and in parts of the harrier's wintering region. The influence of autumn migrant numbers on the abundance of wintering harriers apparently persists through the following spring's migration, as evidenced by the positive correlations between the Cape May and Interior Sum counts in the autumn with those at Derby Hill in the following spring. Individual harriers from the northern part of the summer range tend to migrate to the southern part of the species' winter range (Bildstein 2006). Local winter predation pressure by harriers will further depend on the degree of geographic concentration of harriers in response to winter severity.

Autumn counts of harriers at Hawk Mountain were unrelated to the local passage of atmospheric cold fronts (Allen et al. 1996), which may make the statistical relationship between the season-long harrier residuals and NAOI (Fig. 5) surprising. The atmospheric circulation represented by the NAOI in autumn is detectable as far west as the region between James Bay and the Great Lakes, with anomalous near-surface airflow from the south in the eastern United States when the NAOI is positive (Hurrell and Deser 2009). If there is a more continuous or forceful effect of the NAO on overland circulation closer to the Atlantic, that would be consistent with the NAOI's correlation with the Cape May counts being stronger than its one with the Interior Sum. Evaluating a functional connection between the NAO and harrier movements would require scrutiny of a possible relationship between the September–November NAOI and any interannual differences in autumn weather relevant to the migrants that potentially reach the observation sites. Although the associations between the season-long counts and September–November NAOI are suggestive, the correlations between the counts and NAOIs for individual months are weak except for Cape May in September (Table 4). It is possible that the NAO influences the season-long count at Cape May but correlates with the season-long Interior Sum mostly because these two counts are positively correlated for other reasons. If a scenario that the NAO does influence the season-long count is entertained, detrending the log₁₀ Interior Sum by the September–November NAOI as well as by year ( Supplemental Materials (38_rapt-54-01-04_s03.pdf)) reduces its r_4-yr from 0.697 to 0.572 (both P < 0.05). If the NAO is found to be influential, attempts to make inferences about regional harrier populations from long-term trends of migration counts in northeastern North America may gain some power by adjusting for the NAOI.

More conclusively, the autumn harrier cycle should be considered when inferences about temporal changes in regional harrier abundance are sought from interior migration data. The twofold to fourfold difference between peak and trough counts within each wave of the Interior Sum (Fig. 3) is larger than statistical adjustments of migration counts that have been made on the basis of interannual variation in hours of observation (Farmer et al. 2007). Our detrending of counts was only for the purpose of exploring possible periodicity, and trends expressed in our regressions ( Supplemental Materials (38_rapt-54-01-04_s03.pdf)) were not adjusted for variation in observation effort. Each year-to-year change in the Interior Sum resulted more from the phase shift within a wave than it did from longer-term tendencies. The clearest investigations of simple multiyear trends would use first and last years that both represent, retrospectively, the same position in a wave. Such an approach could be taken for spans as short as 4 or 5 yr.

Supplemental Materials (available online)

(1) Raw counts of Northern Harriers by year and month for each count site. Values of the North Atlantic Oscillation Index by year and month (Excel file). (2) Raw counts of 14 diurnal raptor species, other than Northern Harrier, by year for Interior Sum of Hawk Mountain, Hawk Ridge, and Holiday Beach (Excel file). (3) Results of statistical analyses for each diurnal raptor species and for North Atlantic Oscillation Index (Word file).