Study area

Bonaire covers 28 800 ha and lies between 12°18′N, 68°23′W and 12°01′N, 68°15′W (Fig. 2). Including the planned but not yet officially designated sections (1257 ha), the park covers 6900 ha in the northern part of the island, which is characterized by dry forest-shrub vegetation, and limestone soils with steep hills, such as Mount Brandaris, rising 243 m above sea level (De Freitas et al. 2005, Wells and Debrot 2008). Annual rainfall is highly variable (e.g. January–December 1968–2017: mean ± SD = 465.1 ± 201.2 mm; Fig. 3a). March is the driest (13.0 ± 15.2 mm) and November the wettest month of the year (95.2 ± 86.8 mm; Fig. 3b). In general, the dry season extends from January to September and the wet season from October to December (Fig. 3b). The onset of the wet season coincides with the productivity peak of plants that are important in the diet of the parrot (Dunn 2009, Simal and Bertuol unpubl.), and with fledglings leaving their nests at the end of August (Williams 2009, Williams and Evans 2010).

Figure 2.

Map showing the survey region (A = 17 000 ha) and the systematic sampling scheme (k=400 points) used for surveying the yellow-shouldered parrot on Bonaire in 2009–2017.

Figure 3.

(A) Total annual rainfall and (B) mean monthly rainfall ± SD at Cargill, southern Bonaire, 1968–2017. Consecutive wet years occurred in 2010 and 2011 (solid blue bars), with 2010 being the second wettest year on record. Consecutive dry years occurred in 2014–2016 (solid red bars), with 2015 being the driest year on record. Rainfall data were collected and provided by Cargill Bonaire, Inc. (D. Deanda, Jr. pers. comm.).

The vegetation of Bonaire is typical of dry islands in the eastern Caribbean, including common plants that provide food to the parrot, such as the golden-spined cactus Pilosocereus lanuginosus, apple cactus Subpilocereus repandus, candle cactus Stenocereus griseus, calabash Crescentia cujete, pigeon berry Bourreria succulenta, lignum-vitae Guaicum sanctum, olive wood Capparis odoratissima, gumbo limbo Bursera simaruba, logwood Haematoxylum brasiletto, mesquite Prosopis juliflora, divi-divi Caesalpinia coriaria, West Indian cherry Malpighia emarginata and indigo berry Randia aculeata. Mango Mangifera indica, nasberry Manilkara zapota, tamarind Tamarindus indica, neem Azadirachta indica and other plants with ornamental, medicinal and/or nutritional value are commonly found at residential gardens and farms, and also provide food to the parrot. Small-scale agricultural farming with corn Zea mays and other crops can be found across the island. Most farmers consider the parrot an agricultural pest (Parks 2010). Refer to De Freitas et al. (2005) and Wells and Debtor (2008) for additional information about Bonaire's landscape, vegetation and areas important to bird conservation.

Distance sampling

We used 100-ha grid cells to establish 400 points systematically across the range of distribution of the parrot on Bonaire (Fig. 2). We made abundance inferences (i.e. = × A) from the sampled to the targeted population (i.e. from the estimated mean density of parrots at surveyed points to the estimated mean number of parrots in survey region A = 17 000 ha). A team of two experienced observers conducted 6-min counts at surveyed points to increase the probability of detecting calling parrots visually for accurate distance measurement using binocular rangefinders. When parrots were heard but not seen, we measured detection distances to the nearest locations, and used the following distance categories: 0–15, 16–30, 31–45, 46–60, 61–90, 91–120, 121–180, 181–240, 241–340, 341–440 m (Rivera-Milán et al. 2005).– We did not include moving parrots in density estimates, unless their initial locations were known exactly (i.e. to the nearest meter) or approximately (i.e. within a distance category). To minimize responsive movement, we did not implement a settling period before count initiation. To maximize parrot availability for detection (Burnham et al. 2004, Rivera-Milán et al. 2005, Nichols et al. 2009), we conducted the surveys from sunrise to midmorning and from midafternoon to sunset, and visited the points twice during two-week sampling periods before-after reproduction. Survey effort accounted for the number of visits per point during each sampling period (Buckland et al. 2001).

We modelled detection probability as a function of distance r and other covariates represented by vector z (g[r, z]; Marques et al. 2007). Density was estimated using

