Study area

This study was conducted in the RESEX-TA, a terrestrial extractive sustainable-use protected area in western Pará (Brazilian Amazonia; 02º20′–03º40′ S, 55º00′–56º00′ W; Fig. 1) created in 1998 to protect local people from the impacts of the expansion of logging companies (Oliveira et al. Reference Oliveira, Junior and Chaves2005). Covering a total area of nearly 650 000 ha, the RESEX-TA is currently the most human-populated extractive reserve in Brazil, housing c. 23 000 residents (Silva et al. Reference Silva, Nascimento, Magalhães and Spínola2022). It comprises the first Brazilian protected area of this category to develop a management plan, a technical document that establishes the zoning and norms that guide the use of the protected area.

Figure 1. Location of the Tapajós–Arapiuns Extractive Reserve (RESEX) in the cities of Santarém and Aveiro (Pará), Brazil.

The predominant vegetation is the ombrophilous forest (66.19% of the territory), characterized by large trees, with abundant woody, epiphytic and liana plants (Veloso et al. Reference Veloso, Rangel-Filho and Lima1991, Instituto Brasileiro de Geografia e Estatística 1992) in terra firme (unflooded) forests. Other common plant typologies encompass secondary forest (8.88%), savannah (0.62%) and flooded forest (Carvalho-Junior Reference Carvalho Junior2008, Montag et al. Reference Montag, Prudente, Ferreira, Dutra, Benone, Barbosa and Ruffeil2012). The climate is humid tropical (Ami), with a rainy season from November to May and a dry season from August to October (Instituto Nacional de Pesquisas Espaciais & Centro de Previsão de Tempo e Estudos Climáticos 2023).

The RESEX-TA has 78 communities inhabiting areas along the Tapajós and Arapiuns rivers (Silva et al. Reference Silva, Nascimento, Magalhães and Spínola2022), in addition to a few families living along small streams (Oliveira et al. Reference Oliveira, Junior and Chaves2005), Instituto Chico Mendes de Conservação da Biodiversidade 2014; Fig. 1). All communities are established in terra firme forests, in which the residents directly depend on the exploitation of natural resources (hunting, fishing and extraction of timber and non-timber forest products) and/or family farming and small-animal breeding (Saúde & Alegria 2012). Hunting is allowed for subsistence purposes, and its management depends on internal agreements signed between communities and reserve managers. In mid-2021, there was an attempt at logging through a community forest management plan, but the lack of prior popular consultation resulted in the cessation of activities.

Participatory Biodiversity Monitoring Project

In 2014, the Brazilian federal government’s Instituto Chico Mendes de Conservação da Biodiversidade (ICMBio) implemented the Brazilian Biodiversity Monitoring Program (Monitora Program) in the RESEX-TA. Given that local residents expressed their willingness to understand the current and future population status of terrestrial game fauna (i.e., medium- to large-sized mammals and birds), the monitoring programme in the RESEX-TA focused on species belonging to these groups, as well as reptiles (total of 25 target species, carefully chosen by local people based on their hunting activities; Table S1). This programme was carried out via the Participatory Biodiversity Monitoring Project (PMPB), a partnership between ICMBio and the Instituto de Pesquisas Ecológicas (IPÊ; Ribeiro Reference Ribeiro2018).

Firstly, various meetings were conducted in the RESEX-TA to identify potential monitors belonging to these communities who would become responsible for transect establishment and data collection. In particular, transects needed to be placed in areas accessible on foot or by boat/canoe to the monitors, in areas lacking direct human impacts apart from hunting activities. Furthermore, transects were placed at various distances from the nearest community, with transects nearby being more likely to be affected by hunting pressure (Sampaio et al. Reference Sampaio, Morato, Royle, Abrahams, Peres and Chiarello2023). Prior to effective implementation, potential areas were identified through maps, taking into account elevation profiles, obtaining coordinates from maps and on-site recognition of the areas. Finally, transects needed to be at least 5 km from each other and placed in non-flooded forests, which represent the main land cover type in the RESEX-TA; a total of nine linear transects of 5 km in length were thus established. However, after beginning data collection in 2014, one community left the programme, leaving only eight communities/transects in the programme (Fig. 1). Each target community selected three to four monitors older than 18 years of age based on their willingness to participate combined with vast experience in walking along the forest. Subsequent training courses were provided by ICMBio and partners to all monitors, focused on providing instructions for game fauna surveys. Since 2022, data collection has been performed by 31 previously trained monitors.

Data collection

