{'fraction_non_alphanumeric': 0.028233271621264214, 'fraction_numerical': 0.017587939698492462, 'mean_word_length': 5.025697211155379, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Cast away on Mindoro island: lack of space limits population growth of the endangered tamaraw (Bubalus mindorensis)

Endangered species, despite often living at low population densities, may undergo unexpected density-dependent feedbacks in the case of successful recovery or marked reduction in range. Because density-dependence dynamics can increase risk of extinction, these effects can hamper conservation efforts. In this study, we analyse the dynamics of the largest population of the tamaraw (Bubalus mindorensis), a critically endangered ungulate species endemic to Mindoro island, Philippines. The population is located within a < 3,000 ha area in Mounts Iglit-Baco Natural Park, with limited expansion possibilities. We took advantage of a 21 year time series of tamaraw counts to estimate annual population growth rate and possible density-dependence, accounting for sampling errors in the counts. The tamaraw population has been increasing at an average rate of +5% per year, as would be expected given its protected status by law. Population growth showed strong spatial structuring, with a population growth close to +10% in the core area of protection, and a reduction of abundance of −5% at the periphery of its range, inside the protected area. This range constriction is concerning because our best population dynamics model suggests significant negative density-dependence (Bayes factor = 0.9). The contraction of tamaraw range is likely caused by anthropogenic pressures forcing the species to live at relatively high densities in the core zone where protection is most effective, creating source-sink dynamics. Our study highlights the fact that, despite the continuous population growth over the last two decades, the long-term viability of the Mounts Iglit-Baco Natural Park tamaraw population remains uncertain.


INTRODUCTION
Assessing population abundance is an important component of conservation planning because of the strong link between population size and likelihood of extinction (Boyce 1992). Accordingly, low population abundance is one of the primary criteria used for the evaluation of conservation status (Mace et al. 2008;Neel et al. 2012). Indeed, the effects of small population size on species persistence has been identified as one of the two central paradigms of conservation biology (Caughley 1994). Research has shown that the population dynamics of small populations differ in several ways from what is reported at higher abundance (Caughley 1994;Mugabo et al. 2013). When population size is small, demographic stochasticity -the random death of a few individuals -is more influential on the population dynamics and its viability than in populations that are larger in size (Shaffer 1981;Lande 1993). Similarly, when mating partners are too few to meet and reproduction fails, demographic stochasticity generates the so-called Allee effect (Courchamp et al. 1999), characterized by a decrease of population growth at low population abundance (i.e. demographic component, see Stephens et al. (1999)). Yet demographic stochasticity and the Allee effect are not the only considerations for the conservation of small populations; other density-dependent effects may also put these vulnerable populations at risk.
Classically, negative density-dependence describes situations where population growth rate is altered by increasing animal density (Nicholson 1933), with crowding driving the fate and viability of populations at risk of extinction (Beissinger & Westphal 1998;Henle et al. 2004). Negative density-dependent effects on population dynamics and demographic rates are often overlooked in conservation for a simple reason: they are a priori expected to occur at high population abundance, a counter-intuitive situation for threatened populations.
Yet, even in small populations, increasing abundance can occur rapidly as the result of (i) a range contraction caused by habitat loss, (ii) the reintroduction of individuals for population reinforcement, (iii) the disruption of natural emigration due to habitat fragmentation, or (iv) following successful protection measures without the possibility for the population to expand geographically. In situations where there is low carrying capacity and limited connectivity for animals to seek less crowed areas, negative density-dependence may hinder some conservation goals. Achieving long-term population growth, or the removal of animals for reinforcement are, for instance, deemed to be difficult when negative density-dependence is at work.
There are multiple well-documented density-dependence effects on vital rates (Eberhardt 1977;Fowler 1987;Bonenfant et al. 2009), most of which are the result of behaviour and life history trait changes as population density approaches carrying capacity (Bonenfant et al. 2009;Choquenot 1991; Van der Wal et al. 2014). First, with increasing population density, individuals may die from starvation because of a reduction in the availability and accessibility of food resources by conspecifics. Such antagonistic interactions may generate high stress levels which, in turns, magnify the effect of other mortality sources (Peterson & Black 1988;Hone & Clutton-Brock 2007). For instance, at high population density, epizootics are more likely to emerge and to spread quickly in populations (Langwig et al. 2012). Second, higher densities may also force non-dominant animals to occupy low-quality territories, ultimately decreasing the average productivity and survival of the population (Krüger & Lindström 2001). Given its numerous biological consequences, early detection of density-dependence is important for the management of threatened species because its evidence from data would suggest that the population is approaching the carrying capacity of the environment, and starts to be limited by food resources (Eberhardt 2002).
The tamaraw (Bubalus mindorensis) is a large ungulate endemic to the island of Mindoro, Philippines (Heude 1888;Custodio et al. 1996). The species was originally found across the entire island of Mindoro, with a total estimate of 10 000 individuals in 1900 (Long et al. 2018). Since then, tamaraw populations suffered from land conversion for agriculture and logging, trophy hunting, disease outbreaks spread by domestic cattle (Maala 2001), and traditional hunting conducted by upland indigenous communities who share their living space with the species (Long et al. 2018). Despite being officially protected by law since 1954, as well as the creation of the sizeable Mounts Iglit-Baco Natural Park (MIBNP) protected area in 1970 to protect its key habitat (Maala 2001), tamaraw range has decreased substantially over the last decades. As a consequence, the tamaraw is now found in a few isolated populations scattered across Mindoro, with MIBNP holding the largest population (Matsubayashi et al. 2010;Wilson & Mittermeier 2011;Long et al. 2018), restricted to a single location <3 000 ha hereafter referred to as the core zone of monitoring (CZM).
Given its importance for the conservation of this iconic species, the MIBNP tamaraw population has been the cornerstone of past, present, and future conservation plans (Maala 2001;Long et al. 2018). In an attempt to reinforce tamaraw protection, wildlife managers came to an agreement with residing indigenous communities in 2016, declaring a 1 600 ha no-hunting zone area within the CZM. Unfortunately, one unforeseen consequence of the no-hunting area establishment is an increase in hunting pressure at its boundaries. Park rangers regularly report the setting of snare and spear traps in the border areas, and violation of the no-hunting zone agreement is not uncommon, with several deaths of tamaraws confirmed in the last years (Boyles pers. comm.). Furthermore, field observations by rangers suggest that tamaraw population density is growing within the no-hunting zone, which could increase the risks of crowding. In this context, it is possible that negative density-dependence combined with the source-sink dynamics will affect the population in the near future, undermining long-term conservation strategies.
Here, we took advantage of a 21 year time series of abundance data to provide a detailed analysis of the tamaraw population dynamics in the CZM of MIBNP. We predict (i) that because of the substantial conservation effort in the last decade, overall population growth rate (noted r) should be statistically > 0.
Yet (ii), because the no-hunting agreement area creates a spatial demarcation between risky and safe habitats, we expect local growth rates to differ in space, and predict a decrease of r in areas further from the main rangers' base camps and patrolling routes. Finally (iii), given the intensive conservation efforts in the CZM, and relative safety of tamaraw in the CZM, we predict a negative association between r and tamaraw abundance within the no-hunting zone as a result of negative density dependence effects. To test the first two hypotheses, we computed the population growth rate of the tamaraw population at two spatial scales: over the entire CZM, and at specific locations used in annual tamaraw counts (hereafter referred to as vantage points). To test our third hypothesis of density-dependent demographic responses, we fitted and compared three population dynamics models: the baseline exponential model (no density-dependence), a Ricker model to test for a linear relationship between r and annual abundance, and a particular formulation of Gompertz model to test for a decrease of r through time.

Study site
The Mounts Iglit-Baco Natural Park (hereafter noted MIBNP) is a 106 655 ha protected area of the Philippines in the south-central part of the island of Mindoro (N12 • 54 , E121 • 13 ). The Mounts Iglit-Baco Natural Park harbors the largest population of tamaraw, found in a single 3 000 ha area located on the south-west edge of the protected area. This 3 000 ha area demarcates our study site, and is where most patrolling and observation activities are carried out (the so-called CZM). There are currently around 15 park rangers, which are divided into two or three teams to monitor the area on a regular basis. The teams patrol routes within the CZM, centered around the three permanent base camps. The CZM and no-hunting zone represent approx. 3 and 1.5% of MIBNP, respectively.
The CZM area is a rolling grassland plateau with an average elevation of 800 m a.s.l., dominated by cogon (Imperata cylindrica) and wild sugarcane (Saccharum spontaneum), and interspersed with numerous wooded creeks, secondary forest fragments, and steep hills. Introduced from Borneo to Mindoro in the early 60s (Pancho & Plucknett 1971), the invasive Siam weed (Chromolaena odorata) has been rapidly spreading in the grassland area and now covers an unknown fraction of the tamaraw habitat. The tropical climate is strongly seasonal with a rainy season from June to October (mean temperature 26.5 • C; mean rainfall: 150 mm), and a hot dry season from November to May (mean temperature 29.8 • C; mean rainfall: 5 mm). Apart from humans, the only potential predator to tamaraw is the reticulated python (Malayopython reticulatus) that is capable of preying upon calves or yearlings, though no such predation observations have been made so far.

Tamaraw counts
Every year since 2000, the Tamaraw Conservation Program (TCP), a program of the Philippine Department of Environment and Natural Resources (DENR), has conducted tamaraw counts across the CZM, yielding an annual index of relative population abundance. The tamaraw counting area consists of 18 vantage points, covering nearly 2 200 ha within the CZM (Fig. 2). The tamaraw counting process is rather invasive for the MIBNP ecosystem because it involves the burning of grasslands a few weeks ahead of each count to increase visibility and to retain tamaraw to areas with nutrient-rich young grasses. For the past 21 years, between 1 200 and 1 500 ha of grasslands were intentionally burned each year for the purpose of monitoring tamaraw abundance.
During the annual counts, individual tamaraws were counted simultaneously at the 18 vantage points for one hour and a half at dusk and dawn over a period of four consecutive days in late March or early April. For the 8 sessions and at each vantage point, at least two observers spotted and recorded all animals seen, writing down observation time, approximate location, and sex and age-category of individuals (split into calves, yearlings, sub-adults, and adults). This sampling design remained unchanged for the 21 years, though the number of observers could vary from year to year. Immediately after the completion of the 8 sessions, raw data were cross-checked and cleaned for potential multiple observations. Rangers observing at neighbouring vantage points decided together if a group of tamaraws was previously seen or not based on the timing, moving direction and composition of the group. When timing, moving direction and composition did match, then those individuals were withdrawn once from the total number of sightings. We derived the total number of different animals corrected for multiple counts for each vantage point and year, and constitute the consolidated counts C. These annual counts comprise the only available source of information about tamaraw population trends and dynamics for officials to use for conservation management decisions. Note that, because of the COVID-19 pandemic, the tamaraw counts in 2020 did not take place and we missed detailed counts (vantage point level) for years 2000 to 2003 and 2005 following the loss of field data sheets. We considered all these missing values as not available to be estimated by our statistical models (see below for details).

Data analyses
For populations of long-lived and endangered species, robust and accurate estimations of the population growth rate (r) remain difficult to obtain.
Knowledge about r requires either the estimation of demographic rates combined into a Leslie matrix (Leslie 1945) or the availability of long-term series of population abundance (Royama 1992), preferably on an annual basis. Because we were unable to mark tamaraw, both for legal and ethical reasons (Putman 1995), we could not individually identify animals, and thus it was not possible to apply methods of abundance estimations based on capture-recapture (see Schwarz & Seber 1999, for a review). We therefore used auto-regressive statistical models to analyze the annual counts. Working with abundance index data to derive an empirical estimates of r, however, comes with serious methodological issues (Îto 1972;Lindley 2003). Over the last decades, population counts have repeatedly been associated with a large sampling variance (e.g. Caughley 1977;Morellet et al. 2007), and coefficients of variation for population abundance close to 30% have been reported for various large mammal species (Lubow & Ransom 2016). From a conservation perspective, a large sampling variance in the input data for estimating r is a major pitfall because it leads to over-optimistic results, usually from an over-estimation of both the point estimate and its precision (Lindley 2003 (2002). To do so, we defined N s,t as the unobserved 'true' abundance linked to C s,t by a random effect corresponding to the process variance (σ 2 c ) on the log scale to ensure positive values for abundance. Our baseline observation model was: Accordingly, we computed the population growth rate r s,t of the tamaraw population at each vantage points s, as: from counts corrected for the sampling variance, thus returning an unbiased estimation of the population growth rate provided there is no long-term trend in the detection probability of tamaraws. The overall population growth of the tamaraw population of MIBNP is then given by: To test our third hypothesis we investigated potential density-dependence in the annual growth rate by fitting two models classically used to detect negative density dependence (Dennis & Taper 1994): the Ricker and a Gompertz growth model, to be compared with an exponential growth model. The exponential growth model was fitted implicitly for the estimation of the average population growth rater and served as our baseline model: Then we fitted a Ricker growth model to the total number of tamaraws seen per year N t as follows: where r m is the maximum population growth for the tamaraw, and K l is the Because the pair of parameters r m and K l (Ricker) and β and K G (Gompertz) are not separately identifiable (Lebreton & Gimenez 2013), we used informative priors for r m and β based on the reported maximum population growth for other large bovids (see Traill et al. 2007, for a similar approach). Accordingly, we set a moderately informative prior value of mean 0.3 and a precision of 0.1 in our models (Table 1). We tested the occurrence of density-dependence processes in the MIBNP tamaraw population by computing a numerical approximation of the Bayes factor. We used MCMC samples to calculate the probability (π k ) for our abundance time series to be generated by each of the three candidate demographic models during the model fit. At each MCMC step, a multinomial variable stated which of the exponential, Gompertz and Ricker model was the most likely underlying model given the data by taking the value one and zero otherwise. We estimated the three probabilities π k as the number of times each models received a one divided by the number of MCMC iterations.
We fitted our statistical model to the tamaraw count data with JAGS 4.0 (Plummer et al. 2003), with 3 MCMC chains and a burn-in of 40,000 iterations.
We obtained parameter estimates with an additional 30 000 iterations and a thinning factor of 5 (therefore giving a total of 6 000 MCMC samples). We Note that despite the fact that we analysed the number of tamaraws as described above, we chose not to display such numbers on figures and to report a standardized abundance (C t = C t /max(C t )) instead. This transformation applies to the estimated carrying capacities of the Gompertz and Ricker models, We did not provide raw counts to avoid confusion between the index of relative tamaraw abundance we worked with, and a real population size estimator that was not implemented at MIBNP. With this variable transformation K l and K G are unitless and can no longer be interpreted as the maximum number of tamaraws but the maximum relative abundance (C t ) we can expect to observe at MIBNP. Reported values of the carrying capacities K l and K g can hence take values > 1 because the carrying capacity has not yet been reached. Biologically speaking, we put more emphasis on how far from the carrying capacity the current abundance of tamaraws is (C t /K) than of the number on animals per se. with an associated variance process equal toσ proc = 0.25 0.31 0.38 (Fig. 1). It means the observed abundance increased 2.41 times over the course of the study.

From
As expected for count data, the sampling variance was largeσ c = 7.18 9.21 9.98 , suggesting frequent double counts of the same individuals. For the vantage point specific growth rates (r s ), we observed a marked variability in space ( Fig. 2; see Supplementary Table S1) as shown by the fact that growth rates ranged between We explicitly estimated the probabilities π k that the exponential, Ricker, or Gompertz models would be the underlying process generating the observed variation in tamaraw abundance (proxy of Bayes factor). The estimated probabilities for the three population dynamics models given our count data werê π 1 = 0.10,π 2 = 0.90 andπ 3 = 0.00, respectively, supporting our hypothesis of density-dependence feedbacks in the tamaraw population. Point estimates for the two parameters of the Ricker model werer m = 0.07 0.08 0.09 for the maximum population growth rate andK l = 1.16 1.37 1.62 for the carrying capacity.
Although the Gompertz model received no statistical support, the estimated coefficients were statistically significant and readβ = 0.04 0.05 0.07 , and K G = 1.20 1.41 1.66 for the carrying capacity. Note that the statistical support for the Ricker model held for the dynamics of the core area only, that is, when considering only the 9 vantage points inside the no-hunting area (π 1 = 0.02, π 2 = 0.83 andπ 3 = 0.15).

DISCUSSION
The isolated tamaraw population in MIBNP has been increasing in size since 2000 (Fig. 1). Undoubtedly, a positive growth rate of 0.05 over 21 years of monitoring is a noteworthy success for such a highly endangered species, both for local populations and for global conservation efforts. However, the devil lies in the details, and two additional results darken this otherwise positive picture.
First, we found support for a source-sink dynamic of the tamaraw population linked to poaching and evidenced by a clear spatially structured pattern in the local population growth rates. Second, our results highlight a progressive decrease of the average growth rate in time, suggesting density-dependent effects at the population level. The recent population dynamics of tamaraw we document at Mindoro suggests that new conservation actions might be needed to secure the future of the species.

Long-term Positive Growth of Tamaraws
The average population growth rate ofr = 0.05 -corresponding to a 2.5 fold increase in abundance between 2000 and 2021 -gives strong evidence that successful conservation policies over the past two decades have led to an increase in tamaraw abundance. Nonetheless, an annual growth of 5% is relatively low compared to other similar-sized or even larger cattle species (

Possible Impacts of Human Disturbance and Hunting
While the positive growth of the tamaraw population over the past two decades is evidence of a successful decrease in hunting pressure in the CZM, human activities could also account for the low tamaraw population growth rate we report. Tamaraw have always been traditionally hunted by indigenous people in Mindoro, and despite the "no-hunting agreement" established in 2016, both intentional and unintentional tamaraw deaths still occur. Moreover, the tamaraw is sought after for bushmeat by lowland poachers. Because illegal activities are difficult to control in MIBNP, it has likely led to increasing the mortality of all age-classes of tamaraws to an unknown, but perhaps considerable extent.
In the early years of the tamaraw monitoring, the apparent population growth The marked spatial heterogeneity we document, with negative growth rates at the periphery of the counting area, could indeed indicate a lower protection efficiency in areas further away from the main patrolling routes and ranger base camps. Such a result supports our second hypothesis, suggesting an unexpected contraction of the tamaraw range at MIBNP, and the disappearance of individuals at the periphery of the CZM. In other words, the CZM area may function like a source-sink system (Pulliam 1988), with the more central and better patrolled areas exhibiting positive growth rates, and more peripheral areas with higher hunting pressure serving as population sinks. Again, our findings provide some support for this idea, as local growth rates in the core area of the no-hunting zone have values as high as 14% (Fig. 2), which is more in line with the average growth rates that we reviewed for large bovids (Table 1).
Source-sink dynamics were found in populations of apex predators (wolves

Space Needed for Sustaining Population Growth
The gradual decrease of the tamaraw population growth rate in MIBNP, coupled with the overall abundance in the CZM, suggests a type of density-dependence is occurring, where a change in demographic rates lowers the number of individuals added to the population from year to year (Fowler 1987;Sinclair 1989). Predicted growth rates from the lowest to the highest abundance varied between r = 0.08 and r = 0.01 for tamaraw. From the current data, our best model suggests a carrying capacity of 1.4, meaning that only 40% of growth is left possible before reaching the carrying capacity. If competition for food or mating becomes strong then the main tamaraw population may experience an increased net loss of animals in the near future. Tamaraws may already be experiencing a food resource shortage, which could make the population more prone to parasitism (Stien et al. 2002) or epizootics (Davidson & Doster 1997).
More generally, density-dependence should be of concerned and investigated for large herbivores living in protected environments lacking apex predators, which contribute to limit abundance and reduce annual growth of prey species (Ripple & Beschta 2012;Sinclair et al. 2003).
Another behaviour associated with high population abundance is dispersal (Matthysen 2005). In mammals, sub-adult individuals are more likely to disperse at high than low population abundance (e.g. Loe et al. 2009, for an example on red deer Cervus elaphus). At the high population abundance we documented in MIBNP, we would expect a greater number of dispersing individuals to venture outside of the no-hunting zone, where the risk of mortality is high because of hunting or poaching. In addition to jeopardizing the fate of the population by decreasing survival probabilities of sub-adults that are key to maintaining the future of the population, the empirical evidence for density-dependence that we report would limit the ability of this population to act as a source. With a combination of low recruitment and low survival of sub-adults, old individuals would quickly become over-represented in the population, which could depress growth even further.
Unfortunately, in the absence of individual monitoring of tamaraws to estimate demographic rates independently of abundance, the real cause of the density-dependence we report is uncertain. We are currently unable to determine if the gradual decrease in growth rate is the consequence of an increased mortality or dispersal rate, a reduced reproduction, or some combination of these factors. The difficulty of teasing apart different demographic processes is a well-known limit of the so-called pattern-oriented approach (sensu Krebs 2002) when working with time series of abundance. Further investigation of these possibilities is needed, as each situation would require different management approaches.

Implications for Conservation
The main conservation goals for the tamaraw are (1) to ensure its long-term viability, and (2) support the expansion of the population within MIBNP and potentially in other sites. Current conservation measures assume that tamaraw population growth will continue in the next decade. The detection of density-dependence, however, suggests continued population growth will not be the case. In this situation, the dynamics of the tamaraw will change in the years to come if resources become the main limiting factor inside the no-hunting zone, forcing animals to disperse to, and potentially get killed in the outer zone.
Several actions must be implemented by managers and authorities to address this issue in the short term. The priorities should be to expand the no-hunting zone or to create safety corridors through consensual agreement with indigenous people, both implemented in the context of creating a more sustainable land-use and traditional hunting system. Arriving at such an agreement with indigenous people will be a long process, but these conversations are already ongoing.
Extending the no-hunting area will ultimately increase the sustainable population size, thus making it more viable and maintaining more genetic diversity.
Likewise, allowing for dispersing individuals to leave the CZM and possibly settle in distant areas would relieve density dependence pressure on the core population.
However, one could argue that such conservation actions only represent a short-term solution, as even if the size of the no-hunting area was increased, the tamaraw population would eventually grow towards the carrying capacity. We acknowledge this, but believe that these short-term conservation actions are still important, as they would help buy valuable time that is necessary to devise more long-term strategies. Any long-term strategy will likely involve active population proposition, but a preliminary feasibility study for translocation and captive breeding is a realistic option and an appropriate next step. Ultimately, however, preventing the extinction of tamaraw will require decision makers to balance short-term actions with long-term strategies that involve active population management.

Conclusion
Securing large areas of habitat is often an effective conservation approach, but may not always be a viable strategy to protect highly threatened species, as illustrated by our study on the tamaraw population in MIBNP. More generally, we show that positive growth rates could potentially hide a less optimistic picture for some threatened species. First, the growth of a population that is confined to a small conservation area may reflect a range contraction as a result of anthropogenic pressures, rather than the effect of increased abundance per se.
Second, because of this, the population could face important demographic challenges that, in the long-term, lead to unusual problems for conservation. In this context, there is an increasing need to understand the interaction between anthropogenic activities and the population dynamics of threatened species. Our results suggest that, in some situations, conservation strategies designed to protect small, isolated populations may lead to long-term demographic issues that undermine species recovery. On the positive side, our findings also show how detailed analyses of relative population abundance data, even without absolute population size estimators, can help provide information that assists in the development of conservation measures that are essential for large mammal conservation.

ACKNOWLEDGMENTS
We are grateful to Alvaro Gonzalez for his contribution to the fieldwork and his insights on tamaraw distribution at MIBNP. We also acknowledge the invaluable help of all rangers, students and volunteers who contributed to the annual tamaraw counts. Thanks