First integrative trend analysis for a great ape species in Borneo

For many threatened species the rate and drivers of population decline are difficult to assess accurately: species’ surveys are typically restricted to small geographic areas, are conducted over short time periods, and employ a wide range of survey protocols. We addressed methodological challenges for assessing change in the abundance of an endangered species. We applied novel methods for integrating field and interview survey data for the critically endangered Bornean orangutan (Pongo pygmaeus), allowing a deeper understanding of the species’ persistence through time. Our analysis revealed that Bornean orangutan populations have declined at a rate of 25% over the last 10 years. Survival rates of the species are lowest in areas with intermediate rainfall, where complex interrelations between soil fertility, agricultural productivity, and human settlement patterns influence persistence. These areas also have highest threats from human-wildlife conflict. Survival rates are further positively associated with forest extent, but are lower in areas where surrounding forest has been recently converted to industrial agriculture. Our study highlights the urgency of determining specific management interventions needed in different locations to counter the trend of decline and its associated drivers.

with proximity to newly converted forest to industrial agriculture [31][32][33] . The tendency of village communities to hunt orangutans for bushmeat was found to be driven by complex socio-economic circumstances. Hunting tends to increase with a decrease in forest cover surrounding the village and an increase in area for agriculture in the village but a decrease in income from this sector 32,33 . The proportion of Muslim populations was also found to represent a religious constraint on orangutan hunting for meat consumption 32,36 .
Because of the challenges associated with surveying and modelling the population trends and drivers of population change of Bornean orangutans (or other species), we developed a dynamic abundance modelling methodology. Our integrated dynamic population model was applied within a hierarchical Bayesian framework 37 and can (a) project the density of orangutans based on nest counts, (b) simultaneously integrate multiple types of data (i.e. nest counts from ground and aerial line transect surveys, presence-absence data from line transects and targeted surveys, and observations from interview surveys), and (c) explicitly account for the detection error inherent in each survey methodology due to associated effort. Using this novel approach we assessed the abundance and distribution of the Bornean orangutan through time and determined the contribution of climate and land use dynamics to the changes observed.