On each transect, pairs of monitors conducted linear transect surveys of the transects in their neighbourhood. The methodology widely used in tropical forests (Peres Reference Peres1999, Benchimol Reference Benchimol and TH2016) consisted of walks at a slow and constant speed (average 1 km/h) along a linear transect, searching for both visual and acoustic records of the target species. Surveys were conducted only in the morning (06.30–11.30) by a pair of monitors and did not occur on rainy days. Whenever an individual or group of target species was observed, the following information was recorded: species name, type of record, time of encounter, location of the animal on the transect, number of individuals and perpendicular distance to the first detected individual, subgroup or group (measured with a tape measure). Given that local monitors had extensive experience of walking in the forest and of searching for and identifying animals, and given that appropriate training courses were also provided, the data collected were consistent and robust. Surveys were conducted between 2015 and 2020, with an average of nine (SD = 2.59) survey days per transect per year, totalling 384 sampling days across the eight transects. Although initially the protocol aimed to obtain data from the dry and wet periods for 10 days per year per transect, data were collected in different months depending on the availability of the monitors, resulting in the sampling effort varying amongst transects. Furthermore, there were situations in which the linear transects extended throughout the day or rainy periods that only enabled partial sampling; such data were excluded from the analysis. The cumulative distance sampled in each transect varied from 5 to 65 km per year, with the total distance being 1915 km over the 6 years (mean ± SD = 239 ± 45 km per transect; Table S2).

Target taxa

We used data collected from nine game forest taxa (i.e., nine taxonomic orders or families of hunted medium- and large-sized forest mammals and birds) in the RESEX-TA (Reis et al. Reference Reis, Valsecchi and Queiroz2018, Reference Reis, Soares, Maduro, Spinola and da Rocha2022; Table S1). The selected focal taxa were: those of the Monitora Program monitoring protocol, which considered the target hunting species in the region based on local people’s knowledge; and those recorded in hunting events (not covered in this study) from the Monitora Program (Reis et al. Reference Reis, Valsecchi and Queiroz2018, Reference Reis, Soares, Maduro, Spinola and da Rocha2022). Specifically, for our study, we selected game species from those two categories that could be recorded within the RESEX-TA through the use of the linear transect technique (Table S1). A species-level analysis was not feasible due to the limited number of records for most animals. There were six taxa of medium-sized animals (0.54–5.00 kg for birds, 0.92–5.00 kg for mammals): Tinamidae (Tinamus spp., Crypturellus spp.), Psophiidae (Psophia viridis), Cracidae (Penelope spp., Crax sp. and Pauxi tuberosa), Dasyproctidae (Dasyprocta croconota), primates (Cebus unicolor, Sapajus apella, Pithecia irrorata, Callicebus hoffmannsi) and Procyonidae (Nasua nasua); and three taxa of large animals (>5 kg): Atelidae (Alouatta nigerrima), Cervidae (Mazama americana and Mazama nemorivaga) and Tayassuidae (Dicotyles tajacu and Tayassu pecari). Despite some species not frequently being hunted at other Amazonian sites (e.g., P. irrorata, C. hoffmannsi and S. apella), they were being captured for consumption in the RESEX-TA (YMS Reis, unpublished data 2023). In addition, we did not include Ateles marginatus as a focal species as it has been extirpated from the RESEX-TA by previous severe hunting (Peres et al. Reference Peres, Barlow and Haugaasen2003).

Density and biomass

We estimated the density of each of the nine animal taxa in each transect and year using Equation 1:

(1) $$D = {N \over {2\ \times\ EW\ \times\ L}}$$

where D represents the density (individuals/km²), N represents the total number of sightings, EW represents the effective width of the transect (in km) and L represents the total distance travelled (Bernardo & Galetti Reference Bernardo and Galetti2004). To obtain EW, we grouped the visualization events per taxa and constructed six competing models using the Distance 4.1 software: a uniform key function with either cosine or simple polynomial series expansion; a half-normal key function with either cosine or hermite polynomial series expansion; and a hazard rate key function with cosine or simple polynomial series expansion (Thomas et al. Reference Thomas, Buckland, Rexstad, Laake, Strindberg and Hedley2010; models in Table S3). We selected the best-fit model based on the lowest Akaike information criterion (AIC) values (Burnham & Anderson Reference Burnham, Anderson, Burnham and Anderson2002) and coefficient of variation and the highest goodness-of-fit χ² test probability values. From the estimated density per taxa, we evaluated the effect of hunting on animal densities considering (1) group (medium- and large-sized) density and biomass and (2) ‘individual’ densities of the four most frequent taxa (Tinamidae, Dasyproctidae, Primates and Cervidae) for each transect in each year. Given the low number of records, we were unable to conduct individual analyses for the other taxa. Biomass was calculated as the sum of densities of each species multiplied by the mean body mass of the species in each group (based on Robinson & Redford Reference Robinson and Redford1986, Ayres & Ayres Reference Ayres and Ayres1979, Valsechi Reference Valsechi2013).

Anthropogenic stressors

We obtained two metrics related to anthropogenic pressure. Firstly, we used a proxy of HI based on the distance to nearest community and human population size of its community (adapted from Scabin & Peres Reference Scabin and Peres2021), according to Equation 2:

(2) $$HI = \sum {{S\left( {com} \right)} \over {\sqrt {d\left( {com} \right)} }}$$

