Study area

The Juan Fernández Archipelago is a volcanic group of islands and islets located c. 660 km off the coast of central Chile (Fig. 1). The three main islands of the archipelago – Robison Crusoe (48.0 km²), Santa Clara (2.2 km²) and Alejandro Selkirk (49.5 km²) – are characterized by a rugged topography and a subtropical climate (Haberle Reference Haberle, Gillespie and Clague2009). Forested areas occur in lowland areas and gullies, whereas the upper elevations are dominated by bush and grasslands, interspersed with unvegetated zones (Stuessy et al. Reference Stuessy, Swenson, Crawford, Anderson and Mario-Silva1998, Bernardello et al. Reference Bernardello, Anderson, Stuessy and Crawford2006, Haberle Reference Haberle, Gillespie and Clague2009). The proximity of the archipelago to the Humboldt Current means that the ocean surrounding the islands is relatively productive, which in turn translates into high terrestrial and marine biodiversity (Hahn et al. Reference Hahn, Römer, Vergara and Walter2009). It is not surprising, then, that the Juan Fernández Archipelago is home to important breeding colonies of seabirds and the Juan Fernández fur seal (Arctocephalus philippii) (Colwell Reference Colwell1989, Bourne et al. Reference Bourne, de Brooke, Clark and Stone1992, Hahn et al. Reference Hahn, Römer, Vergara and Walter2009). The natural values of the Juan Fernández Archipelago are recognized and protected as a Chilean National Park (www.conaf.cl/parques/parque-nacional-archipielago-juan-fernandez).