where n is the number of detections of single parrots and clusters of parrots; is the average cluster size, which was used for density estimation when cluster detection was not size biased; and k is the number of surveyed points. For conventional distance sampling, when cluster detection was sizebiased (p < 0.15), we regressed log(s_i) on ĝ(r_i) to estimate the value of expected cluster size [Ê(s)] where ĝ(r_i) = 1, and Ê(s) instead of was used for density estimation (Buckland et al. 2001). After right-truncation of the distance data (ω = 240 m), detection probability of the ith single parrot or cluster of parrots was estimated as

We evaluated the fit of uniform, half-normal, and hazardrate key functions with and without cosine and polynomial series expansions using quantile–quantile plots and goodness-of-fit tests (Burnham et al. 2004). Based on the minimization of Akaike information criterion (AIC_c) and the desired coefficient of variation for the density estimator (CV < 0.20), we selected the half-normal key function without series expansion to explore the influence of the following detection covariates; year (2009–2017), sampling period (pre-reproduction versus post-reproduction), detection mode (audio versus visual), detection time (0–3 min versus 4–6 min), time of day (06:30–10:30 versus 15:30– 18:30), point location (on road versus off road), disturbance level (none-low versus medium-high), vegetation cover (open versus closed), land cover (forest-shrub versus agriculture-urban), predator and competitor presence (detected versus undetected), and cluster size (≥2 parrots within 10 m of each other, showing similar behaviour).

After 6-min counts, we moved around point centres as much as needed to reach a consensus about land covers, vegetation covers, disturbance levels, and food availability levels, which were also recorded as none-low or mediumhigh, depending on the number of plants bearing greenleafed shoots, flower buds, fruits and/or seeds known or assumed to be eaten by the parrot. We checked thoroughly for parrots undetected near point centres. When in doubt, we approached single parrots and parrot clusters to verify detection distances and cluster sizes. We used the binocular rangefinders to determine vegetation covers and food availability levels within 60 m and land covers and disturbance levels within 200 m of point centres. We recorded the presence of introduced mammals, fires, deforestation and other anthropogenic disturbances resulting from hunting, poaching, littering, farming, urbanization and/or recreational activities. We completed the assessment of land covers and disturbance levels around point centres with Google Earth Pro, ver. 7.1 (< www.google.com/earth/explore/products/desktop.html>).

We collected the data in a manner that would allow the use of different count methods for detection and abundance estimation (e.g. 3-min snapshot counts, two 3-min timeremoval counts, or two-week repeated counts; Burnham et al. 2004, Royle 2004, Buckland 2006, Nichols et al. 2009, Rivera-Milán and Simal 2012, Amundson et al. 2014, Buckland et al. 2015; the methods were recently reviewed by Dénes et al. 2018). However, because distance sampling is pooling robust, and our main interest was to estimate abundance across the entire survey region for population monitoring and modelling, we used multiple-covariate distance sampling to explore the influence of detection covariates, and conventional distance sampling to compare parrot mean density estimates at surveyed points after stratification (e.g. inside versus outside the park) or post-stratification (e.g. on road versus off road) of the distance data (Buckland et al. 2001, 2015, Marques et al. 2007). Results are presented as means ± SE with 2.5% and 97.5% quantiles of the bootstrap estimates. For data analysis, we used DISTANCE ver. 7.0.1 (Thomas et al. 2010).

Density hypotheses

We used the two-tailed Z test to compare parrot mean density estimates at surveyed points (Buckland et al. 2001). Because parrots are widely distributed across the survey region, we did not expect to find a significant difference (p > 0.05) between mean density estimates at surveyed points inside and outside the park. For the same reason, we did not expect to find significant differences between mean density estimates at surveyed points along and away from roads, regardless of land covers and disturbance levels. However, because food can be scarce and spatially clumped in a dry island environment, we expected to find a significant difference between mean density estimates at surveyed points with different food availability levels. In addition, because dense vegetation can provide protection against raptors, and because limited foraging and nesting resources can promote negative interspecific interactions, we expected to find significant differences between parrot mean density estimates at surveyed points in open and closed vegetation cover with and without potential avian predators and competitors, such as the peregrine falcon Falco peregrinus, merlin F. columbarius, crested caracara Caracara cheriway, pearly-eyed thrasher Margarops fuscatus bonairensis, troupial Icterus icterus, tropical mockingbird Mimus gilvus, and brown-throated parakeet Aratinga pertinax xanthogenia.