Results
Model diagnostics and performances. Prior to fitting the model to the data, we tested for correlation among the original (unstandardized) variables and among the standardized environmental variables explaining the initial abundance, occupancy and survival rates, and found weak correlations among these variables (with absolute Pearson correlation <0.45, see Supplementary Table 1). The WinBUGS simulation converged well, as confirmed by the value of Rhat (ranged between 1 and 1.1) for all parameters, and the absence of seasonality within each Markov chain Monte Carlo (MCMC) chain plot and overlap between the three chains (Supplementary Figure 1). We also detected no apparent correlations between the posterior distributions of the coefficients of the linear and the quadratic terms for altitude (ALT), the longterm mean monthly rainfall during the dry season (DRY) and wet season (WET) (Supplementary Figure 2), which suggests the reliability of the estimated coefficients obtained for these variables. Our dynamic abundance model performed well with a good correspondence between the simulated nest predictions and the actual observations. The average Pearson correlation coefficient for all time periods is r = 0.828 (with r 1997-2002 = 0.824, r 2003-2008 = 0.818 and r 2009-2015 = 0.841) and the average R 2 is 0.804. The model also has a good correspondence between the simulated orangutan presence-absence and the actual observation obtained from interview surveys, with Sensitivity SN = 0.812 and Specifity SP = 0.726. Survey specific parameters. The probability of detecting orangutan nests from field surveys per km 2 varied depending on respective survey protocol ( Table 2). Aerial transects surveys had the highest probability of detecting orangutan nests (logit(1.516) −1 = 82%), followed by the ground transect surveys (75%). This could be because aerial surveys were usually conducted in areas with prior knowledge of orangutan occurrences due to the cost of operating the helicopter. The occurrence data of the combined aerial and ground line transects and other targeted surveys had a lower probability of detecting the nests (64%).
The probability of detecting orangutans via interview surveys was 15% on average if the respondent entered the forest less than once per month and 21% if they entered the forest more frequently ( Table 2). The reason for low detection rates of orangutans from interview survey, in comparison to the nests from field survey, is twofold: (1) orangutans are much less common than their nests, and (2) nest count surveys are generally targeted at areas with prior knowledge of orangutan occurrences due to cost constraints.
Nest decay rate was estimated to be 228 days on average for Borneo (Table 2). This however varied slightly across different forest types, where mangrove forest had the longest time to decay (266 days), followed by lowland forest (244 days), montane forest (236 days), and peat forest (209 days).
Orangutan abundance by region and land use. The dynamic abundance model estimated that the density of Bornean orangutans has declined by 25% over the last ten years (Fig. 1a). We estimated the overall density of orangutans over Borneo in the period 1997-2002 was about 15 individuals per 100 km 2 , but the density was reduced to 10 individuals per 100 km 2 in 2009-2015 (Supplementary Table 2). We estimated that Central Kalimantan had the highest density of orangutans during 1997-2015, followed by Sabah, West Kalimantan, East Kalimantan, Sarawak, and North Kalimantan (Fig. 1b and Supplementary Table 2).
The distribution of orangutan populations across different land uses varied across regions. In Sabah and Sarawak, most of the orangutan populations resided within the boundaries of protected areas (PA) and logging concessions on natural forests (LOGG) (Fig. 2). In Kalimantan, the population generally resided within the  Table 2. Posterior means and the 95% credible interval (CI) of the mean for each parameter explaining the latent orangutan population density (first level of the orangutan dynamic abundance model), the observed orangutan occurrence (second level of the orangutan dynamic abundance model), latent orangutan nest density (third level), and the observed orangutan nest density and occurrence (fourth level).
boundaries of PA and LOGG and in areas without concessions (or classified as 'OTHER'). Across the whole of Borneo, the proportion of orangutans residing within the boundary of PA has increased through time (Fig. 2), mainly because the orangutan populations have gradually disappeared from other land uses and/or the extent of PA had increased recently 9, 10, 20 , e.g. with the establishment of the Sebangau National Park in Central Kalimantan,  Drivers of changes in orangutan abundance. The initial abundance of orangutans per km 2 was most strongly associated with the amount of rainfall during both wet (WET) and dry seasons (DRY), with the greatest abundance observed in areas of intermediate rainfall during each season ( Table 2 and Fig. 3a). Survival rates also correlated most strongly with the amount of rainfall during both wet (WET) and dry seasons (DRY), however, with the rates being lowest in areas of intermediate seasonal rainfall (Table 2 and Fig. 3b). Natural forest extent (FR) was positively associated with the initial abundance and survival rates. The interactions between natural forest extent and distance to forest recently converted to industrial agriculture (FR × CFA) was positively associated with survival rates, suggesting that survival rates are lowest in areas with fragmented forest and near to new areas of industrial agriculture, as the possibility of human-orangutan conflicts increase (Table 2 and Fig. 3c). Survival rates are also positively associated with proximity to protected areas (DPA), indicating that protected areas are mitigating some threats to orangutans ( Table 2).
Based on variables explaining survival rates, we assessed drivers of orangutan population decline during 1997-2015 in each region, and this includes habitat loss, human-orangutan conflicts, anthropogenic activities, and habitat fragmentation (Fig. 4). For Sabah, we estimated that orangutan population decline is driven by (1) moderate rates of habitat loss within the boundaries of LOGG, and (2) high levels of habitat fragmentation. For Sarawak, the decline is mainly driven by (1) moderate rates of habitat loss within the boundaries of LOGG, and (2) moderate anthropogenic pressure within the boundaries of LOGG and OTHER. For East and North Kalimantan, orangutan population declines were mainly driven by (1) moderate to high rates of habitat loss and (2) moderate to high intensities of human-orangutan conflicts within the boundaries of oil palm plantation concessions (OPP) and OTHER. For West and Central Kalimantan, drivers of decline include (1) moderate to high

Discussion
Our analysis is the first robust population trend analysis for orangutans or other great ape species that includes quantitative assessments of drivers of change. Methodological challenges associated with determining spatial and temporal variation in ape density across large areas have so far made such studies infeasible, but our novel approach has overcome these challenges. Our analysis advances current estimates by providing the underlying population trend through time, with the species estimated to have declined at an alarming rate of 25% over the past 10 years. This contradicts crude population estimates proposed by different authors that have indicated an increasing number of orangutans across the island, reflecting increasingly available data on the species and associated survey efforts and not an absolute increase of orangutans (Supplementary Figure 3). This is mainly because the previous studies were conducted separately for each time period and they failed to take into account the dynamic process affecting the orangutan population change.

