We conducted our study in a forest landscape located on Navarino Island in the southernmost region of Chile (Fig. 1). The study landscape was covered by southern beech forest of Nothofagus betuloides, N. pumilio, and N. antartica (Fig. 1). Open habitats also were present in this landscape and included patches of shrublands, wetlands, peat bogs, meadows, and ponds, with the latter two being produced by the introduced beaver (Castor canadensis) (Fig. 1; [41]). The cover of N. pumilio, N. betuloides, and N. antarctica was 29.2%, 29.5% and 6.4%, respectively (Table 1). Although forest stands in old-growth stage of succession covered 45.3% of the land, forest disturbances (e.g., logging and fires) have resulted in some second-growth stands (20.2%). Shrublands and meadows (upland open areas) covered 10.2% of the study area and peat bogs and pond (lowland open areas) represented 6.5% of the landscape (Table 1).
Table 1
Cover (%) of each habitat type in the study area where the movement of Magellanic Woodpecker was studied. For PSRI, the mean and SE value are shown.
Habitat type | Total cover (%) |
N. antarctica | 6.38 |
N. betuloides | 29.23 |
N. pumilio | 29.52 |
Open upland habitats | 10.21 |
Open lowland habitats | 6.49 |
PSRI* | -0.73 (0.31) |
Old-growth forest | 45.34 |
Second-growth forest | 20.22 |
*mean (SD) | |
Movement data
We acquired GPS locations from 24 tagged male Magellanic Woodpeckers [39, 28] using ATS G10 UltraLITE GPS Logger (Advanced Telemetry Systems, Inc.) devices, each attached to a very high-frequency transmitter (ATS model A2440, 2.3 g) for later recovery. GPS devices were placed on the back of adult male woodpeckers using a small amount of epoxy to six feathers. We chose adult males because males guide family groups by eliciting a dominant social behavior while moving across forest habitat [42, 43]. The locations of woodpeckers were recorded every 5 min between 08:00 to 16:00 and during the 2014–2015 post-reproductive period (Jan to Mar).
We estimated the accuracy (measurement error) of GPS locations as 12.9 ± 2.9 m (mean ± SE), which was obtained by quantifying the Euclidean distances (m) between 12 different GPS measurements and the actual position of a reference point identified on an imagery-based map layer [28]. We calculated the overall speed for each woodpecker by dividing the total traveled distance (i.e. the sum of the distance between GPS relocations) by the total time of a given burst of continuous relocations.
First-passage time
We examined the role of habitat heterogeneity on the movement of woodpeckers by estimating the first-passage time (FPT) from GPS location data. FPT can be used as an indicator of foraging behavior along animal's tracks, with longer FPT areas being interpreted as evidence of ARS at ecologically proper spatial scales [31, 44, 45]. The first-passage time is defined as the time spent by an individual in circles of radius r centered on subsequent GPS positions along each trajectory (i.e., a movement path consisting of a sequence of GPS locations within the home range; Fig. 2A;[46]). As r increases, longer trajectory sections will be included in the circle [31]. We established the proper spatial scale of FPT analysis by searching the value of r that maximizes the relative variance [\(S\left(r\right)\)] of the FPT because the ability of the FPT to detect area restricted search (ARS) increases as the variance takes maximum values [31]. The ARS is a behavior characterized by slow and tortuous movements typically displayed by woodpeckers when selecting trees for foraging [2].
To estimate the FPT, we used the fpt function of the R package adehabitatLT. The last six woodpecker’s data were excluded from further analysis because it was not possible to define a regular trajectory due to the different time lag between relocations (Fig. S1). The first-passage time method is designed for trajectories with three or more relocations, so the trajectory data with less than three observations were eliminated. Because observations of six woodpeckers were discarded from the analysis, we estimated the first-passage time (FPT) from the trajectories of 18 woodpeckers. We maximized \(S\left(r\right)\) for each trajectory and individual by estimating FPT over 50 different radii (r) in a range from 12 to 250 m corresponding, respectively, to the GPS accuracy and a quarter of the calculated net distance displacement of Magellanic Woodpeckers [32]. We determined the proper FPT scale, r value at which \(S\left(r\right)\) reached its maximum value, by examining plots of \(S\left(r\right)\) against r. From \(S\left(r\right)\) plot examination, a set of 36 trajectories were selected out of 62 trajectories (Fig. S2). When we did not observe a maximum value of \(S\left(r\right)\), we assumed that the path traced by a woodpecker was random and did not represent a movement pattern including different behavioral modes, as shown by the eighth trajectory of woodpecker 3 (Fig. S2). However, we considered proper FPT scales for trajectories where a local maximum was observed, for example, the first trajectory of woodpecker 4 (Fig. S2). Based on the proper FPT scales, we defined two spatial scales to quantify properties of the habitat used by woodpeckers (see below), including site (habitat-patch) scale, defined by the area of each FPT circle (Fig. 2A), and home-range scale, defined by the area comprising the union of all the FPT circles along the trajectory (Fig. 2B).
Habitat variables
We quantified habitat variables at the site and home-range scales using a high-resolution (0.50 m) multispectral image from the WorldView-2 sensor (2014). We created geographic information system (GIS) layers of the main categories of habitat types and tree species in the study site as well as a remote sensing index of tree senescence (Table 1). We used digital supervised classification and a Bayesian maximum likelihood algorithm carried out by [39] to classify habitat types based on age and composition of forest stands as well as the cover of different habitat types. We categorized habitats as old-growth forest, second-growth forest, open upland (shrub and meadows), and open lowland (peatlands and beaver ponds). We also determined the composition of dominant tree species by quantifying the cover of N. Antarctica, N. betuloides, and N. pumilio. Tree senescence was estimated with the remote sensing-based Plant Senescence Reflectance Index (PSRI), which distinguishes between tree decay states based on the spectral carotenoid/chlorophyll ratio with increasing values for increasing tree decay [39]. We used an image segmentation algorithm to identify individual trees, estimate their PSRI and classify them by species of tree [37, 39]. The PSRI values of the subpolar forest range from − 2.7 to 0.4 (Table 1, [39]).
We used the GIS maps of habitat type, tree composition and PSRI to derive the habitat variables used in statistical analysis. First, for each FPT circle, we estimated the percentage of old-growth forest, second-growth forest, open upland habitats and open lowland habitats. Second, we estimated the percentage cover of the dominant Nothofagus species. Third, we averaged PSRI values over all cells in each FPT circle and considered those values as estimates of tree senescence at the site level (Fig. 2). Fourth, PSRI values were averaged over all FPT circles, considering those values as estimates of tree senescence at the home-range level (Fig. 2). Tree senescence estimated at the site and home-range scales were not collinear (see below; Table S1), hence those variables were included in the same statistical models. All tree species located at the site and home-range levels (i.e. all trees across the trajectories) were interpreted as the foraging habitat quality at the home-range scale [37].
Statistical modelling
We used a methodological framework based on Linear Mixed Effect models (LMM) to determine the effects of habitat variables on FPT (Fig. 3). This approach was intended to reduce spatial autocorrelation arising from consecutive circles whose areas tend to be overlapped along trajectories ([31]; Fig. 3A). To analyze independent FPT data, we performed a randomization procedure by randomly selecting subsets of not-overlapping circles for each trajectory (Fig. 3). This procedure was repeated 1000 times, which resulted in 1000 sets of trajectories, each contained independent data later used in LME analyses (Fig. 3). LME were fitted to different datasets, which precluded comparing models with a criterion derived from likelihood functions, such as the Akaike Information Criteria (AIC; [47]). Thus, for each of the 1000 sampled datasets, we used the RMark R package to compute model-averaged coefficients based on the AIC weights [48]. The AICc weight quantifies the probability that a given model is the best among a set of candidate models [48, 49]. The distributions of model-averaged coefficients were assumed to represent the effects of predictors. The mean and 95% confidence intervals of those distributions were used to interpret the significance of the coefficients. Model averaging was carried out on a set of 93 models including all possible combination of variables derived from a global model. The global model was built to assess the independent effects of habitat variables. Thus, we checked collinearity with the variance inflation factor (VIF; Table S1), where a VIF > 10 indicates lack of orthogonality between variables [50]. From this analysis, we detected that the percent cover of Nothofagus species and that of old-growth forest were collinear, so we discarded the latter from analyses (Table S1). Therefore, the global model included as predictors the cover percentage of dominant Nothofagus species, second-growth forest, open upland habitats and open lowland habitats, in addition to the PSRI values quantified at site and home range levels. The dependent variable (FPT) was divided by the area of each circle to obtain the time woodpeckers spent in an area of similar size (s/m2), and thus to allow the comparison between trajectories with circles of different radii. The trajectory and individual were included as random effects in LMM, with analysis being performed using the lme4 R package [51].