Lastly, because rainfall can influence food availability, which, in turn, can influence parrot survival and reproductive rates, we expected to find a positive relationship between post-reproduction mean density estimates during the wet season in October–December 2009–2017 and total rainfall during the dry season in January–September 2008–2016. Details of statistical analyses are provided with corresponding results. For data analysis, we used R ver. 3.3.2 (Crawley 2007, < www.r-project.org>).

Bayesian state-space logistic model

We modelled the dynamics of the parrot population, recognizing that distance sampling survey data were only available for 2009–2017, and that abundance estimates were subject to error. We used a Bayesian state–space modelling framework to account for observation and process error variances (Meyer and Millar 1999, Millar and Meyer 2000, Kéry and Schaub 2012, Hostetler and Chandler 2015, Rivera-Milán et al. 2016). The state dynamics of the population were modelled with a discrete form of the standard logistic equation in the presence of human-induced mortality from illegal hunting and other forms of human–parrot conflicts resulting in purposeful or accidental deaths in 2018–2066. We calculated annual changes in population state N_t with

where r_max is the maximum intrinsic rate of population growth, K is the population carrying capacity, N_t is the true unknown state of the parrot population, and M_t is the total number of deaths from human-induced mortality in year t. The total number of deaths from human-induced mortality M_t = N_tm_t, where m_t is the mortality rate between time period t and t+1 (Runge et al. 2009). Because parrot demographic data are insufficient or unavailable at this time, and because we cannot ascertain the relative importance or magnitude of illegal hunting and other human–parrot conflicts on population dynamics, human-induced mortality rates were randomly generated as part of the Markov chain Monte Carlo (MCMC) algorithm using uniform distributions for three potential scenarios; low mortality rates: m ∼ uniform (0.001, 0.100); moderate mortality rates: m ∼ uniform (0.101, 0.250); and high mortality rates: m ∼ uniform (0.251, 0.500). We reparameterized the unknown population state as a proportion of population carrying capacity NJK to reduce the autocorrelation in the MCMC samples used for Bayesian estimation (Meyer and Millar 1999, Millar and Meyer 2000). We assumed that the error of state model predictions ε was lognormally distributed with mean 0 and an estimated standard deviation σ_process. Based on this reparameterization, the state dynamics were projected forward in time according to

We modelled the population proportion in year 1 using a lognormal distribution with mean P₀ and variance . That is,

While the process model in the state–space formulation accounted for incomplete understanding of parrot population dynamics, we also needed to relate true population state to the surveys accounting for observation error. We specified that population size y_t and observation error were directly estimated from the distance sampling survey data (Knape 2008). Because the distribution of distance sampling abundance estimates tends to be positively skewed, we assumed a lognormal distribution for the observation error (Buckland et al. 2001). We transformed the population size estimates to the natural logarithm scale by transforming the bootstrap standard error to the standard deviation of the corresponding lognormal distribution. To complete the observation model of the state-space formulation, true unknown population state N_t = P_tK was related to observed population estimates with

where

Assuming linear density dependence (Runge et al. 2009), we derived the following parameters:

where m_ms is the maximum sustainable human-induced mortality rate, _eq is the equilibrium population size, and M_ms is the maximum sustainable total number of deaths from human-induced mortality. With this model formulation, we assumed that all human-induced mortality occurred after reproduction, and that all age classes (juveniles, subadults, adults) had equal mortality probability. In addition, we assumed additive mortality from natural and anthropogenic disturbances, although the model formulation allowed for compensatory response through density-dependent population growth (Runge et al. 2009).

We used the Bayesian state–space modelling framework to estimate the unobserved population parameters and unknown population states, assuming conditional independence for each time step. We specified uniform prior distributions with wide but biologically realistic bounds for maximum population growth rate (r_max ∼ uniform [0.010, 0.500]), population carrying capacity (K ∼ uniform [1000, 10 000]), and the mean of the initial population proportion on the log scale (P₀ ∼ uniform [–2, 0]). For the process error and initial population proportion standard deviations, we also specified uniform priors (σ_process and σ_P0 ∼ uniform [0, 2]).