Orangutan abundance and competition from humans in area with intermediate rainfall.
Our model indicates that the long-term abundance of orangutans per km 2 is strongly determined by seasonal rainfall, with the species being most abundant in areas receiving intermediate rainfall during the dry season (150-250 mm per month from May to September) (Fig. 3a) and the wet season (200-400 mm per month from November to March). This is comparable to Indonesian agro-climatic zone B with 7-9 consecutive wet months (>200 mm per month) and around two consecutive dry months (<200 mm per month) 27 . This area essentially receives the right amount of rain throughout the year and is likely able to support plenty of wild tropical fruits essential for orangutans, such as Moraceae (figs) and Anacardiaceae (mangos) 38,39 . The extent of the intermediate rainfall zone on Borneo is smaller than the extent of lowlands with altitude <500 m above sea level (Supplementary Figure 4), the range that has long been recognized as the primary niche for Bornean orangutans 11 . The extent is also smaller than the area of low-moderate mean annual rainfall (Supplementary Figure 4) recently suggested by Wich et al. 18   occurs in Central and West Kalimantan, the two provinces with currently the largest orangutan populations outside protected areas. In Sarawak, the zone of intermediate rainfall also occurs around the Batang Ai National Park and Lanjak Entimau Wildlife Reserve, where most of the orangutan populations in this state currently reside.
Besides being important for orangutans, areas with intermediate rainfall are also important for people. The climate in this zone optimally supports plant productivity and agriculture, allowing year-round cultivation of crops, fruits and vegetables 27 . This is supported by the fact that the proportion of agricultural areas, i.e. plantations and agriculture fields and shrublands from abandoned agriculture, outside the government-sanctioned protected areas on Borneo, increases as they are located closer to zones with intermediate rainfall ( Supplementary  Figure 5d). Because orangutans and humans favor the same climate zone and range, orangutans are facing severe competition from humans, as confirmed by our model where the species survival rates were lowest in this zone (Fig. 3b). In this study we were able to include both altitude and rainfall seasonal pattern as predictors explaining abundance and survival rates because there are no strong correlations between these variables (Supplementary Table 1). Altitude (and its quadratic term) by itself was found to be a non-significant predictor, suggesting that altitude indirectly affects orangutan abundance and survival rates, most likely through rainfall.
While the relationship between rainfall and orangutan abundance is relatively easy to understand from the direct impact of intermediate rainfall on the abundance of wild fruits, the connections between rainfall and orangutan survival rates are more difficult to discern and are most likely related to multifaceted consequences of changing rainfall patterns as part of global climate change and anthropogenic land use change in this area, i.e. vast conversion of forest to agriculture 35,40,41 . Forest clearing has led to the loss of orangutan habitat, as well as the loss of livelihood for communities who greatly depend on forest goods. As climate becomes more erratic, periods of wild fruit scarcity may have increased and the intensity and frequency of forest fires (often originating in drained peat swamp areas) and flooding events (due to upstream deforestation) also increased 42,43 . These severe environmental circumstances have most likely led to increased competition between humans and orangutans 20 . Displaced communities who cannot generate sufficient income from agriculture may seek other income opportunities such as hunting and poaching, or are more sensitive to conflicts with orangutans over crop-raiding 44 .
The link between areas with intermediate rainfall and hunting propensity can be explained in light of recent research, suggesting that hunting tends to increase with a decrease in forest cover surrounding settlements and an increase in area for agriculture around settlements but a decrease in income from this sector 32,33 . Based on population census and land cover data among administrative districts in Kalimantan, we found that districts located within the intermediate rainfall zone have the socio-economic features that lead to higher propensities of hunting compared to districts located outside these zones. The proportions of agricultural areas outside the government-sanctioned protected areas are generally higher in districts where large proportions of these areas overlap with intermediate rainfall range (Supplementary Figure 5a). As anticipated, the proportion of forest areas within the same zones is generally lower in these districts (Supplementary Figure 5b Figure 5d). Despite being agriculturally rich, however, the percentage of people living in poverty is generally higher in these districts that derive lots of their income from industrial-scale agriculture (Supplementary Figure 5e). Also, the poverty-gap index is higher in these agriculturally rich districts (Supplementary Figure 5f), indicating that profits from agricultural development accrue to a small section of society. This indicates that the current orangutan hunting activities could be exacerbated by social and economic circumstances with displaced orangutans competing with small-holder farmers that have less and less land for their own agricultural activities. The connection between socio-economic background, particularly poverty, and hunting and poaching, is generally well known based on various case studies from Asia and Africa 39,45 . However, the evidence for claims around poverty as a driver of hunting is weak, mainly because hunting has been overwhelmingly framed exclusively as an issue of conservation and biodiversity loss rather than of poverty and development 46 , but that does not mean that poverty is not an important factor.
Recent studies have also found that hunting tends to increase with a decrease of Muslim populations in the village, suggesting that religious affiliation potentially provides a barrier to current orangutan hunting 39,42 . Based on census data, we found that agriculturally rich districts located within the intermediate rainfall range in Kalimantan generally have a large proportion of non-Muslim people (Supplementary Figure 5g). This is likely because the high agricultural value has long made these areas the primary home for large indigenous communities, most of which are non-Muslims. Thus, a low proportion of Muslim populations is likely confounded within an area's high agricultural value, without necessarily influencing the propensity to hunting and orangutan survival rates. Furthermore, our model found a minimal impact of the percentage of Muslims within districts on orangutan survival, suggesting a weak correlation between religious affiliation per se and orangutan survival rates. Furthermore, earlier study suggests that hunting for bushmeat is not solely carried out by non-Muslims for their own consumptions, but also by various communities for selling the meat 39 , implying that the current hunting practices are also driven by economic incentives such as trade. To inform suitable strategies for abating orangutan hunting requires a better understanding of individual hunter motivations, and the anthropological and economic motives driving them 47 .
Increased contact with humans may also increase the risk of infectious disease in orangutans, which can affect the survival rates of the species in the wild. Previous serological studies suggest that exposure to human pathogens does occur both in free-ranging and semi-captive orangutans 48 . Pathogens, such as intestinal parasites, can be transmitted directly from humans 49 . In rehabilitation centers, overcrowding, abnormality in the population social structure, and dietary imbalances, can exacerbate disease transmission among orangutans 48 . Forest, conversion to industrial agriculture, and climate change. Our model indicates that the long-term abundance of orangutans and survival rates per km 2 are strongly determined by the extent of natural forest. This suggests that the reduction of forest extent alone will decrease orangutan survival rates. The loss of natural forest was found to be an equally important driver of orangutan declines across all regions of Borneo during 1997-2015 (Fig. 4).
When threats from forest clearing are absent, such as in the case of populations within the boundary of protected areas, survival rates can also decline due to decreasing forest carrying capacity, e.g. increased period of wild fruit scarcity due to climate change. Both global climate change, and climatic changes directly driven by deforestation are predicted to impact rainfall patterns on Borneo, with some areas anticipated to experience significant rainfall reductions, such as prolonged consecutive dry months 50 . Isolated forest patches of orangutan habitats are particularly prone to extinction due to this type of disturbance. This is exactly the issue currently faced by orangutan populations in Sabah. Comparison among orangutan habitat networks across different regions of Borneo shows that the average size of forest patches where orangutans currently reside are lowest for Sabah (Supplementary Figure 6a) and the distance between forest patches is also largest for this region ( Supplementary  Figure 6b), suggesting that the populations in this state face the highest risks due to habitat fragmentation (Fig. 4). Hence, although large proportions of orangutan populations in Sabah currently reside within the boundary of PAs, threats from global climate change and other disturbance such as disease, as described earlier, can potentially annihilate orangutan populations within a PA due to relatively small PA size and lack of connectivity among orangutans' habitats within the current PA networks 51,52 .
Our model also found that survival rates were determined by the interaction between forest extent and proximity to forest recently converted to industrial agriculture. This is likely to be directly related to the increased possibility of human-orangutan conflicts, such as crop-raiding, over newly established large-scale industrial agriculture and hence killing of crop-raiding individuals 53 . However, we also found that survival rates increase with proximity to PAs, indicating that forest protection is mitigating some threats to orangutans. Human-orangutan conflicts during 1997-2015 were found to be equally important drivers of orangutan declines across all regions of Borneo (Fig. 4). Although conflicts due to conversion of forest to industrial agriculture appear to occur most intensively in West and Central Kalimantan compared to other regions 31 , this is probably because large orangutan populations are found in these provinces, and thus does not necessarily imply that conflicts have a relatively minimal impact on populations in other regions.
Here, we addressed human-orangutan conflicts by assessing the interaction between forest cover and proximity to forest that has been recently converted to industrial agriculture. Conflicts become less frequent with time either because orangutans become less common or adapt to the new landscape 54 . This is what likely happened in extensive areas of lowland forests in Sabah that had high densities of orangutans prior to the 1960s when the forests were converted to oil palm. However, we did not take into account the possibility that the frequency of conflicts may also vary depending on fruit scarcity. As rainfall is predicted to be more extreme in the future, increased periods of wild fruit shortages are anticipated and this could potentially affect orangutan crop-raiding behavior.