The Juan Fernández Archipelago is sparsely populated by humans, with a total population of 926 people in 2017 (https://reportescomunales.bcn.cl/2017/index.php/Juan_Fernández). The entire population lives in San Juan Bautista on Robinson Crusoe Island (the main township in the islands) and seasonally on Alejandro Selkirk, while Santa Clara Island is uninhabited.

Monitoring of pink-footed shearwater breeding colonies

We monitored a representative sample of breeding burrows in six colonies on two islands in the Juan Fernández Archipelago (Fig. 1 & Table 1): Robinson Crusoe Island (three colonies) and Santa Clara Island (three colonies). The ruggedness of the terrain and the difficulty of access to the breeding colonies imposed logistical and safety constraints to monitoring the breeding colonies. Moreover, many breeding burrows were located in soft volcanic soil, making them prone to being altered or destroyed by human activity. Therefore, we took care to ensure that our monitoring activities did not disturb or affect the breeding burrows. Our monitoring regime included reproductive samples chosen non-randomly to maximize spatial representativeness within colonies while ensuring accessibility for the researchers and a low risk of disturbance to the breeding burrows. In each colony, we monitored 30–100 burrows each year (Table 1). We excluded from the sample any burrows that were too long for the infrared camera to confirm nest contents (the 3-m-long camera was sufficient to see the back of the majority of burrows). In the smaller colonies on Santa Clara Island, all viable burrows were included.

Table 1.

Characteristics of the six pink-footed shearwater colonies monitored in the Juan Fernández Archipelago. Number of burrows monitored each year and summary statistics of the burrow occupancy and hatching success per year are shown (mean ± standard deviation, and range; excluding missing data for hatching success).

The breeding burrows were monitored using a miniature camera with infrared illumination through two stages of the annual reproductive season: (1) burrow occupancy – the viable burrow being occupied by a breeding pair with an egg; and (2) hatching success – transitioning from the presence of an egg in a burrow to the presence of a chick. All of the burrows were monitored during the first stage (burrow occupancy), but not all of them were monitored during the second stage (hatching success) within the same reproductive season. Unoccupied burrows were not monitored for hatching success, and we optimized our monitoring efforts by selecting a sample of occupied burrows to be subsequently surveyed for hatching success. Burrow occupancy was determined by two checks during the incubation period between late December and mid-January. Burrows were considered occupied if an egg was seen once or if an adult was seen in the burrow during both checks (it was not always possible to see the egg under the incubating bird). A previous study determined that >99% of burrows with adults observed during both incubation-stage checks contained an egg (Hodum et al. unpublished data). We determined hatching success by checking burrows once during mid- to late February, following the conclusion of the hatching period. Presence of a chick during this check indicated successful hatching, while absence of a chick indicated failure to hatch. Burrows with no egg or an unattended egg during the hatching check were considered to have failed to hatch a chick. Similarly to other procellarid seabirds, pink-footed shearwaters are not known to re-lay eggs after the first attempt, so failure to hatch indicated reproductive failure for a breeding pair for the entire year.

Three colonies were monitored on Robinson Crusoe Island (see Fig. 1 & Table 1 for further details): (1) Piedra Agujereada had a total of 87 burrows monitored from 2014 to 2016 (total 64 burrows with information on hatching success, 73.6%); (2) Reserva had a total of 44 burrows monitored from 2005 to 2007 (30 burrows with information on hatching success, 68.2%); and (3) Tierras Blancas had 67 burrows monitored from 2015 to 2016 (all burrows with information on hatching success). Another three colonies were monitored on Santa Clara Island (see Fig. 1 & Table 1 for further details): (1) Cerro Alto had a total of 376 burrows monitored from 2002 to 2016 (221 with information on hatching success, 58.8%); (2) Refugio had a total of 760 burrows monitored from 2002 to 2016 (386 burrows without information on hatching success, 50.8%); and (3) Volcán had a total of 694 burrows monitored from 2002 to 2016 (360 burrows with information on hatching success, 51.9%). Note that the number and spatial location of the breeding burrows monitored were not necessarily constant across the whole study period.

The spatial and temporal autocorrelation induced by our monitoring design was accommodated in our modelling framework that accounts for the year- and island-specific effects (see details below). The between- and within-island variation in the presence of alien mammals allowed us to develop and fit a statistical model to evaluate the impacts of different alien mammal assemblages on the burrow occupancy and hatching success of the pink-footed shearwater (Table 2). The presence of each of the seven alien mammals in each of the six pink-footed shearwater colonies was obtained from the literature (Hahn & Römer Reference Hahn and Römer2002, Croxall et al. Reference Croxall, Butchart, Lascelles, Stattersfield, Sullivan, Symes and Taylor2012, Carle et al. Reference Carle, Beck, Colodro and Hodum2016) and from our personal observations during the fieldwork. The burrow monitoring data are publicly available from FigShare: 10.6084/m9.figshare.12017814.

Table 2.

Calculation of the probabilities of burrow occupancy and hatching success as a function of the five alien mammal assemblages in the six colonies monitored.

p _–rabbit = the probability of occupancy (po) or hatching success (ph) given rabbits are absent (i.e., the probability of burrow occupancy or hatching success when only rabbits are absent); p _–cow = the probability of occupancy (po) or hatching success (ph) given cows are absent (i.e., the probability of burrow occupancy or hatching success when only cows are absent); p _–coati = the probability of occupancy (po) or hatching success (ph) given southern coatis are absent (i.e., the probability of burrow occupancy or hatching success when only southern coatis are absent); p _–others = the probability of burrow occupancy (po) or hatching success (ph) given alien mammals other than cows, coatis and rabbits are absent, but cows and rabbits are present (i.e., the probability of burrow occupancy or hatching success when cows, coatis and rabbits are absent, but other alien mammals are present).

Statistical analysis

We developed a model to separate the effects of alien species assemblages (Table 2), the island where the monitored burrows were located (Santa Clara or Robinson Crusoe) and the monitoring year. We assumed that these three factors affected the observed burrow occupancy and hatching success independently. This means that the observed burrow occupancy and hatching success in each monitored burrow were assumed to be the outcomes of the combined, but independent, effects of the three factors. This yielded a model whose components can be interpreted as the probabilities of burrow occupancy and hatching success directly attributable to the alien mammal assemblages, the island and the monitoring year alone, regardless of the effects of the other factors. Therefore, to obtain the probability of burrow occupancy and hatching success in any given breeding burrow on any of the two islands, it was necessary to use the three independent probability estimates together and not separately (see model description below).

We modelled the transitions between the two subsequent stages, defining burrow occupancy and hatching success as a function of five alien species assemblages, monitoring year and island through a Bayesian hierarchical model, whereby one stage depends on the state of the previous stage (Hobbs & Hooten Reference Hobbs and Hooten2015, Kéry & Royle Reference Kéry and Royle2016):

(1) $${O_{j,z,y}}\sim Bernoulli\left( {p{o_{j,z,y}}{V_{j,z,y}}} \right),$$

(2) $$p{o_{j,z,y}} = 1\!-\!\left( {\left( {1-p{o_z}} \right)\left( {1-p{o_{sp}}} \right)\left( {1\!-\!p{o_y}} \right)} \right),$$

(3) $${\rm{logit}}\left( {p{o_y}} \right) = \beta oy\,{\rm{year}},$$

where O _j,z,y ϵ {0, 1} is the occupancy of breeding burrow j of a colony on island z during year y, V _j,z,y is the viability of breeding burrow j of a colony on island z during year y, po _j,z,y is the probability of occupancy of breeding burrow j of a colony on island z during year y, po _z is the island-specific probability of burrow occupancy (one parameter for each island), po _sp is the alien mammal assemblage-specific probability of burrow occupancy (see Table 2 for equations to derive these probabilities), po _y is the year-specific probability of burrow occupancy and βoy is the effect, on the logit scale, of the monitoring year on po _y. The hatching success was modelled conditional on burrow occupancy (a burrow needs to be occupied for hatching to occur):

(4) $${H_{j,z,y}}\sim Bernoulli\left( {p{h_{j,z,y}}\,{O_{j,z,y}}} \right),$$

(5) $$p{h_{j,z,y}} = 1\!-\!\left( {\left( {1-p{h_z}} \right)\,\left( {1\!-\!p{h_{sp}}} \right)\,\left( {1\!-\!p{h_y}} \right)} \right),$$

(6) $${\rm{logit}}\left( {p{h_y}} \right) = \beta hy\,{\rm{year}},$$

where H _j,z,y ϵ {0, 1} is the hatching success in burrow j of a colony on island z during year y, ph _j,z,y is the probability of hatching success in burrow j of a colony on island z during year y, ph _z is the island-specific probability of hatching success (one parameter for each island), ph _sp is alien mammal assemblage-specific probability of hatching success (see Table 2 for equations to derive these probabilities), ph _y is the year-specific probability of hatching success and βhy is the effect, on the logit scale, of the monitoring year on ph _y.

Note that the effect of the monitoring year was measured by equating the first year when any monitoring occurred (2002) to zero, and subsequent years were re-coded accordingly as one (2003), two (2004) and so on. Therefore, the year 2002 provided the baseline for assessing the temporal trends. The year 2002 was also the reference for those colonies for which the first monitoring took place later than 2002. For instance, this means that the first year that Piedra Agujereada was monitored (2014) was re-coded as 12 in our analyses. By following this approach, together with the cross-validation procedure described below, we reduced the potential biases in the modelling induced by the unbalanced sample size corresponding to each year and colony (Table 1). Moreover, we did not estimate year-specific parameters (i.e., one parameter for each year), but instead a trend across years and colonies. This will further minimize any potential effects of the unequal sample sizes.

Our hierarchical model was fitted to the breeding burrow monitoring data using Bayesian Markov chain Monte Carlo methods implemented in the package NIMBLE (de Valpine et al. Reference de Valpine, Turek, Paciorek, Anderson-Bergman, Lang and Bodik2017) as interfaced to the R statistical environment (R Development Core Team 2019). Given the imbalance in the monitoring data (different sample sizes for the two islands and the five invasive species assemblages) and the missing hatching success data, we used a 10-fold cross-validation approach to fit and validate the model (Hooten & Hobbs Reference Hooten and Hobbs2015, Roberts et al. Reference Roberts, Bahn, Ciuti, Boyce, Elith and Guillera-Arroita2017, Conn et al. Reference Conn, Johnson, Williams, Melin and Hooten2018). In each fold, we randomly selected 30 burrow occupancy observations of each alien mammal assemblage (150 observations in total) and randomly assigned 70% of those observations to a within-sample data set (a training data set used to fit the model; 105 observations) and 30% to an out-of-sample data set (a test data set used to validate the model; 45 observations). This procedure produced a total of 150 burrow occupancy observations used in each fold, such that for each alien mammal assemblage there were 21 observations in the within-sample data set and nine observations in the out-of-sample data set.

We ran our Bayesian hierarchical model using one chain with 42 000 iterations and no thinning in each of the 10 folds. We discarded the first 2000 iterations of each fold as burn-in time after visually checking for the appropriate mixing and convergence of the 10 chains (one chain per fold). This procedure produced a total of 400 000 draws of the posterior distribution of all of the parameters in the hierarchical model. We used relatively uninformative Beta priors for all the probability parameters in the hierarchical model ~Beta (1,1) and an uninformative Normal prior for the effects of the monitoring year, βoy and βhy, ~Normal (0, σ = 3.16). The script for fitting the Bayesian hierarchical model is available in Supplementary Appendix S1 (available online). We assessed the goodness of the fit and the predictive abilities of the hierarchical model by estimating a Bayesian p-value based on the deviance of the model (Hobbs & Hooten Reference Hobbs and Hooten2015, Kéry & Royle Reference Kéry and Royle2016, Conn et al. Reference Conn, Johnson, Williams, Melin and Hooten2018). A Bayesian p-value compares the estimate of a goodness-of-fit test – the deviance in this case – with the test expected value if the model is adequate. Extreme Bayesian p-values (i.e., <0.05 or >0.95) indicate a poor-fitting model (Hooten & Hobbs Reference Hooten and Hobbs2015, Kéry & Royle Reference Kéry and Royle2016, Conn et al. Reference Conn, Johnson, Williams, Melin and Hooten2018).