We assessed demographic sustainability (i.e. births ≥ deaths in a geographically-isolated population) in terms of the probability that predicted abundance would exceed 2000 parrots in 2066, given low, moderate and high humaninduced mortality rates (i.e.P[N > 2000|m = 0.001–0.500]; Reed and Hobbs 2004, Traill et al. 2010). We used the 2.5% percentile of carrying capacity to establish a populationbased conservation objective for baseline and evaluation monitoring (Sanderson 2006, Rivera-Milán et al. 2016). To estimate the posterior distributions of population parameters, we used MCMC methods by running JAGS, ver. 3.4.0 (Plummer 2003) within R2JAGS (Su and Yajima 2015). We conducted 250 000 iterations and used the first 50 000 iterations as a burn-in period. We generated three Markov chains with different initial parameter values and used trace plots and node summary statistics to check for MCMC algorithm convergence (e.g. Brooks–Gelman–Rubin diagnostic statistic = 1.00; Gelman et al. 2004). Markov chains were thinned by 25 to obtain samples of 8000 points. Results are presented as means and MCMC SD and medians with 2.5% and 97.5 percentiles.

Distance sampling

We surveyed 303 points (k) and made 542 detections (n) of single parrots and clusters of parrots within 240 m of point centres before reproduction in March 2010 and 2012 and after reproduction in October–December 2009–2017. The half-normal key function without series expansion provided the best fit to the distance data (Kolmogorov–Smirnov test: D_n = 0.17, p = 0.47; Cramer–von Mises family tests: W² = 0.15, p > 0.40, C² = 0.09, p > 0.40; AIC = 3949.26). Mean density estimates (i.e. parrots ha^-1) were similar for the top ranking conventional distance sampling detection models. That is, estimated mean density was 0.171 ± 0.020 bootstrap SE (2.5% and 97.5% quantiles = 0.136–0.215) for the half-normal key function without series expansion, 0.174 ± 0.022 (0.136–0.222; ΔAIC_c = 1.24) for the uniform key function with two cosine series expansions, and 0.178 ± 0.045 (0.107–0.295; ΔAIC_c = 2.38) for the hazardrate key function without series expansion.