Conclusion
Orangutan populations on Borneo have declined at a rate of 25% over the last 10 years. Pressure on orangutan populations in the same period of time varied substantially among regions, with the populations in Sabah, Sarawak, East and North Kalimantan experiencing relatively moderate pressure, as opposed to high pressure in West and Central Kalimantan. The co-occurrence of orangutan populations with areas most suitable for human activities has led to an enhanced risk of human-wildlife conflicts. Unless threats from climate change, land use change and other anthropogenic pressure are abated, we predict that most populations of the Bornean orangutan will be severely impacted by human activities.
Poor connectivity among orangutan habitats between the boundaries of PAs is currently the predominant threat to orangutan populations in Sabah. Orangutan populations in Sarawak, East and North Kalimantan face the same threats as West and Central Kalimantan due to habitat loss from continuing forest conversion to industrial agriculture and human-orangutan conflicts, but the latter two areas also suffer additionally from anthropogenic activities.
As the populations in different regions face different threats, specific abatement plans should be implemented to ensure the long-term persistence of the species. This includes (1) maintaining high forest cover in orangutan habitats and improving the connectivity among the remaining forest patches where orangutans live through better spatial planning for all regions of Borneo, (2) close cooperation with plantation companies, smallholder farmers and wider communities in managing conflicts with orangutans in Kalimantan, and specifically in West and Central Kalimantan (3) improving the effectiveness of anti-hunting efforts and education and (4) developing a better understanding of the underlying socio-economic motivations of hunting.