where S represents the human population size in each community (com; Table S4) and d represents the Euclidean distance from the centre of each transect to the centre of each nearby community (Table S4). For this, we considered all communities located up to 10 km from the centre of the surveyed transect. This distance was based on the evidence that hunters forage in areas located up to 10 km from their homes in the RESEX-TA (YMS Reis, unpublished data 2018). Distances were calculated in QGIS 2.18.9 (QGIS Development Team 2017) and the human population size of each community was obtained from the local ICMBio database from 2022. We assume that population size did not substantially change amongst surveyed years based on a comparison between the most current list of protected area beneficiaries (from 2022) and 2018, the previous year for which we had data (Spearman correlation, r = 0.95, p = 0.001).

We also measured the distance to the nearest road, defined as the distance from the middle of each transect to the nearest dirt road (Table S4). Because dirt roads are the main access ways in terra firme in the RESEX-TA, we assumed that transects closer to roads would be more accessed by humans than isolated transects (Benítez-López et al. Reference Benítez-López, Alkemade, Schipper, Ingram, Verweij, Eikelboom and Huijbregts2017, Beirne et al. Reference Beirne, Meier, Mbele, Menie, Froese, Okouyi and Poulsen2019).

Data analyses

We used generalized linear mixed models to relate variation in HI, distance to the nearest road and year of data collection (with the random term ‘transect’) to overall densities and biomasses of both medium-sized and large-sized groups separately and to the densities of the four selected taxa individually. We emphasize that it was only possible to analyse diurnal species that were commonly recorded (i.e., those that had greater detectability and therefore provided sufficient visual records for analysis); we excluded those game species that are rarer or more difficult to record using linear transects. Spearman tests between our predictor variables (i.e., HI, distance to the nearest road and year) showed that none was correlated at (>|0.50|) with any other, and thus we maintained all of these variables in the global model (mean, standard deviation, minimum and maximum values of anthropogenic variables in Table S5).

For each response variable, we constructed different models: the null model, models presenting a single or two predictors, models with all possible additive combinations of them and models with the interaction effects between HI and distance to the nearest road. We selected the best model(s) based on the AIC values (Burnham & Anderson Reference Burnham, Anderson, Burnham and Anderson2002), considering as parsimonious those presenting ΔAICc ≤ 2.0. When the null model appeared amongst the parsimonious models, we considered that no other model best explained the specific pattern than chance. We used the 95% confidence interval (CI) of the coefficient of each variable included in the best models and considered its influence when the 95% CI did not include 0. Due to convergence problems, we were unable to consider the interaction model for primates and Tinamidae. We transformed the data using log10 for groups or square root for taxa, and we confirmed data normality using the Shapiro–Wilk test. We used the Gaussian distribution family (‘identity’ linkage function) in all models. The analyses were performed using the package ‘glmmTMB’ (Brooks et al. Reference Brooks, Bolker, Kristensen and Maechler2022) in R 4.2.1 (R Core Development Team 2022).

To assess population trends over the years, we built state-space models that describe the stochastic and deterministic relationships between the observable and unobserved values using a Bayesian approach (Royle & Kéry Reference Royle and Kéry2007, Kéry Reference Kéry2010). We used as a response variable the annual mean density values of the medium-sized and large-sized groups and nine taxa considering all transects. The data were log transformed. The assessment was made by estimating the exponential growth rate r of density, in which the r value is a measure of change in population size that assumes positive, negative and zero values in increasing, declining and stable populations, respectively, and is calculated using Equation 3:

(3) $$r = {{log\left( {{N_t}/{N_0}} \right)} \over t}$$

where N ₀ is the population size at the beginning of the period and N _t is the population size after t years (Caughley & Sinclair Reference Caughley and Sinclair1994). A time-series mean r value ≥ 0 could suggest that populations are stable or growing, but values r values < 0 could be interpreted as evidence of a decline in the evaluated parameter (de Paula et al. Reference de Paula, Carvalho, Lopes, de Alencar Sousa, Maciel and Wariss2022). Thus, the r value was considered as indicative of the behaviour of density over the years. If the r values over the time series were equal to 0 or positive for the medium-sized and large-sized groups or the nine ‘individual‘ taxa, we assumed that populations were not declining; declining r values were taken to indicate otherwise. We used R version 4.2.1 (R Core Development Team 2022) and the ‘R2jags’ package (Su & Yajima Reference Su and Yajima2012) to fit the state-space model in JAGS (Plummer Reference Plummer2015). Non-informative priors were used for the density, and 200 000 Markov chain Monte Carlo iterations were run in two independent chains with a burn-in of 100 000 and a thinning factor of 0.06. Furthermore, the Gelman–Rubin (Rhat) diagnosis was used to evaluate parameter convergence (Gelman & Shirley Reference Gelman, Shirley, Brooks, Gelman, Jones and Meng2011), where values of 1.001 indicated a satisfactory model (the closer to 1, the better the model) and 1.1 was considered the acceptable limit. If the Bayesian 95% CI did not include 0, the r value was deemed substantially different from 0. The results are expressed as means with standard deviations.