Based on the half-normal key function without series expansion, mean cluster size was 2.417 ± 0.102 (2.225– 2.625) and expected cluster size was 2.077 ± 0.067 (1.950–2.213). Cluster detection was size biased (Pearson's correlation coefficient r = –0.12, p = 0.01). Mean encounter rate (n / k) was 0.656 ± 0.055 (0.557–0.773). Mean detection probability was 0.442 ± 0.014 (0.415–0.470) within 240 m and mean effective detection radius (w×) was 159.559 ±4.308 (151.336–168.229 m).

The half-normal key function without series expansion and vegetation cover and land cover defined as two-level factor covariates were the top ranking detection models (ΔAIC_c < 2.00; Table 1). Mean detection probability tended to be higher at surveyed points in land covers with open vegetation cover (agriculture–urban: = 0.492 ± 0.026, 0.441–0.543; forest–shrub: = 0.445 ± 0.027, 0.392–0.498) than closed vegetation cover (agriculture– urban: = 0.434 ± 0.027, 0.381–0.487; forest–shrub: = 0.391 ± 0.028, 0.336–0.445; Fig. 4). The halfnormal detection models with other covariates did not receive support from the distance data (e.g. point location; ΔAIC_c ≥ 18.28), but generated similar mean density estimates (Table 1). Based on model averaging of the top ranking detection models (Table 1), mean density estimate was 0.172 ± 0.020 (0.137–0.216), and mean population size estimate was 2924 ± 340 (2329–3672) parrots in the survey region.

Table 1.

Half-normal detection models with Akaike information criterion values (truncation distance w = 240 m) and mean density estimates and standard errors with 2.5% and 97.5% quantiles of bootstrap estimates based on the yellow-shouldered parrot distance-sampling survey data before-after reproduction on Bonaire in 2009–2017. Sample size n = 542 detections.

Figure 4.

Half-normal key function without series expansion and with vegetation cover and land cover defined as two-level factor covariates for the yellow-shouldered parrot distance-sampling survey data on Bonaire in 2009–2017. The blue line is for open vegetation cover and the green line is for closed vegetation cover in urban–agriculture land cover. The yellow line is for open vegetation cover and the red line is for closed vegetation cover in forest–shrub land cover. Sample size n = 542 detections.

Parrot detection was not influenced by sampling period (ΔAIC_c = 60.25) nor year (ΔAIC_c = 63.74; Table 1). Mean density and population size estimates were 0.165 ± 0.026 (0.121–0.224) and 2805 ± 438 (2057–3808) in March 2010, and 0.192 ± 0.029 (0.136–0.250) and 3264 ± 492 (2306–4242) in October 2010. That is, mean population growth rate before–after reproduction was 0.164 ± 0.036 (0.107–0.251) in 2010. Mean density and population size estimates were 0.161 ± 0.025 (0.120–0.215) and 2737 ± 405 (2040–3655) in March 2012, and 0.182 ± 0.027 (0.128–0.235) and 3094 ± 464 (2184–4001) in October 2012. That is, mean population growth rate before-after reproduction was 0.130 ± 0.028 (0.086–0.197) in 2012.

Density hypotheses

Parrot distribution was clumped (i.e. dispersion parameter b > 1), particularly at surveyed points with medium-high food availability levels (b = n×[CV D]² = 321 × 0.158² = 8.021; Table 2). As expected, parrot mean density estimates were significantly higher at surveyed points with medium-high than none-low food availability levels (Table 2). Mean density estimates also were significantly higher at surveyed points with closed than open vegetation cover (Table 2). Contrary to expected, mean density estimates were significantly higher at surveyed points inside than outside the park (Table 2). However, mean density estimates did not differ significantly at surveyed points along and away from roads, regardless of land covers and disturbance levels, or the presence of potential avian predators and competitors (Table 2).

Also as expected, there was a positive relationship between post-reproduction mean density estimates in October–December 2009–2017 and total rainfall in January– September 2008–2016 (simple linear regression: r² = 0.57, β₁ = 0.146 ± 0.048, 0.078–0.274, t=3.05, p = 0.02). During wet years in 2010 and 2011, parrot mean density estimates were 0.192 ± 0.029 (0.136–0.250) and 0.193 ± 0.029 (0.136–0.251), and mean population size estimates were 3264 ± 492 (2306–4242) and 3281 ± 492 (2306– 4262; Fig. 3a). However, during dry years in 2014–2016 (Fig. 3a), parrot mean density and population size estimates decreased to 0.158 ± 0.039 (0.080–0.234) and 2679 ± 670 (1358–3979). That is, the parrot population declined significantly between consecutive wet and dry years (simple linear regression: r² = 0.70, β₁ = –0.023 ± 0.006, 0.011– 0.035, t = –4.16, p = 0.006; Fig. 3a, 5a).

Bayesian state-space logistic model

Realizations of Markov chains and node summary statistics showed convergence of the MCMC algorithm (Brooks– Gelman–Rubin diagnostic statistic = 1.00–1.01). The mean for maximum population growth rate was 0.179±0.129 MCMC SD (median r_max= 0.146, 2.5% and 97.5% percentiles = 0.017–0.468). Mean population carrying capacity was 5623 ± 2043 parrots (median K=5113, 2853–9667). That is, the population-based conservation objective (or desired target level) was above 2800 parrots in the survey region. Mean maximum sustainable human-induced mortality rate was 0.090 ± 0.064 (median m_m = 0.073, 0.008–0.234). Mean maximum sustainable total number of deaths from human-induced mortality was 219 ± 135 parrots (median M_ms= 200, 26–503). Mean equilibrium population size was 2811 ± 1022 parrots (median N_eq = 2557, 1427–4833). Mean process error variance was 0.012 ± 0.025 (median = 0.004, 0.00001–0.078).

Table 2.

Mean density estimates and standard errors with 2.5% and 97.5 quantiles of the bootstrap estimates and Z scores with p-values after stratification or post-stratification of the yellow-shouldered parrot distance-sampling survey data on Bonaire in 2009–2017.

Figure 5.

(A) Yellow-shouldered parrot distance-sampling population size mean estimates ± bootstrapped SE in October–December 2009–2017 (blue crosses with vertical lines), and Bayesian state– space logistic model means ± MCMC SD (red circles with vertical lines) and medians with 2.5% and 97.5 percentiles (red lines) and low human-induced mortality rates (m ∼ uniform [0.001, 0.100]) on Bonaire in 2018–2066; and (B) frequency distribution of model-based abundance predictions in 2066.

With low to moderate human-induced mortality rates (0.001– 0.100, 0.101–0.250), predicted population sizes were 2963 ± 668 parrots (median N = 2879, 1925–4474) and 2754 ± 690 parrots (median N = 2717, 1557–4141) in 2018, and 2703 ± 1660 parrots (median N = 2688, 243–6576) and 2297 ± 1301 parrots (median N = 2352, 29–4874) in 2066 (Figs. 5a, 6a). With low to moderate human-induced mortality rates, the probabilities of population size exceeding 2000 parrots in the survey region were 0.717 and 0.666 (Fig. 5b, 6b). With high human-induced mortality rates (0.251–0.500), predicted population sizes were 1780 ± 1160 parrots (median N = 1547, 431–4443) in 2018 and 26 ± 139 parrots (median N = 2, 0–191) in 2066 (Fig. 7a). With high human-induced mortality rates, the parrot population probably would not exceed 500 parrots in the survey region (Fig. 7b).

Distance sampling

We used conventional and multiple-covariate distance sampling to estimate parrot detection probability and abundance in 2009–2017, and population growth rate before-after reproduction in 2010 and 2012. Vegetation cover and land cover were the most important detection covariates, with parrots being harder to detect in closed forest–shrub than in open agriculture–urban vegetation cover. However, detection probability was similar in open forest–shrub and closed agriculture–urban vegetation cover. None of the detection covariates, including point location (on road versus off road) and detection mode (audio versus visual), caused extreme heterogeneity in the detection function (Marques et al. 2007; Fig. 3). In addition, conventional and multiplecovariate detection models generated similar mean density estimates, suggesting that model selection was of secondary importance for abundance inferences in the survey region (Buckland et al. 2001).

Figure 6.

(A) Yellow-shouldered parrot distance-sampling population size mean estimates ± bootstrapped SE in October–December 2009–2017 (blue crosses with vertical lines), and Bayesian statespace logistic model means ± MCMC SD (red circles with vertical lines) and medians with 2.5% and 97.5 percentiles (red lines) and moderate human-induced mortality rates (m ∼ uniform [0.101, 0.250]) on Bonaire in 2018–2066; and (B) frequency distribution of model-based abundance predictions in 2066.

Figure 7.

(A) Yellow-shouldered parrot distance-sampling population size mean estimates ± bootstrapped SE in October–December 2009–2017 (blue crosses with vertical lines), and Bayesian statespace logistic model means ± MCMC SD (red circles with vertical lines) and medians with 2.5% and 97.5 percentiles (red lines) and high human-induced mortality rates (m ∼ uniform [0.251, 0.500]) on Bonaire in 2018–2066; and (B) frequency distribution of model-based abundance predictions in 2066.

The purpose of having two observers, with one observer mainly measuring detection distances and the other mainly recording the data, was to meet the basic model-based assumptions of conventional and multiple-covariate distance sampling. That is, certain detection at point centres (g[0] = 1); accurate measurement of detection distances to initial locations; and correct estimation of cluster sizes (Buckland et al. 2001, 2008, 2015, Burnham et al. 2004, Marques et al. 2007). Detection probability remained high near point centres (e.g. p > 0.729 within 60 m in closed forest–shrub cover). We did not likely miss parrots at point centres during 6-min counts; and detection distances to initial locations were measured with little error, as evidenced by quantile–quantile plots and goodness-of-fit tests (Burnham et al. 2004; Fig. 11.9). We did not include moving parrots in density estimates, unless detection distances or distance categories were determined prior to responsive movement (Rivera-Milán et al. 2005), which was usually away from the observers with loud and easy to locate squawks. Cluster detection was size biased (i.e. large clusters were more detectable than small ones at longer distances from point centres); but parrot clusters tended to be small and had a negligible effect on density estimation. In addition, we surveyed on repeated occasions a systematic grid of points across the parrot range of distribution to increase encounter rate and justify abundance inferences in the survey region. Similar mean density estimates were obtained with snapshot, time-removal, and repeated counts (Rivera-Milán and Simal 2012); and independently by Saunders (2011) using conventional distance sampling. Therefore, we suggest that the design-based and model-based assumptions of conventional and multiple-covariate distance sampling were met, and that unbiased abundance estimates were obtained for monitoring and modelling the dynamics of Bonaire's parrot population.

We estimated population growth rate before–after reproduction in 2010 and 2012. Our abundance estimates suggest that the parrot population remained relatively stable, with births offsetting deaths during average and above-average rainfall years in 2010–2012, resulting in a mean population size estimate of about 2800 parrots pre-reproduction and 3200 parrots post-reproduction. That is, mean population growth rate was about 0.15 before-after reproduction under favourable conditions in 2010–2012. This population growth rate seems biologically plausible, despite the fact that parrot populations typically contain large proportions of nonbreeding adults (Renton and Salinas-Melgoza 2004, Wiley et al. 2004, Stahala 2005, Beissinger et al. 2008); and that the number of breeding pairs may be limited by the availability of suitable nesting sites on small islands (Williams 2009, Richards 2010, Roberts et al. 2014). For example, on the basis of life-history characteristics typical of the family Psittacidae (i.e. delayed maturity and long lifespan; Forshaw 2010) and an average reproductive rate of 1.89 young per breeding pair in 2006–2010 (Williams and Evans 2010), and using a modified version of Moffat's equilibrium model (Moffat 1903, Hunt 1998, 2003), population growth rate ranged from 0.10 to 0.28, assuming survival rates were 0.60–0.85 for juveniles and 0.75–0.90 for sub-adults (i.e. hatch-year and 1–3 years old), and 0.80–0.95 for breeding and nonbreeding adults (4–20 years old). However, these demographic rates can be lower or higher and vary temporally by age or sex, depending on rainfall variability and food availability, intraspecific and interspecific negative interactions, and human-parrot conflicts (Dunn 2009, Martin 2009, Williams 2009, Parks 2010, Richards 2010, Williams and Evans 2010, Roberts et al. 2014). In addition to longterm count data, telemetry and nest-monitoring data are needed to estimate demographic parameters and determine their relative importance with respect to realized population growth rates before–after dry and wet years (e.g. breeding and nonbreeding female survival rates; Stahala 2005). Our survey-based and model-based estimates of population growth rate suggest that the parrot can have high survival and reproductive rates during average and above-average rainfall years.

Density hypotheses

We tested specific hypotheses about parrot mean density estimates at surveyed points, and the relationship between post-reproduction mean density estimates and total rainfall during the first nine months of the previous year. The parrot behaved as a habitat generalist, feeding on green-leafed shoots, flower buds, fruits and/or seeds of a wide variety of plants, including columnar cacti and leguminous and nonleguminous trees and shrubs across the survey region. We suggest that habitat use, behavioural plasticity, and daily movements mainly reflected adaptive foraging strategies when parrots searched for scarce and spatially clumped food (e.g. social interactions among neighbours to inform foraging decisions and adapt according to temporal, energetic and cognitive constraints; Stephens and Krebs 1986). These adaptive foraging strategies may explain why mean density estimates did not differ at on-road and off-road surveyed points, regardless of land covers and disturbance levels, or the presence of potential avian predators and competitors.

Mean density estimates differed at surveyed points inside and outside the park, as well as at surveyed points with open and closed vegetation covers and different food availability levels. In general, estimated mean density was highest at surveyed points with medium-high food availability and closed forest–shrub vegetation cover inside the park. Parrot distribution was clumped, particularly at surveyed points with medium-high food availability levels; and there was a significant positive relationship between post-reproduction mean density estimates during the wet season in October— December 2009–2017 and total rainfall during the previous year dry season in January–September 2008–2016. Based on these results, we suggest that rainfall during the dry season plays a primary function as a proximate factor that can influence vegetation cover and food availability, which, in turn, can influence parrot survival and reproductive rates, and ultimately determine before–after reproduction population growth rate and abundance during the wet season. We also suggest that mean density estimates were higher inside than outside the park because of the greater diversity and abundance of plants that provide cover and food to the parrot during the breeding and nonbreeding periods. However, additional data are needed to quantify and compare plant species diversity, abundance and seasonal productivity inside and outside the park (Simal and Bertuol unpubl.); and to determine whether rainfall and food can mediate intraspecific and interspecific negative interactions that ultimately limit and regulate parrot abundance (Dunn 2009, Martin 2009, Williams 2009, Williams and Evans 2010).

Bayesian state-space logistic model

The model was useful for estimating the posterior distributions of population parameters and predicting changes in abundance in 2018–2066. Although we made a number of simplifying assumptions (e.g. additive and equal mortality among age classes after reproduction, linear density dependence), our primary focus was on modelling changes in parrot abundance across the entire survey region, and therefore restricting the model to essential parameters sufficed for estimation and prediction, while explicitly accounting for observation and process error variances. Model-based abundance predictions and distance sampling abundance estimates were relatively precise (mean CV = 0.140) and quite close in 2009–2017 (mean ΔN=111 parrots), showing a significant population decline between consecutive wet and dry years. However, the posterior distributions of maximum population growth rate (mean CV = 0.718), population carrying capacity (mean CV = 0.363), and abundance predictions in 2066 (mean CV = 2.185) were relatively imprecise mainly as a result of process error variance (mean CV = 1.724). Precision can be gained by modelling longer time series using distance sampling abundance estimates (Rivera-Milán et al. 2016); and by combining count, telemetry and nest-monitoring data using integrated population models (Johnson et al. 2010). In addition, combining multiple data types can help us identify the most important demographic parameters (e.g. survival versus reproductive rates) and drivers (natural versus anthropogenic disturbances) of parrot population dynamics in a changing environment (Oppel et al. 2014).

With low to moderate human-induced mortality rates the parrot population would likely recover from the decline in 2009–2017; and probably would exceed 2000 parrots in the survey region in 2018-2066 (i.e. p[N > 2000|m = 0.001–0.250] = 0.666–0.717). In contrast, with high human-induced mortality rates (m > 0.250), the population probably would decline and become too small to be self-sustainable. Because mortality rates from human–parrot conflicts may be unsustainable (i.e. m > m_ms), we consider Bonaire's parrot population vulnerable to the risk of extinction. However, because data about human-induced mortality are lacking, our abundance predictions should be taken only as suggestive of demographic sustainability during the modelled time horizon. Therefore, given the current state of knowledge and concerns, we consider quantification and control of human–parrot conflicts a conservation priority (e.g. poaching of nestlings inside and outside the park, Roberts et al. 2014; and hunting at agricultural farms, Parks 2010).

We used a simple modelling framework and focused on abundance estimation and prediction in the survey region rather than on elucidation of the relative importance of the factors and mechanisms driving population dynamics (e.g. scramble versus contest competition for limited resources; Henson and Cushing 1996). Ideally, we would like to use multiple data types to develop integrated population models associated with specific hypotheses about changes in abundance and demographic rates during the breeding and nonbreeding periods. However, demographic data are insufficient or unavailable at this time, and additional model complexity would require greater speculation, regardless of the framework used (Johnson et al. 2012, Kéry and Schaub 2012, Sekeris 2012). Long-term monitoring using multiple data types is needed to better understand population dynamics. For example, an increase in the number of consecutive dry years in the Caribbean can cause significant local changes in vegetation structure and composition, reducing food availability, promoting negative intraspecific and interspecific interactions, and exacerbating human–parrot conflicts. This, in turn, can increase the amplitude of population fluctuations at lowered parrot numbers, with the strength of density dependence determining realized population growth rates and the extent to which the population will remain below carrying capacity (Reed and Hobbs 2004, Abadi et al. 2012, Sæther et al. 2016).

Because model-based abundance predictions can be updated regularly with distance sampling survey data, we can continue learning from the comparison of predicted and estimated abundances before–after dry and wet years, as well as before-after conservation actions taken inside and outside the park. Therefore, we recommend long-term monitoring and modelling for assessing changes in abundance and the results of conservation actions taken to keep the population above 2800 parrots in the survey region (i.e. N > 2.5% percentile K). We also recommend using multiple data types and an adaptive management approach (e.g. experimentation based on before–after-control–impact design, Stien and Ims 2015) to accelerate learning about demographic rates and factors influencing changes in abundance from conservation actions (e.g. removing goats and fencing to increase vegetation cover and food availability; and creating or restoring nesting cavities to increase the number and productivity of breeding pairs; Dunn 2009, Williams 2009, Williams and Evans 2010, Roberts et al. 2014). Lastly, our monitoring and modelling approach is flexible and has been used to assess other parrot and landbird populations of conservation concern in the Caribbean (Rivera-Milán et al. 2005, 2015, 2016, Haakonsson et al. 2017).