Methods
Study area. Borneo is the third-largest island in the world (approximately 740,000 km 2 ) and is shared by the Malaysian states of Sabah and Sarawak and the sultanate nation of Brunei in the north, and by Indonesian provinces in the south (i.e. West, Central, South and East Kalimantan; the latter was recently divided to establish North Kalimantan province) (Fig. 5a). The island is largely mountainous, with mountains branching westward from the central core along the border between Sarawak and West Kalimantan, and a discontinuous series of mountain ranges running parallel to the east and southeast coasts of the island) (Fig. 5a). Borneo's interior is largely mountainous but extensive lowlands and swamps occur along the coasts. A large part of Borneo is drained by navigable rivers, which represent the principal and sometimes only routes for trade and commerce, but also present barriers to orangutan dispersal 55 We divided Borneo into grid cells with a spatial resolution of 1 × 1 km 2 , and excluded Brunei and South Kalimantan as they are outside the known orangutan range. This resolution allows us to simulate orangutan dispersal from each focal cell (100 ha) to eight neighboring grid cells, resulting in a 3 × 3 km 2 dispersal block (900 ha). This resolution conforms roughly to the home ranges of female Bornean orangutans, which vary between 150 and 850 ha 57 .
Orangutan data. We utilized two types of orangutan data: nest counts and presence-absence data. The nest count data were obtained from line transect surveys (aerial and ground) (Fig. 5b). The presence-absence data were derived from two survey approaches: (1) line transect (aerial and ground) and targeted surveys of nest observations, and (2) interview surveys of direct orangutan sightings (Fig. 5b). For each survey method, we divided the data into three time periods: (1)  The aerial survey data mainly cover Sabah and were collected between 1999 and 2012 using helicopters following different flight routes, as described in Ancrenaz et al. 8,9 , giving a total route length of approximately 2,200 km. The ground surveys were carried out sporadically between 1997 and 2015 across Borneo by various orangutan research teams and non-governmental organizations, giving a total transect length of approximately 1,200 km. The targeted surveys mainly include the reconnaissance walks, i.e. a walk following a predetermined direction through the survey area. These surveys followed a standard established methodology to detect and record the nests of great apes 3 .
To facilitate the use of nest count data collected from various methods of line transect surveys, we standardized the metric of orangutan nests to obtain a nest density estimate for each 1 × 1 km 2 grid cell. For the ground surveys, we calculated the density of orangutan nests using the Distance sampling method, based on the perpendicular distance of each nest to the transect 12 . For the aerial surveys, the data were mainly in the form of an aerial index value (AI) describing the number of nests detected per km of flight. Following Ancrenaz et al. 8 , the density of orangutan nests per km 2 , i.e. gnest, can be estimated via: log(gnest) = 4.7297 + 0.9796 log(AI). Density estimates for each 1 × 1 km 2 grid cell were then obtained by averaging the estimate across all aerial surveys conducted within the grid cell, giving approximately 6,500 of 1 × 1 km 2 grid cells where orangutan nest surveys had been conducted across Borneo. These data were then used to form a matrix array of orangutan nest density Y i,j,t comprising three matrices of survey period (t), with each matrix consisting of 6,500 rows of grid cells (i) and 2 columns of survey protocol (j), i.e. ground and aerial transects.
To derive the occupancy of nests in each 1 × 1 km 2 grid cell from the ground and aerial transect and targeted surveys for each time period, we first divided the grid into sub-cells with the resolution of 200 × 200 m 2 . This is to avoid duplicated reports of the same clusters of nests. If at least one survey reported the occurrence of a nest within a sub-cell, we defined that orangutan nests were observed in this sub-cell. If no orangutan nests were recorded within the sub-cell in any of the surveys, we defined that orangutan nests were unobserved in this sub-cell. We then constructed a matrix array Znest i,k,t comprising three matrices of survey period (t), and with each matrix comprising 6,500 rows of grid cells (i) and 25 columns of nest observations within sub-cells (k).
The interview surveys of orangutan sightings were conducted in 540 villages across Kalimantan and Sabah in 2008 and 2009, and verification surveys in 2011, with 10 respondents in each village, as described in Meijaard et al. 10 . Each respondent was asked how frequently he or she entered the forest around the village (i.e. more than once per month or less than once per month) and the last time they had seen an orangutan either in the forest or in the village (i.e. within this year or more than a year ago). Additionally, personal details of each respondent were recorded, including their age and how long they had resided in the village. Based on this information, we derived the occurrence (observed or unobserved) of orangutans in each 1 × 1 km 2 grid cell and constructed a matrix array Zou i,m,t , comprising three matrices of survey period (t) with each matrix consisting of 540 rows of grid cells (i) and 10 columns of respondent observations (m). Because the chance of any respondent sighting an orangutan would likely depend on that respondent's frequency of entering the forest, we also constructed a corresponding binary matrix FE i,m , coded as '1' when respondent m entered the forest around the village in grid cell i more than once a month and '0' when less than once a month.
Dynamic abundance model. The model. We adapted a dynamic population model developed by Chandler & Clark 37 for integrating count data and presence-absence data of a species. Our model generalizes the negative binomial model for open populations and assumes that abundance patterns are determined by an initial territory establishment process followed by gains and losses resulting from births, mortalities and dispersal. It also accounts for varying detection errors inherited from different survey data. Our model requires both spatial and temporal data and consists of four broad levels: (1) latent orangutan population density, (2) observed orangutan occurrence, (3) latent orangutan nest density, and (4) observed orangutan nest density and occurrence. The first level (latent orangutan population density) can be described as: The second level (observed orangutan occurrence) as: The third level (latent orangutan nest density) as: Finally, the fourth level (observed orangutan nest density and occupancy) as: O i,t is the latent occurrence of orangutan at grid cell i in survey period t, Nou i,t is the latent number of orangutans at grid cell i in survey period t, S i,t is the latent number of survivors at grid cell i that do not emigrate between period t and t + 1, R i,t is the latent number of recruits (including births and immigrants) at grid cell i between period t and t + 1, Ñou i,t is the latent number of orangutans at grid cell i in survey period t, as a result of individuals survived and recruited in the previous survey period (S i,t−1 and R i,t−1 , respecitively), Zou i,m,t is the observed orangutan occurrence at grid cell i in survey period t from respondent m Nnest i,t is the latent number of orangutan nests at grid cell i in survey period t, Onest i,t is the latent occupancy of orangutan nests at grid cell i in survey period t, derived as a binary value of Nnest i,t Y i,j,t is the observed nest count at grid cell i in survey period t from survey type j, Znest i,k,t is the observed nest occurrence at sub-grid cell k and grid cell i in survey period t.
The parameters estimated from the model are the initial abundance rate at grid cell i (λ i ), survival probability and recruitment rate at grid cell i between survey period t and t + 1 (θ i,t and δ i,t ), the orangutan occupancy rate at grid cell i and survey period t (ϕ i,t ), the scaling factor of the nest and the orangutan density at grid cell i and survey period t (ψ i,t ), the probability of detecting orangutan individuals from the interview survey at grid cell i and survey period t for respondent m (ρou i,m,t ), the probability of detecting orangutan nests from the line transects at grid cell i and survey period t for survey type j (ξ i,j,t , where j ∈ {aerial, ground}), and the probability of detecting orangutan nests from the line transects and other targeted surveys at sub-grid cell k and grid cell i and survey period t (ρnest i,k,t ).
These parameters can be modeled by including site-specific covariates. We modeled the initial abundance rate at grid cell i, i.e. λ i , as a function of altitude (ALT i ), mean annual monthly rainfall during the dry season from May to September (DRY i ), mean annual monthly rainfall during the dry season from November to March (WET i ), the quadratic term of ALT i , DRY i and WET i , nearest distance to protected areas (DPA i,1 ), the proportions of Muslims per district (MS i ), natural forest extent (FR i,1 ), and the interaction between natural forest extent and nearest distance to forest recently converted to industrial agriculture (FR i,1 × CFA i,1 ) that all occurred prior to 2003, i.e.  We included the quadratic term of ALT, DRY and WET to test the preference of orangutan to occupy areas with intermediate values for altitude and rainfall during the dry and wet season. We also tested whether or not proximity to protected areas (DPA) increases survival rates by reducing the risk of orangutan killings. Descriptions of the covariates used to explain the initial abundance, occupancy and survival rates are given in Supplementary Method 1.
The recruitment rate at grid cell i between period t−1 and t, i.e. δ i,t , was modeled as the number of individuals in site i and the neighboring sites at the previous survey period 62  . This is to simulate the possible dispersal routes taken by an orangutan from grid cell i to the surrounding grid cells. We then intersected this line with the river barrier layer. We assumed w k = 0 if at least one intersection was found within grid cell k ∈ n j (i.e. rivers prevent orangutan dispersal from grid cell i to grid cell k) and w k = 1 if no intersection was found.
In earlier studies, the density of orangutans at grid cell i, i.e. gou i , has typically been estimated by the following equation where b i is the proportion of nest builders, i.e. juveniles less than around 3 years of age are unlikely to build nests 64 , q i is the daily rate of nest production, and d i is the nest decay rate or the number of days a nest remains visible. Based on previous studies in Borneo, the proportion of nest builders has been estimated at around 0.9 4, 7, 23 .
The average daily rate of nest production for Bornean orangutans has been estimated to range between 1 and 1.2 4,7,23 , but this can fluctuate depending on the level of forest disturbance, i.e. between primary and logged over forest 23 . Generally, the multiplication of b i and q i results in a value around 1. The nest decay rate is much more uncertain, however, ranging between 85 to over 800 days [21][22][23] and has been shown to vary across different forest types and with altitude 4,7,23 . Hence, to take into account the variability in the total denominator of Eq. (5) across different grid cells i and survey periods t, we modeled ψ i,t as where MGV i,t is a binary variable denoting whether or not the majority of forest at grid cell i and time t are mangrove forest, and similarly PT i,t for peat forest, LOWL i,t for lowland forest (altitude < 500 m), MONT i,t for montane forest (altitude ≥ 500 m), and FRGM i,t for highly fragmented forest (< 25 ha per km 2 ). The probability of detecting orangutans from the interview surveys at grid cell i and time t for respondent m, i.e. ρou i,m,t , was modeled as a function of respondents' frequency for entering the forest around the village (1 for more than once a month and 0 for less than once a month), i.e. FE i,m , such that i m t i m , , The probability of detecting orangutan nests at grid cell i and time t and for survey j (j ∈ {aerial, ground}), i.e. ξ i,j,t , was modeled constant for each survey type, such that Finally, the probability of detecting orangutan nests at sub-grid cell k and grid cell i and time t for line transects and other targeted surveys, i.e. ρnest i,k,t , was modeled constant, such that Model fitting and evaluation. We used WinBUGS Version 1.4.3 65 to estimate the parameter posterior distributions and the regression coefficients for λ i , ϕ i,t , θ i,t , δ i,t , ψ i,t , ρou i,m,t , ξ i,j,t , and ρnest i,k,t . The WinBUGS code for the dynamic abundance model is provided in Supplementary Method 2. We assumed a vague prior for each parameter, as described in Table 2.
We ran three Markov chain Monte Carlo (MCMC) chains, where each chain consists of 100,000 iterations and the first 50,000 were discarded as burn-in. To improve convergence and to reduce the autocorrelation in the MCMC chain, we standardized all variables prior to model fitting. Prior to fitting the model to the data, we tested the correlation among the original (unstandardized) environmental variables explaining λ, ϕ t and θ t , i.e. variables ALT, DRY, WET, DPA, MS, FR and CFA, and also among the standardized variables. Convergence for each model parameter was assessed from the values of Rhat statistics and visualization of the chain plot of the MCMC iterations. Rhat values around 1 and the absence of seasonality within each chain plot and overlap among the chains indicate convergence. We also tested for correlations among posterior distributions of the coefficients, especially between the linear and the quadratic terms of variables ALT, DRY and WET, to ensure correct functional forms were specified for these variables and the coefficients were not biased.
The goodness-of-fit of the model was assessed by comparing the simulated nest abundance predictions for each time period with the observed nest counts. For each simulated prediction and time period, we calculated the Pearson's correlation coefficient r and also fitted a linear regression between the predicted values and the observed values to calculate the R 2 value 66 . We also validated the simulated orangutan presence-absence predictions for each time period against the actual observations based on interview surveys. In the validation dataset, we defined "presence" in a village if at least one respondent reported the occurrence of orangutan, and we defined "absence" if more than 50% of the respondents who enter the forest more than once a month had never seen the species. We used the proportions of correctly predicted presence or Sensitivity (SN) and the proportions of correctly predicted absence or Specifity (SP) as the measure of performance. SN and SP values close to one indicate high accuracy.
Assessing orangutan abundance change among regions and land uses. We assessed orangutan population trends by measuring the change in the number of individuals obtained from the simulated predictions. We investigated how the trends vary across different regions (states and provinces), as well as across different land uses. We considered five land use categories: (1)  Assessing drivers of orangutan population decline among regions and land uses. To inform orangutan conservation planning, we assessed the drivers of orangutan population decline in each region. This was achieved mainly by relating the environmental covariates explaining survival rates in Eq. (3) across 1 × 1 km 2 grid cells where orangutans are predicted to occur with known actual threats observed on Borneo. These threats includes: (1) habitat loss, i.e. the loss of natural forest of orangutan habitats, (2) human-orangutan conflicts, (3) anthropogenic human activities, such as hunting and poaching, and (4) habitat fragmentation, i.e. breaking up intact forest habitats into small forest patches.
The decline of orangutan population due to habitat loss in grid cell i at time period t, i.e. HLOSS i,t , was related specifically to forest cover covariate FR i,t (i.e. the 10 th additive component in Eq. (3)). We measured habitat loss based on counterfactual analysis, i.e. the discrepancy between the survival rates under the 'counterfactual assumption of no forest loss, or forest cover remains the same as in the previous time period (FR i,t−1 )' versus 'the actual forest cover in that period (FR i,t )' , such that High HLOSS i,t implies low orangutan survival rate, or high contribution of habitat loss to population decline in grid cell i at time period t. The decline of orangutan population due to human-orangutan conflicts in grid cell i at time period t, i.e. CONFL i,t , was related specifically to the interaction between forest cover FR i,t and the distance to newly converted forest to industrial agriculture CFA i,t (i.e. the 11 th additive component in Eq. (3)), such that , , Low CONFL i,t implies low orangutan survival rate, or high contribution of human-orangutan conflicts to population decline in grid cell i at time period t. For measuring the decline of orangutan population due to anthropogenic activities in grid cell i at time period t, i.e. ANTH i,t , we used monthly rainfall during the dry DRY i and the wet seasons WET i and proximity to protected areas DPA i,t as proxy (i.e. 4-8 th additive components in Eq. (3) where η η η η ηˆˆˆ, , , , and 4 5 6 7 8 are the estimated coefficients obtained from WinBUGS simulations. This is because seasonal rainfall patterns determine socio-economic structure and livelihoods on Borneo 34 . Additionally, protected areas were assumed to provide a refuge for the species against hunting and poaching 11 . Low ANTH i,t implies low orangutan survival rate, or high contribution of anthropogenic activities to population decline in grid cell i at time period t.
To obtain the relative influence of habitat loss as a driver of orangutan population decline for each region, we averaged HLOSS i,t across all grid cells where orangutans are predicted to occur within the respective region. To obtain the relative influence of human-orangutan conflicts and anthropogenic activities as drivers of population decline for each region, we applied similar procedure to CONFL i,t and ANTH i,t , respectively. We also assessed how these drivers vary across different land uses (i.e. LOGG, ITP, OPP and OTHER) within each region.
For habitat fragmentation, we assessed this as a driver over the entire orangutan distribution range across different landscapes within the region. Because territorial ranges of orangutans, especially the females, are generally restricted to a maximum of 850 ha 55 , the species' dispersal opportunities between habitat fragments are generally limited. This implies that landscapes with isolated forest patches of orangutan habitats (i.e. fragmented habitats) have a higher risk of orangutan decline due to lower colonization rates than landscapes with better habitat connectivity. The relative influence of habitat fragmentation as a driver of orangutan population decline in a region, i.e. FRAG, was estimated as the interaction between the mean size of contiguous forest where orangutan occurred and the mean reciprocal distance of each contiguous forest to the nearest forest patch. Low FRAG implies low orangutan survival rate, or high contribution of habitat fragmentation to population decline.