Puerto Rico plain pigeon, scaly-naped pigeon and red-tailed hawk: population dynamics and association patterns before and after hurricanes

: Since the 1980s, 3 major hurricanes have made landfall on Puerto Rico: Hugo in September 1989 (Saffir-Simpson scale, category 4), Georges in September 1998 (category 3) and María in September 2017 (category 4). María was the most devastating hurricane since the 3 major hurricanes that occurred in 1899−1932. Major hurricanes can cause severe abundance declines and population bottlenecks by de creasing survival and reproductive rates and increasing predation and competition for limited resources. In April to June 1986−2021, we used distance sampling to estimate abundance and monitor the population dynamics of the endangered Puerto Rico plain pigeon Patagioenas inornata wetmorei and the abundant scaly-naped pigeon P . squamosa and red-tailed hawk Buteo jamaicensis . Here, we fit a Bayesian state-space logistic model with distance sampling abundance estimates to generate posterior estimates of maximum population growth rate and population carrying capacity, and predict abundance in April to June 2020−2030. In addition, we used N -mixture and 2-species models to assess association patterns in April to June 2015−2019. The scaly-naped pigeon and red-tailed hawk populations did not decline, or recovered faster from their declines than the plain pigeon population after the hurricanes. The association patterns between species were positive but variable for the 2 pigeon species and negative but variable for the plain pigeon and red-tailed hawk. At lowered abundance (i.e. mean ± SE estimates N ˆ = 1043 ± 476 island-wide and N ˆ = 522 ± 157 at the centre of abundance in the east-central region in April to June 2018−2021), the plain pigeon may become extinct if another hurricane with the path and intensity of María makes landfall on the island during the current decade.


