Spatial aggregation of fruits explains food selection in a neotropical primate (Alouatta pigra)

The availability and spatial distribution of food resources affect animal behavior and survival. Black howler monkeys (Alouatta pigra) have a foraging strategy to balance their nutrient intake that involves mixing their consumption of leaves and fruits. The spatial aggregation of food items should impact this strategy, but how it does so is largely unknown. We quantified how leaf and fruit intake combined (here termed food set selection) was spatially aggregated in patches and how food aggregation varied across seasons. Using variograms we estimated patch diameter and with Generalized Least Square models determined the effect of food spatial aggregation on food selection. Only fruits were structured in patches in the season of highest availability (dry-season). The patches of food set selection had a diameter between 6.9 and 14 m and were explained by those of mature fruit availability which were between 18 and 19 m in diameter. Our results suggest that the spatial pattern of food selection is influenced by patches of large fruit-bearing trees, not by particular species. Fruit also occur along spatial gradients, but these do not explain food selection, suggesting that howlers maximize food intake in response to local aggregation of fruit that are limiting during certain seasons. We demonstrate how the independent spatial modelling of resources and behavior enables the definition of patches and testing their spatial relationship.

The availability and spatial distribution of food resources affect animal behavior and survival. Black howler monkeys (Alouatta pigra) have a foraging strategy to balance their nutrient intake that involves mixing their consumption of leaves and fruits. the spatial aggregation of food items should impact this strategy, but how it does so is largely unknown. We quantified how leaf and fruit intake combined (here termed food set selection) was spatially aggregated in patches and how food aggregation varied across seasons. Using variograms we estimated patch diameter and with Generalized Least Square models determined the effect of food spatial aggregation on food selection. Only fruits were structured in patches in the season of highest availability (dry-season). the patches of food set selection had a diameter between 6.9 and 14 m and were explained by those of mature fruit availability which were between 18 and 19 m in diameter. Our results suggest that the spatial pattern of food selection is influenced by patches of large fruit-bearing trees, not by particular species. Fruit also occur along spatial gradients, but these do not explain food selection, suggesting that howlers maximize food intake in response to local aggregation of fruit that are limiting during certain seasons. We demonstrate how the independent spatial modelling of resources and behavior enables the definition of patches and testing their spatial relationship.
Food resources generally have an aggregated spatial distribution 1 . Therefore, the spatial distribution of diversity, abundance, and behavior of animals that consume such resources are neither uniform nor random 1 . Typically, spatial modeling of animal behavior and feeding patterns have been carried out at a landscape scale 2,3 . At this scale, the effects of spatial variation (e.g., fragment size, forest cover) influence habitat selection and animal behavior 4,5 . However, smaller scales variation can also have an effect. For example, Razgour et al. 6 documented that habitat selection by a bat (Plecotus austriacus) has different drivers on broad scales (limited by inadequate climatic conditions) than on small scales (limited by the availability of preferred foraging areas). Jedlikowski et al. 7 studied habitat selection of two rallid birds at three spatial scales, landscape (50-300 m), territory (14 m), and nest (3 m), and found that bird occurrence was most associated to the territory scale as they make foraging decisions at that scale. In socioecological studies of agonism in primates, it has been suggested that the feeding area that an individual can defend is dependent on body size at small scales 8 . Small species typically defend small foraging areas 8 .
Primates are ideal to examine the effect of spatial scale on behavior and distribution as observations can be carried out on each individual in a group and can include the behavioral variability among groups and populations. The time spent resting and eating leaves by spider monkeys (Ateles geoffroyi) were mostly correlated with the forest cover at the 126-ha scale; these correlations were lower when broader scales were considered 3 . The tree species. This is calculated by adding the density, relative abundance, and dominance of each tree species (sp) for each fragment (modified from Agostini et al. 38 ; Salomao et al. 39  where N sp is the number of trees of sp species, and S f is the area of fragment f; N t is the total number of trees in the fragment. BA is the basal area of each i species (ranging from 1 to n marked trees of the same species), or of tf, all trees in the fragment added.
Activity and feeding of black howler monkeys. We followed two groups of howler monkeys from June 2016 to May 2017 completing 106 days of ≈12 hr of daily observations. Group-I (G-I: 426.9 hours) had 1 adult male, 2 adult females, 2 infants. Group-II (G-II: 433.1 hours) had 1 adult male, 2 adult females, 1 infant. We followed six adult focal animals, each during two days per individual per month for each group. We recorded the time invested: feeding, resting, traveling, and other (i.e., playing, vocalization, and agonistic encounters). We recorded the time devoted to the consumption of each of the five food categories that constitute >90% of howlers' diet 40 : mature leaves, young leaves, mature fruits, immature fruits, and flowers. We estimated the dry weight intake in grams of the items 41 . For this, we calculated a feeding rate for each item by counting the number of food units eaten per minute. To calculate dry weight intake, we collected food items from the tree crowns using 15-m poles and weighed 100 units of each item to estimate the fresh weight mass consumed within a maximum of five days from the first observed consumption of a new food. Subsequently, we collected between 150 and 250 g of fresh material per food item, and dried the samples in a cool, dark room at room temperature (ranges between 21.7 °C and 33.5 °C). To avoid the fermentation and molding of fruits, we used a food dehydrator (NESCO-FD-2000) at 55 °C until the sample reached a constant weight.
food availability and selection. We calculated two indices of food availability.
Intraspecific index of Food Availability (IFA) measures the availability of food items weighted by the size of the tree: where the IFA value of an I food item of the tree a is equal to the average of the PC I phenological category of the item in a given season for the tree a, multiplied by its DBH a . Interspecific index of Food Availability (IDA) measures the availability of food items per tree given the importance value of the species (equivalent to FAI index used in Agostini et al. 38 and Chaves and Bicca-Marques 28 ):  where the IDA value of an I food item is equal to the average of the PC I phenological category of the item in a given season for the tree a, multiplied by the importance value index (IVI) of the species (sp) to which that tree belongs 39 . Note that the estimation of dominant species through the IVI index uses the absolute tree density. In our study this is constant because all tree species in a plot share the same area. However, we have left this component to make IDA applicable to situations where sites with different basal areas and tree species composition are being compared in terms of animal choices. While IFA gives greater weight to the contribution of large trees in the availability estimate, IDA magnifies the contribution of trees if they belong to one of the dominant species in the fragment. Using both indexes allowed us to consider the effect on food item availability by individual tree variation from the effect due to the importance of some species in the fragment. Finally, we evaluated the food selection of each food item type separately and by food set (leaf and fruit intake combined, i.e., the sum of all items consumed) for of each tree by adding feeding bouts: where SF (selected food) intake of the item I in a given season is the sum of the dry weight intake in all feeding bouts FB (from i=1 to n bouts) of this item from the tree a.
Data analyses. Spatial structuring of food availability and selection. The scale of the study is defined by grain, extension, and range, but only extension differed between the fragments 16,42 . Extension was 0.8 ha for G-I site, 1.0 ha for G-IIA site and 1.1 ha for G-IIB site (calculated using geometry option of QGIS 43 ). For all three sites, the grain was 3 m corresponding to the average diameter of tree crowns. The interval was 0.5 m, the minimum distance between trees. We evaluated the spatial structuring of food availability modeling the following variable combinations (96 models): 2 availability indexes (IFA, IDA) x 4 food items (YL, ML, MF, IF) x 4 seasons (Rainy, Nortes, Dry, Annual) x 3 Sites (G-I, G-IIA, G-IIB). The spatial structuring of selected food (SF) was evaluated using the following combinations (60 models): 5 food items (YL, ML, MF, IF, Food-Set) x 4 seasons (Rainy, Nortes Dry, Annual) x 3 Sites (G-I, G-IIA, G-IIB). The spatial modeling method (and its R code) used in this study is described in detail in Negrete-Yankelevich and Fox (2015) 16 ; however, we briefly summarize it here. Spatial modelling was done in two steps, modeling of gradients followed by patches 16 . We log or root-transformed the variables to comply with the assumption of normality of the residuals of linear modeling (see Table 1 for variable-specific transformations).
To detect possible gradients (i.e., monotonic changes in a given direction), we fitted linear models using north-south (y) and east-west (x) coordinates of each tree as explanatory variables. Note that these models are descriptive (i.e., not hypothesis tests), therefore the spatial independence in the data is not required 16 . To consider different forms of gradients (monotonic changes to a given direction or combinations of directions) we started the modeling of variables including the combination of different polynomial forms of their position in space as explanatory variables (x + y + x 2 + x × y + y 2 ). We simplified the models by successively eliminating the components of x and y from the polynomial equation. Variables were removed in increasing order of the estimated slope. The change Akaike information criterion (ΔAIC = AIC null -the AIC model ) was calculated for the model with and without the variable being examined. AIC null was the AIC of the model where the response variable is solely explained by its mean. In each step, the model with the highest ΔAIC was kept before proceeding to the elimination of the next variable. We considered the existence of a gradient if at the end of the simplification the model included an x or y explanatory variable and had a ΔAIC greater than five units 44 .
To model the patch structure of the variables, we used variography. This is a method used to estimate the average size of patches present in the residuals of the linear models in the first step, once the variance explained by the gradient has been eliminated. Experimental variograms were constructed, which are graphs of semi-variance (an inverse measure of the autocorrelation) as a function of separation distance between pairs of observations 16,17 . We calculated: where γ (h) is the semi-variance of all pairs of observations located at h distance (lag), W(h) is the total number of observation pairs separated by h, and y i is the value of the sample at each location, tree position in this case. We followed the recommendation to construct the variograms using pairs of observations that were <2/3 of the observed maximum distance (60 m) between observed points 16 . We did this to avoid including distances represented by a small number of pairs which could bias the estimates of the theoretical models that were subsequently fitted 16 . Variograms that increase and reach an asymptote indicate the existence of patches. The distance at which this asymptote is reached represents the distance at which the autocorrelation ceases and the average diameter of the patch in which the variable is aggregated. To assess if there was a significant patchy structure and to estimate the average patch diameter, we fitted the experimental variograms with theoretical variograms using weighted least squares 16,17 . This strategy is recommended for irregular spatial distributions of sampling points due to the variability that results in the number of pairs available per lag 16,45 . The theoretical models fitted to the variograms were: In the exponential and spherical models, the sill is the semi-variance reached when the variogram reaches the asymptote denoted by (C 0 + C 1 ). C 0 represents the nugget effect or ordinate at the origin and C 1 is the increment in semi-variance from C 0 up to the point at which the variogram reaches the asymptote. For spherical models a is the interval, and for the exponential models the interval is taken when the semi-variance reaches 95% of the asymptote 17 . We only considered those models that explained at least 30% of the variance (R 2 ≥ 0.3) as plausible and chose the theoretical model with greater explanatory power (R 2 ). We took this proportion of explained variance as a proxy of the degree of space structuring for each variable 16 .
Linear modeling of food availability and food set selection. We analyzed the explanatory power of mature fruits and young leaves availability on food set selection (per group) using linear models. We only modeled those seasons when food set was structured in space. Position in space (x + y) of each observation was included in the initial models as explanatory variables to test for the presence of gradients not explained by either availability index. We simplified these models comparing the ΔAIC of models including both explanatory variables together and separately, and chose the model with the lowest ΔAIC, which had a reduction of at least five units with respect to the null model 44 . During this simplification, we built experimental variograms with the residuals of each linear model and fitted theoretical models in the same way we did when we analyzed for spatial structures. We did this to determine (1) if they fulfilled the assumptions of spatial independence of linear models by not presenting spatial autocorrelation (nugget model as best fit) and (2) if the spatial structures (gradients and patches) found in SF disappeared in the residuals when using either of the availability indexes as explanatory variables. If the residuals did not show spatial patterns when the model included the variable, but they did when it excluded it, we considered that the variable explained both the mean variation (and slope) and the spatial variation of SF. If the residuals had spatial structures, we considered that the variable only explained the mean variation of SF 16 . If we still detected a spatial signal in the residuals of the adequate minimum linear models, we proceeded to construct generalized least squares (GLS) models. GLS are linear models that directly include in the residuals a spatial covariance structure through a variance-covariance matrix that considers the range and sill of a variogram function fitted to the residuals 16 . When considering the spatial structuring of the residuals, GLSs allow an adequate estimation of the parameters associated to the explanatory variables, without overestimating degrees of freedom, while estimating the size of the patches not explained by the model.
We additionally constructed linear models of the selected food explained by food availability for each item individually, but the saturated models were not simplified. Only the residuals where examined to confirm that there was no spatial structure being masked by a strong influence of explanatory variables 16 . All the analyses were carried out in R, v. 2.15.3 (http://R-project.org/). We used the geoR package for variography 46 , and the nlme package for linear models and GLS 47 . ethical approval. All applicable international, national, and/or institutional guidelines for the care and use of animals were followed. All procedures performed in studies involving animals were in accordance with the ethical standards of the institution or practice at which the studies were conducted. This research was approved by SEMARNAT (SGPA/DGVS/10426/14), the governmental institution that regulates animal research, and complied with the laws of the country of Mexico.  Table S1). Leaves were more available than fruits throughout the year at both sites ( Fig. 2; see median, min and max values of IFA and IDA in Supplementary Table S2). Mature leaves were constantly available, while young leaves and fruit exhibited seasonal variation, with fruit exhibited a peak in availability during the dry season.

Results
Spatial structuring of food availability and selection. Foods were aggregated in patches nested in gradients in seven combinations out of the 21 with proven spatial structuring, suggesting two scales of spatial structuring for a few availability combinations. Two combinations showed gradients without patches (only in site G-IIB), and 12 combinations had patches without gradients. The degree of spatial structuring, denoted by R 2 of the theoretical model of the variogram, depended on the food item and season (parameters estimated in Table 1 and variograms in Supplementary Fig. S3, parameters of unstructured variables in Supplementary Table S4). In seven cases out of nine combinations with gradients, food items were along a north-south direction (y) and the remainder with east-west gradients (x). The intraspecific availability index (IFA) in G-I site showed east-west gradients (x) in mature leaves in the "nortes" season and young leaves in the rainy season, and x and y gradient (2019) 9:19452 | https://doi.org/10.1038/s41598-019-55932-y www.nature.com/scientificreports www.nature.com/scientificreports/ in mature fruit in the dry season (see the bubble maps representing a gradient in Fig. 3). Only a north-south gradient was found for the interspecific availability index (IDA) in immature fruits in the dry season. For the G-II site, the IFA showed a north-south gradient for mature fruits in the dry season (Fig. 3), young leaves in "nortes" season, and immature fruits for the whole year. The IDA showed this same gradient in mature fruits in the dry season and for the whole year.
The variographic analysis illustrated that food availability was frequently aggregated in patches (19 combinations, 16 for fruits). Neither mature nor immature leaves showed a consistent aggregation pattern between sites, and leaves were aggregated in patches only at site G-I in two seasons. Availability indices of young leaves were not aggregated during the season of higher availability, "nortes" season, and were only aggregated in the rainy season (mean patch diameter of 15 m, R 2 = 0.5). According to IFA, availability of mature leaves was only structured in the  Table 1. Parameters of gradient models and variograms describing the spatial structures of the studied variables of black howler monkey groups. IFA-intraspecific index of food availability = × PC DBH, where PC is the average of the phenological scores. IDA: Interspecific index of food availability = PC IVI × , where IVI is the importance value index of trees. SF = Selected food (grams of dry weight). Food set = leaf and fruit intake (i.e. the sum of all items consumed). α, β and γ denote the estimate parameters for intercept, the north-south and east-west coordinates, respectively. ΔAIC = AIC null -the AIC selected model (AIC null is the AIC of the response variable explained by its mean). ni = not included in the model after model selection. A site and B site correspond to the G-II´s fragment division (North and south, respectively). See variograms in Supplementary  Fig. S3 www.nature.com/scientificreports www.nature.com/scientificreports/ G-I site in patches of 11.9 m in the "nortes" season, and according to IDA in patches of 17.5 m in the rainy season, with a similar aggregation intensity (R 2 = 0.3 and 0.4, respectively). Mature fruit patches occurred only in the dry season for both fragments and three sites ( Table 1). The dry season had the highest availability of mature fruit (Fig. 2), we found high degree of aggregation according to both indices (IFA: R 2 = 0.3-0.7, IDA: R 2 = 0.3-0.5). According to IFA, the mean mature fruit patch diameter was 18 m for the G-I site, 19.5 m for the G-IIA site, and 40.9 m for the G-IIB site. In contrast, IDA showed a mean patch diameter of 26.1 m for the G-I site, 18.9 m for the G-IIA site, and 37 m for the G-IIB site. Immature fruits were aggregated according to IFA in the G-I site during the dry season in patches of 57 m and of 6.7 m according to IDA. Furthermore, when we analyzed immature fruit availability throughout the year, they were aggregated according to IFA in patches of 20 m only in G-IIA, and according to IDA in patches of 17 m in G-IIA and 41 m in G-IIB. In summary, food was structured in patches more than in gradients and mature fruits in the dry season were very patchily distributed. For mature fruits, we found similar mean patch size values according to both availability indices.
We detected that during the dry season food set selection was aggregated with a mean patch diameter of 14 m (R 2 = 0.    (Table 1).
Linear modeling: food set selection explained by mature fruits availability. Food set selection in G-I and G-IIA was aggregated and explained by mature fruit aggregation, but only according IFA. As described above, using variograms we first detected that the food set selection in the dry season was aggregated in patches in G-I and G-IIA sites (mean diameter of 14 and 6.9 m; Fig. 4a,b, respectively), but not in G-IIB, as well as mature fruit availability (mean diameter of 18 and 19.5 m, Fig. 4c,d). Subsequently, in the G-I site, linear modeling showed that once the mean variation of the food set variable was explained by mature fruit availability ( Fig. 4e; GLS model parameters in Table 2), smaller patches persisted in the residuals (3 m mean diameter, Fig. 4g). For the G-IIA site, mature fruit availability also explained the patch aggregation associated to food set (Fig. 4f, parameters of the LM model in Table 2); however, the residuals of the model did not show any additional spatial structure (nugget model: Fig. 4h).

Discussion
We confirmed Hypothesis 1 by documenting that mature fruit productivity was aggregated in patches and this occurs in the peak season of availability. Furthermore, as suggested by the presence of gradients, mature fruits are also structured in a spatial pattern broader than the fragment size. The mean patch diameter associated to food set, leaf and fruit intake combined, during the dry season shows a diameter (6.9-14 m), which is in the same range of that associated to the mature fruits (18-19 m). This confirms Hypothesis 2. As predicted in Hypothesis 3, the spatial structure and mean variation of food set selection was explained by mature fruit patch structure, and not by their gradients. These results show that food selection by this primate is driven by the spatial distribution of fruits during certain seasons at small scale.
Trees showed intra and inter-specific variation in item productivity among seasons. The spatial modeling of both availability indices indicated that mature fruit was structured in space in the fragments during the same season. Hence, fruit productivity was aggregated, both due to the effect of large trees and to the influence of dominant plant species in the fragments. Supporting our first hypothesis, fruit production was structured in patches in the months of greatest production. However, when the availability was analyzed collapsing the seasonal values on a year basis, the intensity of aggregation was higher. This suggests that the aggregation of fruit availability sustains a certain degree of location stability throughout the year but fluctuates in intensity of production. It is also important to analyze both tree size and the contribution of dominant species in the analyses of resource spatial structures. The most dominant species were not equally dominant suppliers of fruits, except for Sabal mexicana (Supplementary Table S1). The most important species according IVI index in the G-I site were Cupressus lindleyi www.nature.com/scientificreports www.nature.com/scientificreports/ and Albizia leucocalyx; however, they did not represent a large part of the monkey's diet (3.1%.and11.9%, respectively). Likewise, the most important species in the G-II site, were Haematoxylon campechianum and Sabal mexicana, but only the second one was a strongly used as food resource. It is possible that the fruit patch structuring is due more to an aggregated distribution of the individual trees bearing fruits and less to the fruiting synchrony among dominant species or even among individuals of the same species. This is supported by the frequent IFA spatial aggregation and subtler and less frequent aggregation of IDA. Moreover, some species have an asynchronous production and are uncommon, such as Maclura tinctoria at G-I site, but they can be a key resource for these monkey groups in periods of scarcity 41 . Likewise, the patchy spatial pattern associated to fruit availability may be a response to habitat loss 1 . In our study site, the conversion of forests to pastures and crops in the last 50 years reduced native vegetation by 90%. This modified the spatial distribution of large trees, and thus the availability of food resources for howler monkeys 48 .
It has been suggested that the availability, distribution, and quality of the feeding patches have an important effect on primate social organization, group cohesion, and food choice 34,49,50 . Many primate species exhibit spatial knowledge of the environment and seemingly plan and travel routes to find a given resource area 51 . The concept of food patches is therefore central to most models of primate socioecology 52 . The definition of a food patch has always been controversial. This is largely because it has been done using different criteria and based on the observations of animal behavior 31,32 . Several studies have supported the use of DBH as an indicator of plant productivity and demonstrated that it is a good estimate of fruit abundance, hence it is used as an estimate of food patch size 20 . Alternatively, patches are defined as a continuum of food arranged in such a way that the forager can feed without interruption and may occur in an isolated tree or in a group of trees with contiguous canopies 31,53,54 . These are convenient operational definitions; however, these definitions result in a patch being defined according to the specific primate being studied and the specific questions addressed 8,55 . This is especially true when primates tend to spread out and simultaneously feed from multiple trees 55 . Consequently, other studies have attempted to improve upon traditional patch definitions by incorporating the individual perspective 8 , such as the focal-tree method. Some of them include spatial analyses to define food patches using grid cells 55 . Expanding the way in which a feeding patch is evaluated using the animal's perspective and linking this information to the spatial arrange of the resource, results in a more accurate assessment of the relationship between foragers and food distribution.
Under the traditional definition of food patch (DBH), it has been observed in black howlers that time spent feeding in a patch did not correlate with patch size 35 . In contrast, in mantled howler monkeys (A. palliata), there is a strong relationship between feeding time and the size of the fruit foraging patch 31,56 . However, this evidence does not allow to address whether the limits associated to foraging patches and its aggregation coincide with those of food items because these concepts have not been spatially and explicitly operationalized before 32 . The present study estimates patch diameters independently, one for the food resource and other for the food selection, based on modeling spatial autocorrelation. In support of our second hypothesis, the average patch diameter associated to local fruit availability (IFA 18-19 m) may represent several tree crowns and was assessed independently of the patch size associated to the food selection (6.9-14 m). This finding supports previous evidence that food selected by howlers often occur in clusters or patches and that they exploit larger fruit patches possibly due to its grouping behavior 53 . Similar conclusions have been described for other primates, such as black spider monkeys (Ateles paniscus chamek) and muriquis (Brachyteles amchnoides) 57,58 . Our study shows how an independent spatial modeling strategy of animal resources and behavior is able to operationalize the definitions of food patches and foraging, estimating the scales in which both are structured, and allowing their independent variation. This approach also subsequently allows testing hypotheses on the explanatory power of resource spatial distribution on behavior. www.nature.com/scientificreports www.nature.com/scientificreports/ The spatial structuring and slope associated to food set selection by howler monkeys is explained by the availability of one of the most preferred and limiting foods for primates, mature fruits (Hypothesis 3). This pattern occurred for both groups when fruit production was high, but not for young leaves. The howler's diet is not typically rich in lipids, but some fruits have high lipid content, which likely plays an important role in their food choice 59 . Although the spatial relationship has not been addressed in this way, previous studies have also observed that the consumption of fruits and flowers shows a positive correlation with the availability of fruits for howler monkeys 28,60 and in other primates such as gibbons (Hylobates lar) 61 and chimpanzees (Pan troglodytes) 62 .
Changes in the timing and intensity of fruit production have effects on the structure and function of the ecosystem 19 . For example, the density of frugivorous animals is negatively affected by a strong seasonality in the availability of fruits 18 . However, large animals (>5 kg) may be more tolerant to highly seasonal resources than smaller ones 63 . Studying the spatial structuring of food resources and evaluating whether the forager follows the preferred or limiting food could help explain this tolerance. Our findings suggest that for black howlers the tolerance to resource seasonality is due in part to the fact that they exploit the preferred or limiting food when it is available. However, the mixed consumption of leaves and fruits 59 , along with their small gut capacity 64 , prevent these primates from consuming large quantities of a food in a single day. Instead, they concentrated in consuming other foods in small quantities in the area where the preferred food is aggregated, probably according to its nutrient content, hence the spatial association in a food set basis. Several primates have a mixed foraging strategy in which some foods are encountered when animals travel. However, most of this movement is driven by the need to exploit a known set of patches in a near-optimal order 51 . It is noticeable that the food set selection exhibits spatial structuring, but when each item consumed is analyzed separately there is no spatial structuring pattern. The spatial statistics methods used here are very data demanding 16 and possibly the effect size when separating fruits and leaf data was not enough to detect the expected itemized spatial structures.
It has been suggested that the intraspecific variation in the nutritional content of the available food can generate spatial aggregation on a very small scale (i.e., between trees). This could further influence the spatial structuring of food selection by primates within available patches 26 . Such micro-structuring could be illustrated in our study by the residuals of the GLS model of the analysis of spatial dependence between food set selection and fruit availability were structured in smaller patches (3 m) in the G-I site. The residual spatial autocorrelation of a GLS model suggests the potential relevance of an explanatory process structured in space that was not previously considered 65 . This facilitates the generation of new ecological hypotheses. For example, the fact that the patch diameter associated to the food set was smaller than those of mature fruits, suggest that the monkeys select the trees, or crown zones, with the desired nutritional content within the patches where that resource is available.
It is essential to find the relevant scale to spatially study the relationships between animals and their environment because different processes are often at work, depending on the scale 66 . In some primates, biological processes are associated with specific scales 67,68 . For example, the population density of howlers studied at two spatial scales, 100 and 500 ha, is better explained by fragment size on a smaller spatial scale than on larger ones 9 . Consequently, conclusions derived from a spatial scale often cannot be extrapolated to others 66 . Our results suggest that feeding behavior of howlers operates at a small scale within the food patches and not at broader scales since the food selection did not follow the spatial pattern in gradients presented by the resources. It is known that fruiting changes are more pronounced at small scales due to seasonal and interannual variability, or to abnormally high fruiting of individual plants 21 . In addition, there is evidence that seasonal variation in the spatial availability of food supports environmentally driven changes in primate movement patterns 69 .
It is worth noting that neither food availability, nor selection of young and mature leaves, during the seasons of higher production, were structured in space. However, a spatial aggregation associated to the mature leaf availability in the rainy season was found. In general, leaves are more abundant and more evenly distributed than fruits, therefore, folivorous primates generally spend less time searching for food, rest more, and occupy smaller areas than frugivorous and insectivorous species 11,70 . Although there are reports in other howler monkey species (e.g., A. seniculus) that food availability and the consumption of leaf parts are correlated 71 , our results show that the howlers preferentially feed on the trees that have marked fluctuations in resource availability. That may be the reason for their differing patterns of leaf ' ingestion between localities, seasons, and years 35,64 . For example, for mantled howler monkeys (A. palliata) in coffee plantations in Nicaragua, the seasonal pattern of food selection did not follow the production pattern of the more abundant phenophases (i.e., mature leaves), nor of the most limiting resources (i.e. fruits) 56 . Hence, foods are consumed selectively, even though there are abundant potential resources in the habitat which are ignored.
In conclusion we found that there is a spatial aggregation of food resources at two scales, in patches nested in gradients, mainly fruits, and that this spatial structuring depends on the season. The spatial distribution in availability gradients had no explanatory power over the spatial distribution of food selection by howlers. However, the variance of food set selection was explained by the availability of mature fruit and the spatial structure of the food set selection was explained by the patchy distribution of fruit within the fragments during the season of higher production. Thus, howler monkeys probably base their foraging decisions on local spatial scales influenced by clusters of fruiting trees. This may maximize leaf and fruit intake in response to spatial and seasonal aggregation of the preferred resource.

Data availability
The datasets generated during and/or analyzed during the current study are available into the manuscript and the supplementary material.