INTRODUCTION
Hurricanes can cause abundance declines, give rise to negative interspecific interactions and in crease the risk of extinction of Caribbean island birds with restricted distribution ranges and small population sizes (e.g. see Beissinger et al. 2008). Because of climate change, the frequency of major hurricanes (i.e. Saffir-Simpson scale, categories 3−5) is increasing in the warming waters of the Caribbean region (Goldenberg et al. 2001, Webster et al. 2005, Biasutti et al. 2012). The 2017 Atlantic hurricane season had 6 major hurri-canes, including Irma and María, which at maximum intensities reached category 5, with winds exceeding 252 km h −1 (US National Oceanic & Atmospheric Administration 2020). On 6 September 2017, Hurricane Irma passed about 100 km off the northeast coast of Puerto Rico, bringing as much as 300 mm of rainfall and category 2 winds in the range of 155−177 km h −1 affecting mostly the northeastern and east-central regions (US National Oceanic and Atmospheric Administration 2020). Two weeks later, Hurricane María made landfall as a category 4, with as much as 1500 mm of rainfall and winds in the range of 209− 251 km h −1 affecting the whole island (US National Oceanic and Atmospheric Administration 2020).
The most evident direct effect of hurricanes on Caribbean island forests is structural vegetation damage from heavy rains and high-speed winds (Boose et al. 2004, Uriarte et al. 2019, Hall et al. 2020). However, the widespread reduction of vegetation cover from up rooted trees, stem breaks and defoliation can limit food supplies during extended periods of time (e.g. see Wunderle 1995, their Fig. 3), particularly for fruitand seed-eating birds, such as the endangered Puerto Rico plain pigeon Patagioenas inornata wetmo rei. Reduced vegetation cover and limited food supplies, in turn, can cause abundance declines and population bottlenecks by de creasing decreasing survival and reproductive rates and increasing predation and competition for limited resources (e.g. see Beissinger et al. 2008, their Figs. 5 & 8).
Habitat loss and degradation from agriculture, predation and competition in remaining forest fragments, hunting and the effects of hurricanes San Ciriaco (category 4) in 1899, San Felipe (category 5) in 1928, San Ciprián (category 4) in 1932 and Betsy (category 3) in 1956, may have caused the plain pigeon population bottleneck during the 1920s− 1970s (Wetmore 1916, Danforth 1931, Koenig 1953, Pérez-Rivera 1978, Wiley 1985, Miyamoto et al. 1994, Boose et al. 2004). Because of the discovery of a small number of individuals in the east-central region of the island, the US Fish and Wildlife Service listed the plain pigeon as endangered in 1970 (USFWS 1982). Captive breeding was conducted during 1983−2001, with the collaboration of state and federal government agencies and the University of Puerto Rico at Humacao Campus, but reintroduction efforts were mostly unsuccessful for a variety of reasons, including predation by the red-tailed hawk Buteo jamaicensis (Ruiz-Lebrón et al. 1995, USFWS 2011. Concerned about the lack of long-term monitoring and reliable abundance estimates, we initiated island-wide columbid and raptor surveys once per year beginning in April to June 1986(Rivera-Milán 1992, 1995a,b, Rivera-Milán et al. 2003a. Since then, we have surveyed columbid and raptor populations annually in April to June, except in 2020 due to the coronavirus pandemic. During that time, 3 major hurricanes have made landfall on Puerto Rico: Hugo (category 4) in September 1989, Georges (category 3) in September 1998 and María (category 4) in September 2017. The plain pigeon population recovered after hurricanes Hugo and Georges (see Rivera-Milán et al. 2016, their Fig. 3). However, because María has been the most devastating hurricane to make landfall on the island since the 3 major hurricanes that occurred in 1899−1932 (Boose et al. 2004, Uriarte et al. 2019, Hall et al. 2020, we ex pected plain pigeon abundance to decline below the lowest estimate since monitoring started, which had occurred in April to June 1990(Rivera-Milán et al. 2016). Here, we used distance sampling to (1) estimate plain pigeon, scaly-naped pigeon P. squamosa and red-tailed hawk abundance island-wide and in the east-central region in April to June 1986−2021, (2) fit a Bayesian state-space logistic model with abundance estimates, (3) generate posterior estimates of maximum population growth rate and population carrying capacity and (4) predict abundance in April to June 2020−2030. We used the Bayesian state-space logistic model to project abundance forward in time to cover a decade, based on the precision of short-term abundance predictions and the time period needed for two 5 yr population status reviews (USFWS 2011, Rivera-Milán et al. 2016). In addition, following recommendations provided in the most recent population status review (USFWS 2011), we repeated the surveys 3 times per year at the centre of abundance of the plain pigeon in the east-central region in April to June 2015−2019. We used Nmixture and 2-species models to estimate detection, abundance and occupancy, and to assess association patterns in the east-central region before and after hurricanes Irma and María in September 2017.
While the endangered plain pigeon and the abundant scaly-naped pigeon are similar-sized resident congeners that may compete for limited resources (Pérez-Rivera 1978, Rivera-Milán 1996, 2001, Rivera-Milán et al. 2003b, the red-tailed hawk is an abundant resident raptor that preys on both pigeon species (Santana C. & Temple 1988, Rivera-Milán 1995b, Ruiz-Lebrón et al. 1995, Vilella & Nimitz 2012, Vilella et al. 2013, Gallardo et al. 2019. Therefore, mainly because of reduced vegetation cover and limited food supplies after the 2017 hurricanes, we expected plain pigeon detection, abun-dance and occupancy to be negatively associated with the presence and density of the scaly-naped pigeon and red-tailed hawk, resulting in 'detection δ and species ϕ interaction factors' (sensu MacKenzie et al. 2004MacKenzie et al. , 2006 below what would have been the case under independent coexistence (i.e. pre-hurricane δ and ϕ = 1, post-hurricane ϕ and δ < 1). However, although predator−prey and competitive interactions may influence community structure and population dynamics (Strong et al. 1984, Begon et al. 1996, Newton 1998, we recognize that alternative explanations about association patterns cannot be excluded with survey-based count data (i.e. pattern does not necessarily relate to process, and correlation or association does not necessarily imply relationship or causation; James & Boecklen 1984, Poling & Haysl ette 2006, Cooper et al. 2007). Here, we used the term 'interaction' in a statistical sense, meaning that δ and ϕ = 1 would suggest that the 3 species tended to occur independently of each other, based on estimates of detection, abundance and occupancy for the surveyed areas in the east-central region of Puerto Rico in April to June 2015−2019.

Distance sampling
We surveyed 1400 points island-wide once per year in April to June 1986−2021, including 247 points that were surveyed 3 times per year in the east-central region during 2015−2019 (see Fig. 1 in Rivera-Milán et al. 2016). The surveys covered 7490 km 2 island-wide and 1100 km 2 in the east-central region. We did not survey heavily urbanized and industrialized areas or primary roads (e.g. expressways). We located on-road points systematically with random starts along secondary and tertiary roads (e.g. paved and unpaved roads in rural areas). Off-road points were located systematically with random starts at a minimum distance of 200 m from the nearest roads. Minimum distance was 400 m between off-road points and 800 m between on-road points. We used 2-observer teams to meet the following distance sampling assumptions: (1) certain detection of pigeons and hawks at the point centres (i.e. g [0] = 1), (2) de tection at their initial locations, (3) correct distance measurements and (4) unbiased cluster size estimation (Buckland et al. 2001).
Two observers remained side by side at the point centres for 6 min, with 1 observer mainly recording the data and the other observer mainly measuring the distances to single birds or the geometric centre of clusters. We defined a cluster as 2 or more individuals of the same species within 10 m of each other, showing similar behaviour (e.g. perching). We used 6 min counts to increase the probability of pigeons and hawks becoming available for detection (e.g. calling). Whenever possible, we measured distances to the nearest metre with rangefinders. However, when we were not able to measure the distances to the nearest metre, distances were measured to the nearest locations and detections were allocated to categories 0− 15, 16−30, 31−45, 46−60, 61−90, 91−120, 121−180, 181−240, 241−340 or 341−440 m (Rivera-Milán et al. 2003a. We used distances to the nearest metre and category midpoints for the analysis of grouped and ungrouped distance data (Buckland et al. 2001). The second and third distance sampling assumptions (see above) were relaxed with 2 observers operating as a team and the use of distance categories (Buckland et al. 2001, Burnham et al. 2004. Regarding the first and fourth assumptions (see above), it was unlikely that 2 experienced observers missed pigeons or hawks at the point centres during 6 min counts, and detections were mostly of single birds or small clusters during the surveys (Rivera-Milán et al. 2003a. Flying pigeons and hawks were not included in density estimates, unless we measured the distances to their initial locations, before any movement. To minimize responsive movement, which usually was away from the point centres, we did not implement a settling period before count initiation. We conducted the surveys during 06:30−10:00 and 15:30−19:00 h. Survey effort accounted for the number of visits per point each year (Buckland et al. 2001).
We modelled detection probability as a function of distance r and covariates represented by vector z (i.e. g[r, z]; Marques et al. 2007). Density was estimated as: ( 1) where D is the number of individuals km −2 , n is the number of single and cluster detections, s ⁻ is the average cluster size, which was used for density estimation when cluster detection was not size biased, and k is the number of points. When cluster detection was size biased, we used the size-bias regression method (Buckland et al. 2001), or included cluster size as a covariate in the models. After exploratory data analysis with distances grouped and ungrouped for each species separated and years separated or combined (Buckland et al. 2001, Marques et al. 2007), we right-truncated the distance data at w = 120 m. Detection probability was estimated as: We estimated population size by extrapolating density estimates from the surveyed area (i.e. kπw 2 ) to the survey region (i.e. N = D × A, where A = 7490 km 2 island-wide or 1100 km 2 for the east-central region). We used nonparametric bootstrapping to estimate detection and abundance standard errors (SEs) (Buckland et al. 2001).
After the 6 min counts, the 2 observers moved around the point centres as much as needed to determine if nearby pigeons or hawks were missed, verify any problematic distance measurements and reach a consensus about vegetation cover and food availability within 60 m. For more information about the diet of the 2 pigeon species, see Pérez-Rivera (1978), Rivera-Milán (1996, 2001 and Rivera-Milán et al. (2003b). We post-stratified the distance data of each species by year (e.g. prehurricane 0 = 2015−2017, post-hurricane 1 = 2018− 2021), vegetation cover and food availability (0 = none−low, 1 = medium−high), and used the 2tailed z-statistic to test for significant differences (p < 0.05) in detection and density estimates (Buckland et al. 2001). We used the programs DIS-TANCE version 7.3.1 (Thomas et al. 2010) and R version 3.6.3 (R Development Core Team 2020), and the R package 'DISTANCE' version 1.0.2 (Miller et al. 2019). Results are presented as means with bootstrapped SEs.

N-mixture and two-species models
Based on AIC c (Burnham & Anderson 1998), we selected the Poisson-binomial mixture model (Royle & Nichols 2003, Royle 2004, MacKenzie et al. 2006) to estimate detection probability P ij and abundance λ at 247 points R that were surveyed 3 times per year in the east-central region during 2015−2019. The points were surveyed once in 2021. We defined the repeated counts C ij as conditionally independent binomial random variables with detection and abundance at point i and time j (i.e. P ij , N ij ), and integrated likelihood: (3) We used an additive normal random effect to account for extra-Poisson variation in abundance. For example, log(λ ij ) = β 0 + β 1 (x ij ) + β 2 (y ij ) + ε ij for covariates x ij and y ij , with error ε ij ~ Normal (0, σ 2 λij ). We used a logistic regression model to account for extra-binomial variation in detection probability. For example, logit(P ij ) = β 0 + β 1 (x ij ) + β 2 (y ij ) + ε ij , for covariates x ij and y ij , with error ε ij ~ Normal (0, σ 2 Pij ). We expected plain pigeon abundance and detection estimates to be negatively associated with the pointlevel density estimates of the scaly-naped pigeon and red-tailed hawk (i.e. , where w = 120 m). We standardised the density estimates (i.e. In addition, we treated the point-level count data as non-detection or detection (i.e. 0 or 1), and used 2species models to estimate species ϕ and detection δ interaction factors: (4) and: ( 5) where ψ A is the probability of occupancy by the plain pigeon (i.e. species A), regardless of occupancy by the scaly-naped pigeon or red-tailed hawk (i.e. species B); ψ B is the probability of occupancy by B, re - gardless of occupancy by A; and ψ AB is the probability of AB co-occupancy (MacKenzie et al. 2004(MacKenzie et al. , 2006. P A j is the probability of detection for A during the j th survey, given AB presence; P B j is the probability of detection for B during the j th survey, given AB presence; and P j AB is the probability of AB co-detection during the j th survey, given AB presence (MacKenzie et al. 2004(MacKenzie et al. , 2006. While ϕ and δ = 1 would suggest independent cooccupancy and co-detection (i.e. ψ AB = ψ A × ψ B and P j AB = P A j × P B j ), ϕ and δ < 1 would suggest negative and ϕ and δ > 1 would suggest positive co-occupancy and co-detection. We also estimated the probability of occupancy and detection of A during the j th survey, given A presence but B absence (denoted 'b') (i.e. ψ Ab and P j Ab ). With ϕ and δ < 1, we would expect ψ Ab and P j Ab > ψ AB and P j AB ; that is, higher plain pigeon occupancy and detection with scaly-naped pigeon and red-tailed hawk absence than presence at the surveyed points. For additional information about the formulation of 2-species models, see MacKenzie et al. (2004MacKenzie et al. ( , 2006. We used the single-season Poissonbinomial mixture and 2-species models as implemented in the program PRESENCE version 2.13.6 (Hines 2020) and the R package 'UNMARKED' version 1.1.0 (Fiske & Chandler 2011). Results are presented as means with bootstrapped SEs.

Bayesian state-space logistic model
We modelled the population dynamics of the plain pigeon, scaly-naped pigeon and red-tailed hawk with the standard logistic equation (Runge et al. 2009, Rivera-Milán et al. 2014: where r m is the maximum intrinsic rate of population growth, K is the population carrying capacity, N t is the true unknown abundance state of the population, and M t is the total number of human-induced deaths in year t. M t = N t × m t , where m t is the mortality rate between year t and t + 1. We generated mortality rates randomly as part of the Markov chain Monte Carlo (MCMC) algorithm with uniform distributions (e.g. m ~ Uniform [0.01, 0.10] for the plain pigeon and red-tailed hawk to account for illegal hunting; and m ~ Uniform [0.05, 0.50] for the scaly-naped pigeon to account for legal and illegal hunting; Rivera-Milán et al. 2014. We re-parameterized the unknown abundance state as a proportion of population carrying capacity (i.e. N t /K) to reduce the autocorrelation of MCMC samples. We assumed that the error of state model predictions ε was lognormally distributed with mean 0 and estimated standard deviation σ p . Based on this reparameterization, we projected abundance forward in time with: We modelled the abundance proportion of the first year using a lognormal distribution with mean P 0 and variance σ 2 While the process model in the state-space formulation accounted for our incomplete understanding of pigeon and hawk population dynamics, we related true abundance state to the distance sampling abundance estimates and SEs to account for observation variance (Knape 2008). That is, we specified that abundance y t and observation variance σ 2 t,o were estimated from the surveys. Because the distribution of abundance estimates tends to be positively skewed, we assumed a lognormal distribution for σ 2 t,o (Buckland et al. 2001). We transformed abundance estimates to the natural logarithm scale by transforming the bootstrapped SEs to the standard deviations (SDs) of the corresponding lognormal distribution. To complete the observation model of the state-space formulation, we related true unknown abundance state (i.e. N t = P t × K) and estimated abundance with: where: We assumed equal mortality for all individuals of each species, regardless of age and sex class. We also assumed additive mortality from natural and humaninduced disturbances, although the model formulation allowed for compensation through linear densitydependent population growth (Runge et al. 2009, Rivera-Milán et al. 2014. We estimated unobserved population parameters and unknown abundance states under the assumption of conditional independence for each time step. We specified uniform prior distributions with wide but biologically realistic ranges. For example, in the case of the plain pigeon island-wide (see Rivera-Milán et al. 2016, their Table 1), r m ~ Uniform (0.01, 1), K ~ Uniform (10 000, 100 000), and the mean of initial abundance on the log scale P 0 ~ Uniform (−2, 0). For the SDs of process error and initial annual abundance proportion, we also used uniform priors (e.g. σ p and σ 2 t,o ~ Uniform [0, 2]).
To generate parameter posterior distributions, we ran the program JAGS version 4.3.0 (Plummer 2003) with the R package 'R2JAGS' version 0.5−7 (Su & Yajima 2015). We generated 250 000 iterations and used the first 50 000 iterations for a burn-in period. We thinned 3 Markov chains by 25 to obtain samples of 8000 iterations. We checked for convergence of the MCMC algorithm with trace plots and node summary statistics (e.g. Brooks-Gelman-Rubin diagnostic statistic R), and made posterior predictive checks of model fit with Bayesian p-values (Gelman & Hill 2007). Results are presented as means ± SDs and medians with 2.5 th and 97.5 th percentiles.

N-mixture and two-species models
The most basic Poisson-binomial mixture model fit the plain pigeon count data collected at the 247 points that were surveyed 3 times per year in the east-central region in April−June 2015−2019 (range of Pearson's χ 2 = 0.62−2.13, df = 7, p = 0.95−0.99; range of AIC c values = 709.46−837.37). More complex models (e.g. λ [pigeon + hawk density], P ij [pigeon + hawk density]; range of Pearson's χ 2 = 0.61−1.98, df = 7, p = 0.96−0.99; range of ΔAIC c values = 0.52−0.62) were informative but increased the SEs of detection and abundance estimates. As a result, there was a positive but variable association between the plain pigeon and scaly-naped pigeon point-level density estimates (range of beta coefficients β = 0.24−0.70, SE = 0.20−0.50) and detection estimates (β = 0.26−0.56, SE = 0.22−0.41), and there was a negative but variable association between the plain pigeon and red-tailed hawk point-level density estimates (β = −1.06 to −1.22, SE = 0.97−1.11) and detection estimates (β = −1.02 to −1.27, SE = 0.84− 1.10). The Poisson-binomial mixture models and dis-tance sampling hazard-rate models generated similar abundance estimates for the east-central region during 2015−2019. For example, based on the Poisson-binomial mixture models, population size averaged 543 ± 237 for the plain pigeon, 42 807 ± 1417 for the scaly-naped pigeon and 1353 ± 96 for the redtailed hawk in 2018 and 2019. For comparison with the distance sampling abundance estimates, see Tables 1−3. Based on the 2-species models, plain pigeon occupancy averaged 0.32 ± 0.03, scaly-naped pigeon occupancy averaged 0.91 ± 0.04, and red-tailed hawk occupancy averaged 0.68 ± 0.14 at the points that were surveyed in the east-central region during 2015−2019. Co-occupancy ψ AB averaged 0.48 ± 0.05 for the plain pigeon and scaly-naped pigeon, and 0.19 ± 0.08 for the plain pigeon and red-tailed hawk. That is, the species interaction factor ϕ averaged 1.66 ± 0.25 for the plain pigeon and scaly-naped pigeon, and 0.90 ± 0.43 for the plain pigeon and redtailed hawk. Plain pigeon detection probability, given both pigeon species presence, averaged 0.37 ± 0.11). Scaly-naped detection probability, given both pigeon species presence, averaged 0.69 ± 0.03. Redtailed hawk detection probability, given plain pigeon and red-tailed hawk presence, averaged 0.70 ± 0.15. Co-detection (P j AB ) averaged 0.37 ± 0.06 for the plain pigeon and scaly-naped pigeon, and 0.17 ± 0.06 for the plain pigeon and red-tailed hawk. That is, the detection interaction factor δ averaged 1.42 ± 0.49 for the 2 pigeon species, and 0.64 ± 0.33 for the plain pigeon and red-tailed hawk. The probability of occupancy and detection of the plain pigeon, given plain pigeon presence but scaly-naped pigeon absence, averaged 0.13 ± 0.06 and 0.20 ± 0.09; that is, ψ Ab and P j Ab < ψ AB and P j AB . However, the probability of occupancy and detection of the plain pigeon, given plain pigeon presence but redtailed hawk absence, averaged 0.40 ± 0.25 and 0.71 ± 0.20; that is, ψ Ab and P j Ab > ψ AB and P j AB . The association patterns were similar during 2015−2019. For example, before the 2017 hurricanes, ϕ and δ averaged 1.16 ± 0.45 and 1.42 ± 0.96 for the 2 pigeon species (i.e. positive but variable associations), and 0.89 ± 0.44 and 0.36 ± 0.19 for the plain pigeon and red-tailed hawk (i.e. negative but variable associations).

Bayesian state-space logistic model
Markov chains and node summary statistics showed convergence of the MCMC algorithm (e.g. Brooks-Gelman-Rubin diagnostic statistic, range of R = 1.00−1.01), and Bayesian pvalues showed no evidence of model lack of fit (range of p = 0.29−0.57). Posterior estimates of maximum population growth rate r m and carrying capacity K for the 2 pigeon species and red-tailed hawk are provided in Table 4. Posterior mean ± SD estimates of r m

DISCUSSION
The plain pigeon population recovered after the 1989 and 1998 hurricanes, but the most severe abundance decline occurred after the 2017 hurricanes. In contrast, the scaly-naped pigeon and red-tailed hawk populations did not decline, or recovered faster from their declines, than the plain pigeon population after the 1989, 1998 and 2017 hurricanes. That is, the scalynaped pigeon and red-tailed hawk populations showed more 'resistance' and 'resilience' (sensu Pimm 1991) than the plain pigeon population. Hugo in 1989 was a stronger hurricane but affected a smaller portion of the island than Georges did in 1998 (see Uriarte et al. 2019, their Fig. 1). However, based on the abundance decline in April to June 1990, we suggested that a hurricane with the intensity of Hugo but the path of Georges would have been de vastating for the plain pigeon population (Rivera-Milán et al. 2003a,b). In September 2017, hurricane Irma brought heavy rains and category 2 winds that affected mostly the northeastern and east-central regions. However, a hurricane with the path and intensity of María had not occurred since the 3 major hurricanes that occurred in 1899−1932 (see Boose et al. 2004, their Table  2 and Fig. 4), and the structural vegetation damage from recordbreaking rainfall and category 4 winds was extensive in lowland and highland forests across the island, including the centre of abundance of the plain pigeon in the east-central region (see Hall et al. 2020, their Figs. 1 & 3). Because of reproductive adaptations and strategies typical of the columbids (e.g. crop-milk production, multiple brooding and ex tended nesting seasons; Blockstein & Westmoreland 1993), successful reproduction is of key importance for population recovery after major hurricanes (Rivera-Milán 1996, 2001, Rivera-Milán et al. 2003b). However, contrary to the scaly-naped pigeon, which was actively calling and nesting during the 2018−2021 surveys, the plain pigeon was rarely calling and the few nesting at tempts we observed were un successful most likely due to nest predators, such as the black rat Rattus rattus and pearly-eyed thrasher Marga rops fuscatus, which were abundant in second-growth forests of the eastcentral region (Rivera-Milán 1996, 2001, Rivera-Milán et al. 2003b, Ventosa-Febles 2019. In addition, because the detection probability of the plain pigeon and scaly-naped pigeon was higher at the points with none−low than medium−high vegetation cover, the extensive structural vegetation damage caused by the 2017 hurricanes may have increased pigeon exposure and predation rates (Santana C. & Temple 1988, Ruiz-Lebrón et al. 1995, Beissinger et al. 2008, Vilella & Nimitz 2012, Vilella et al. 2013, Gallardo et al. 2019. Detection, abundance and occupancy estimates suggested positive but variable associations between the plain pigeon and scaly-naped pigeon (i.e. contrary to what we expected, ϕ and δ > 1), and negative but variable associations between the plain pigeon and red-tailed hawk (i.e. as we expected, ϕ and δ < 1). Although alternative explanations cannot be ex cluded  Table 4. Posterior estimates of maximum population growth rate r m and population carrying capacity K for the plain pigeon, scaly-naped pigeon and red-tailed hawk island-wide and in the east-central region of Puerto Rico in April to June 1986−2021 with our survey-based count data (e.g. competitive interactions from exploitation, interference or exclusion; Pérez-Rivera 1978, Poling & Hayslette 2006, Cooper et al. 2007, we suggest that a simple explanation for the positive associations is that both pigeon species responded similarly to hawk predation and foraging resources, converging at localities with medium−high vegetation cover and food availability (James & Boecklen 1984, Rivera-Milán 1995a, 1996, 2001, Rivera-Milán et al. 2003b. As a result, plain pigeon detection and occupancy estimates were higher with scaly-naped pigeon presence than ab sence at the points surveyed in the east-central region during 2015−2019 (i.e. ψ Ab and P j Ab < ψ AB and P j AB ). In addition, because both pigeon species were observed flying and feeding together on plants such as the Puerto Rico royal palm Roystonea borinquena, day-blooming jasmine Cestrum diurnum and princess vine Cissus verticillata before−after the 1989, 1998 and 2017 hurricanes, and because similar re sults were obtained with different datasets (e.g. transect surveys of columbid nests during 1986−1999; see Rivera-Milán 2001, their Fig. 2), we suggest that individual-level competitive interactions with the scaly-naped pigeon had a negligible effect on the plain pigeon population dynamics.
Regarding predator−prey interactions, because the red-tailed hawk was observed attacking the plain pigeon and scaly-naped pigeon before−after the 1989, 1998 and 2017 hurricanes, we suggest that the pigeons tended to avoid localities with none−low vegetation cover and medium−high hawk density. As a result, plain pigeon detection and occupancy estimates were higher with red-tailed hawk absence than presence at the points surveyed in the east-central region during 2015−2019 (i.e. ψ Ab and P j Ab > ψ AB and P j AB ). Although the pigeons can respond to hawk audio and visual cues with avoidance and other antipredator behaviours (Burger et al. 1989, Cook et al. 2005, Stephan & Bugnyar 2013, we suggest that additive mortality from predation and humaninduced disturbance can influence the population dynamics of both pigeon species (Rivera-Milán et al. 2014, and increase the risk of extinction of the plain pigeon at lowered abundance levels (i.e. distance sampling mean ± SE estimates N = 1043 ± 476 island-wide and N = 522 ± 157 at the centre of abundance in the east-central region in April to June 2018−2021; Table 1, Fig. 4).
In the Bayesian state-space logistic model, the parameter r m represented the exponential rate of in crease of the pigeon and hawk populations from low numbers and under favourable environmental conditions (i.e. in the absence of natural and human-induced disturbances and with plenty of food and other resources needed to maximize reproductive and recruitment rates; Niel & Lebreton 2005, Runge et al. 2009). In the model, the parameter K represented the abundance levels above which the pigeon and hawk populations tended to decline due to the effects of density dependence (e.g. the size of home ranges and limited space with suitable resources to survive and reproduce; Vilella & Nimitz 2012, Vilella et al. 2013. Although parameters r m (range of CVs = 0.08−0.63) and K (CVs = 0.08−0.49) and abundance predictions N p (CVs = 0.30−1.60) were difficult to estimate with precision (e.g. at CVs < 0.20) from survey-based count data collected annually under varying environmental conditions (i.e. future population dynamics may or may not repeat previous ones due to uncertainty from climate change, extreme weather and variation in density-independent and density-dependent responses, resulting in high process variance; Knape 2008, Knape & de Valpine 2011, Oppel et al. 2014, Saether et al. 2016, we monitored the recovery of the pigeon and hawk populations from low numbers during a period of natural reforestation of abandoned pasturelands and croplands (Koenig 1953, Brandeis & Turner 2013, Yuan et al. 2017, and before−after severe abundance de clines from the landfall of 3 major hurricanes during 1986−2021 (Boose et al. 2004, Uriarte et al. 2019, Hall et al. 2020. Because parameter posterior estimates and abundance predictions were updated annually with survey-based count data, we continued learning from monitoring and modelling, despite imperfect observations of population state (e.g. due to detection probability), partial control over human-induced disturbance (e.g. due to illegal hunting) and incomplete understanding of population dynamics (e.g. due to environmental stochasticity).
In addition, although we made a series of simplifying assumptions (e.g. linear density dependence and equal mortality probability, regardless of age and sex class), the Bayesian state-space logistic model was useful for estimation and prediction, despite differences in life-history characteristics and demographic traits (i.e. columbid early maturity with low survival and high reproductive rates as opposed to raptor delayed maturity with high survival and low reproductive rates). For example, using the demographic invariant method (see Niel & Lebreton 2005, their Eq. 15), and assuming annual survival rates were in the range of 0.40−0.70 for hatch-year and after-hatch year pigeons (Vilella et al. 2013), r m was in the range of 0.24−0.60; and assuming annual survival rates were in the range of 0.70−0.95 for after-hatch year and after-second year hawks (Gallardo et al. 2019), r m was in the range of 0.03−0.24 (for comparison, see the medians and 2.5−97.5 th percentiles in Table 4). Likewise, using life-stage survival probabilities and a seasonal matrix model, Gallardo et al. (2019) ob tained an intrinsic rate of increase r in the range of 0.05−0.24 for the red-tailed hawk. The model provided an adequate representation of the pigeon and hawk survey-based count data (i.e. Bayesian p-values were within the range of 0.1−0.9; Gelman & Hill 2007), and generated reasonable posterior estimates of K and N p , including 2020 when surveys were not conducted due to the coronavirus pande mic. The posterior estimates of K overlapped with the distance sampling abundance estimates, and have been used to set population management ob jectives (Sanderson 2006, Traill et al. 2010, Rivera-Milán et al. 2014. The posterior estimates of N p , which accounted for observation and process variance, also provided valuable information about the expected magnitude of the pigeon and hawk population fluctuations during 2020−2030. Based on our monitoring and modelling results, we suggest that the abundance of the plain pigeon, scaly-naped pigeon and red-tailed hawk increased from low numbers and approached carrying capacity levels during the late 1990s in response to increased foraging and nesting resources in second-growth forests across the island. However, while the abundance of the scaly-naped pigeon and red-tailed hawk continued fluctuating around carrying capacity, we suggest that the abundance of the plain pigeon declined during 2002−2017 as a result of habitat loss and degradation from development, which was widespread in the east-central region (see Martinuzzi et al. 2007, their Fig. 8), and the combined effects of predation and human-induced mortality (Rivera-Milán 1996, Rivera-Milán et al. 2003b. The plain pigeon has a clutch size of 1 egg instead of the typical columbid clutch size of 2 eggs (Blockstein & Westmoreland 1993), which can limit annual reproductive output and lower maximum population growth rate. The 2.5−97.5 th percentiles of r m ranged from 0.50 to 0.74 for the scaly-naped pigeon and from 0.10 to 0.56 for the plain pigeon. That is, population doubling time (t = ln[2]/r m ) ranged from 0.94 to 1.38 yr for the scaly-naped pigeon and from 1.24 to 7.30 yr for the plain pigeon. Therefore, the plain pigeon population had less 'demographic vigour ' (sensu Caughley 1977) and was less resilient (Pimm 1991) than the scaly-naped pigeon population to natural and humaninduced disturbances.
The plain pigeon population did not recover during the 2018−2021 surveys and may undergo a pro-longed bottleneck (Miyamoto et al. 1994, Beissinger et al. 2008, which may result in extinction, particularly if human-induced disturbance continues unabated, reproduction remains unsuccessful, and another hurricane with the path and intensity of María makes landfall on Puerto Rico during the current decade. Therefore, in addition to long-term monitoring, which is needed to update the posterior estimates of population parameters and increase the precision of abundance predictions, we recommend to (1) continue predator removal efforts at nesting localities to increase and maintain reproductive success rates above 0.40 ± 0.05 (SE) (Rivera-Milán et al. 2003b, Ventosa-Febles 2019), (2) use radio and satellite telemetry to estimate annual survival rates and quantify movements and resource use (Meyer & Wilmers 2007, Vilella et al. 2013, (3) combine count, telemetry and nesting data using an integrated population modelling framework (Johnson et al. 2010, Oppel et al. 2014) and (4) initiate a captive breeding programme with reintroduction methods similar to the Puerto Rican parrot Amazona vittata (White et al. 2005, Earnhardt et al. 2014, Clubb et al. 2015. We also need to minimize human-induced mortality from illegal hunting (Rivera-Milán et al. 2016); and, because most of the island forests are under private ownership (Brandeis & Turner 2013), we need to engage landowners in habitat conservation through programmes such as the Partners for Fish and Wildlife (see www.fws.gov/caribbean/ES/ and www. fws.gov/partners/). We consider the plain pigeon sub species Patagioenas inornata wetmorei a unique conservation unit (Miyamoto et al. 1994, Russello et al. 2010 in need of targeted, question-driven monitoring and evidence-based management (Sutherland et al. 2004, Nichols & Williams 2006, Lindenmayer & Likens 2010, aiming to increase and maintain abundance above the carrying capacity 2.5 th percentiles (i.e. > 30 000 island-wide and > 5000 in the east-central region; Sanderson 2006, Traill et al. 2010, Rivera-Milán et al. 2016.