Study site and subjects
From January 2017 to September 2018, we collected non-invasive observational data on two habituated groups of wild Southern pig-tailed macaques (Macaca nemestrina) at the Segari Melintang Forest Reserve (SMFR) and the oil palm plantations bordering its south-western edge (4°19–20′N, 100°34–36′E). To date, our study groups are the only wild population of this species that is successfully habituated, probably relating to the generally shy and elusive nature of M. nemestrina. The size and composition of the groups slightly changed during the study period, either due to male immigration or dispersal 70, animals dying or being born, or juveniles reaching sexual maturity. During the study period, group 1 (named AMY) consisted of 5 to 8 adult males, 12 to 15 adult females and 18 to 23 immature individuals. Group 2 (named VOL) consisted of 11 to 14 adult males, 19 to 21 adult females and 16 to 18 immature individuals. The age and sex of individual macaques were known from long-term observations 36,37. Both groups visited the plantation area bordering their forest habitat almost daily (mean ± SD (AMY/VOL) = 3.1 ± 1.8/2.7 ± 1.8 hours per day 37). The annual home ranges of group AMY and VOL were 92.7 and 96.6 hectares, respectively, with used plantation areas accounting for approximately one third of the total home range areas 37. As group VOL has not been fully habituated before the start of 2018, assessments of the macaques’ social network and the mother-infant relationship were performed only on group AMY.
Habitat types
In order to assess the impact of anthropogenic environments on the macaques’ social activities, we divided the home range areas of our study groups into three habitat types, i.e. forest, plantation edge and plantation (Supplementary Figure S1). SMFR comprises 2742 hectares of which 408 hectares are strictly protected Virgin Jungle Reserve. Its main vegetation types are dipterocarp lowland forest and alluvial fresh-water swamp 36. The 420-hectare sized oil palm plantation bordering the reserve was established between 1980 and 1990 and is managed by a federal authority. The oil palm estate was accessible to macaques, with encounters between wildlife and plantation workers being occasionally observed, yet these did not involve regular conflicts, such as hunting or chasing. We defined ‘plantation edge’ as the plantation area which is located within 50 metres from the forest border, whereas ‘plantation’ refers to all plantation areas further than 50 metres from the forest border. This distinction was made to account for the fact that plantation areas in close proximity to the forest provide additional shelter and protection for the macaques through close-by forest vegetation. So-called ecotones that form transitional areas between two distinct ecological habitats were reported to be of great environmental importance, potentially serving as speciation and biodiversity centres 71. We chose the distance of 50 metres according to the average diameter of the macaque groups’ dispersion (edge-centre-edge).
Behavioural data collection
We collected data using 30-minute focal animal sampling 72 based on a species-specific ethogram established for the study species (adapted from 73) in the forest, at the plantation edge and inside the plantation. We observed a total of 50 individually recognizable macaques (36 of group AMY, 14 of group VOL). Focal individuals were chosen to represent all age-sex classes. The order of focal observations was randomized, aiming at sampling each individual only once per day. If a focal animal entered another habitat type during a 30-minute sampling protocol or went out of sight for more than ten minutes, this observation was stopped. Incomplete protocols were considered for multivariate analysis if they lasted at least 15 minutes. Total observation time was 724 hours (mean ± SD = 14.5 ± 3.6 hours per subject).
To assess socio-ecological risks posed by forest and plantation habitats, we recorded the frequency and duration of flight behaviour shown by the macaques in response to the presence of potential predators, i.e. humans, feral dogs, birds of prey, snakes or monitor lizards. To provide an overview on shifts in the macaques’ overall time budget, we continuously recorded their activity in the different habitats. Recorded activities included feeding (i.e. ingesting food), foraging (i.e. searching for or manipulating food), locomotion (i.e. walking, running, climbing or travelling without other activity), resting (i.e. standing, sitting or lying without other activity, eyes may be open or closed), socializing (i.e. all positive social interactions, e.g., grooming, being groomed and groom presenting, social play and huddling) and others (e.g. sexual and agonistic interactions or self-grooming). To assess individual differences in affiliation, we recorded the frequency, measured as rate per hour, and duration of all bouts of grooming and social play between the focal subject and other group members. As measures of frequency and duration were highly correlated (Pearson’s r (grooming/social play) = 0.78/0.81, P < 0.001), we considered only frequencies for analyses. Further, we recorded aggressive behaviour exchanged between the focal subject and other group members, considering both physical (i.e. attack, bite, grab, hit, push) and non-physical aggression (i.e. charge, chase, lunge, stare and vocal or open mouth threat). Social data were complemented by ad libitum data 72 on aggression, displacement and submission among adult males and adult females for the purpose of constructing dominance hierarchies (see below). Data on social interactions included information on both the initiator and the recipient. Following previous studies 74, a repetition of a behaviour was scored as a new bout if more than ten seconds had elapsed between occurrences or at least one partner had switched to a mutually exclusive activity (e.g., from grooming to aggression). During an aggressive event in which a number of different agonistic patterns occurred in quick succession, only the most intense kind of aggression was considered for analyses 74.
To assess affiliative social networks across different habitats, we recorded data on spatial proximity between macaques. We took point time scans 72 every three minutes within the 30-minute sampling protocol, recording all group members in body contact with the focal individual. We further recorded whether or not this contact resulted from an affiliative interaction (e.g., during grooming, play or huddling). This was the case for 98.3% of our observations. The total number of scans recorded was 14,205 (mean ± SD = 284 ± 71 scans per subject).
To assess the mother-infant relationship, we additionally observed eleven mother-infant pairs from group AMY in the three different habitats for the first six months after infant birth. Total observation time was 240 hours (mean ± SD = 21.8 ± 9.4 hours per mother-infant pair). We continuously recorded maternal behaviour promoting infant independence 46. Specifically, we recorded the number of contacts broken (i.e. any movements disrupting body contact between mother and infant), increases of distance (i.e. movements increasing the distance between mother and infant from within arm’s reach (about 60 cm) to outside of arm’s reach) and maternal rejection (i.e. attempts by the infant to make contact that were rejected by the mother, e.g., by turning, running away, or holding the infant at a distance with an arm 46). To ensure independence between these measures, increases of distance were only recorded if at least five seconds elapsed since contacts were broken. To assess spatial proximity in mother-infant pairs, we took point time scans 72 every minute during focal observations, recording whether or not mothers and their infants stayed in body contact, including ventro-ventral contact, nipple holding, side-by-side contact, and grooming.
Individual location data
We collected individual location data with a Garmin GPSMAP62s daily during behavioural observations, with the coordinates of each focal observation being taken at half-time of the respective focal protocol. Annual home range areas of group AMY and VOL were adopted from Holzner et al. 37, showing point Kernel Density Estimations (KDE) with 95% probability of use 75. To provide an overview on the occurrence of affiliative and aggressive social interactions across different habitats within the macaques’ home ranges, we created interpolation maps (see Fig. 2) based on mean behavioural rates occurring during focal observations per 50 m x 50 m grid cell using the Inverse Distance Weighting (IDW) tool of the software QGIS version 3.12 (QGIS Development Team 2020).
Dominance hierarchy
From 948 dyadic agonistic interactions with a clear winner and loser outcome collected during focal and ad libitum observations, we estimated rank orders using David’s scores 76. These were obtained in R version 3.4.4 (R Core Team 2020) using the function DS from the package ‘EloRating’ (version 0.46.8 77). We set argument ‘prop’ to ‘Dij’, calculating dyadic win proportions corrected for chance 78. As in macaques rank acquisition and function typically differ between sexes, with non-natal males fighting for dominance, while females socially inherit the rank of their mothers 79, we estimated rank orders separately for males and females. Following Kaburu et al. 80, we controlled for differences in group size and sex ratio by standardizing dominance rank as ‘(Rank-1)/(N-1)’, where N represents the number of focal animals per group and sex class. Standardized dominance rank values range from zero (top-ranking individual) to one (bottom-ranking individual). Referring to literature 81, immature males and females got assigned the same rank as their biological mother, or, if their mother already died, the rank of their closest adult female relative.
As David’s scores do not account for potential temporal fluctuations in rank orders, we further assessed rank stability using the function stab.elo of the package ‘EloRating’ (version 0.46.8 77). Indices range between zero (unstable) and one (stable), reflecting to what extent rank positions of individuals change over time. Males and females of both groups displayed stability indices close to one (males (AMY) = 0.9950, females (AMY) = 0.9956, males (VOL) = 0.9909, females (VOL) = 0.9943), suggesting highly stable dominance hierarchies during the sampling period in both sexes. Therefore, David’s scores seem to be appropriate for estimating dominance ranks in our study groups.
Social network analysis
Based on affiliative interactions observed during individual focal sampling, we constructed the social network of group AMY separately for forest and plantation habitats. Following Lehmann et al. 29, we assessed dyadic affiliation as the proportion of scans two individuals were in social contact (i.e. grooming, social play or affiliative body contact). We created social networks in R version 3.4.4 (R Core Team 2020) using an undirected data structure with the function graph_from_data_frame from the package ‘igraph’ (version 1.2.5 82). For each individual, we extracted the binary degree and eigenvector centrality (EC), two commonly used network parameters to quantify individual social connectedness 39,29. The binary degree reflects an individual’s overall number of interaction partners, while EC is a measure of both direct and indirect network ties, reflecting a node’s importance while considering the importance of its neighbours. Thus, a high value of EC suggests that an individual has many social partners who themselves also have many partners. While considering raw counts for binary degree, with regard to EC we were particularly interested in an individual’s connectedness in relation to other group members. We therefore rescaled the obtained values of EC in each habitat to get percentile scores lying between zero (minimum) and one (maximum).
Statistical Analysis
Multivariate statistical analyses assessing the impact of different habitats on the frequencies of affiliation and aggression (models 1 to 4), social network connectedness (models 5 and 6), and the mother-infant relationship (models 7 to 9) were conducted in R version 3.4.4 (R Core Team 2020) using Generalized Linear Mixed Models (GLMM 83). Models were fitted using the functions lmer and glmer of the package ‘lme4’ (version 1.1.19 84) with the optimizer bobyqa. To facilitate model interpretation and convergence, we standardized all continuous predictors before model fitting to a mean of zero and a standard deviation of one 85. Full-null model comparisons were derived using likelihood ratio tests (LRT 86) using the R function anova with argument ‘test’ set to ‘Chisq’. Tests of individual fixed effects were derived using the R function drop1 with argument ‘test’ set to ‘Chisq’, yet control predictors were not interpreted. Confidence intervals were derived using the function bootMer of the package ‘lme4’ (version 1.1.19 84), using 1,000 parametric bootstraps.
Frequencies of affiliation and aggression (models 1 to 4): To investigate the impact of the habitat on the display of grooming, juvenile social play, and physical and non-physical aggression, we constructed four GLMMs 83 with Poisson error structure and log link function. As response variables we used the number of grooming bouts (model 1), bouts of juvenile social play (model 2), bouts of physical aggression (model 3), and bouts of non-physical aggression (model 4) per focal observation (N = 1,535 focal observations of 50 individuals for models 1, 3 and 4, and N = 510 focal observations of 16 immature individuals for model 2). The models included the habitat (forest, plantation edge or plantation) as fixed effect test predictor, while controlling for individual dominance rank and age-sex class (adult male, adult female, immature male or immature female), as both rank and age-sex class have previously been shown to affect the occurrence of affiliative and agonistic interactions in macaques 27,57,87. To account for changes in the overall group activity throughout the day, inter-group variation and repeated observations of the same individuals, we further included the time of the day, divided into four time blocks (early morning: 7:00–09:59 am, late morning: 10:00 am − 12:59 pm, early afternoon: 13:00–15:59 pm or late afternoon: 16:00–18:59 pm) and macaque group (AMY or VOL) as fixed effect control predictors, and the focal individual ID and sampling date as random effects. Additionally, we included the random slopes of habitat and time of the day within focal individual in models 1 to 4, and the random slope of rank within sampling date in models 1, 3 and 4 86,88. Controlling for differences in the sampling effort, we further included the duration of each focal observation as an offset term into the models 89. To test the effect of different habitats, we compared the full models with respective reduced models lacking only our test predictor (habitat) using a LRT 86.
Social network connectedness (models 5 and 6): To investigate the impact of the habitat on two common network parameters, i.e. the binary degree and EC defined above, we constructed two GLMMs 83 with Poisson and Gaussian error structure, respectively. As response variables we used the individuals’ binary degree (model 5) and scaled EC (model 6) in each habitat (N = 68 observations of 34 individuals). We included the habitat (forest or plantation edge) and its interactions with individual dominance rank and age-sex class (as defined above) as fixed effects and the focal individual ID as random effect. Controlling for differences in the sampling effort, we further included the total number of point time scans per individual as an offset term into model 5 89. To test the effect of different habitats on the binary degree and EC, we performed LRTs 86, comparing the full models with respective reduced models lacking our test predictor (habitat) and its interactions with dominance rank and age-sex class.
To account for interdependency of network measures used as outcome variables in our GLMMs, we used a node-swapping permutation approach, based on 1,000 permutations of the outcome variables 90. This included recalculating the network parameters after randomly swapping the nodes of the original networks. We used node-swapping (as opposed to generating random graphs or using pre-network ‘edge-swapping’ randomizations) since this approach seemed better suited for our purposes of testing regression-based null hypotheses in a taxon with a largely stable group composition and relatively low observation biases 44,45. Specifically, node-swapping preserves the overall size, number of connections, and structure of the network, thereby also preserving the overall distribution of node-level measures such as degree and EC. It is therefore a more conservative approach that may be less susceptible to Type I errors, compared to random graph generation or edge-swapping 44,45. After each permuted swapping event, we refit the same GLMM using these newly created parameters as response variable. Comparing the observed model coefficients with the distribution of the 1,000 permuted coefficients, we calculated p-values as the number of times the coefficient value of the observed network is smaller than a randomized network, divided by the number of randomizations 90.
Mother-infant relationship (models 7 to 9): To investigate the impact of the habitat on the mother-infant relationship, we constructed three GLMMs 83, with the proportion of body contact between mothers and offspring (model 7), the number of maternal breaks of contact (model 8), and the number of maternal increases of distance (model 9) per focal observation being the response variables (N = 491 focal observations of 11 mother-infant pairs for models 7–9). Investigating effects on the proportion of time spent in contact (model 7), we used a GLMM 83 with binomial error structure and logit link function. In R, this analysis of proportions is possible by using a two-columns matrix with the number of contacts and non-contacts per individual as the response 83. Models 8 and 9 were created using a count response with Poisson error structure and log link function. Here, we controlled for differences in the sampling effort by including the duration of each focal observation as an offset term 89. In all three models, we included the habitat (forest, plantation edge or plantation) as a fixed effect test predictor, while controlling for infant and maternal characteristics, which were previously shown to affect the mother-infant bond, i.e. infant age 46 and sex (male or female 33), as well as maternal rank and parity (primiparous or multiparous 31,32). As in models 1–4, we accounted for changes in the overall group activity over the day by including the time of the day as fixed effect control predictor. Further, we included the mother-infant pair and sampling date as well as the combination of these two as random effects, as mother-infant pairs were frequently observed more than once on a given day. Additionally, we included the random slopes of habitat, infant age and time of the day within the mother-infant pair 86,88. As we expected infant age to have a non-linear effect on the rates of maternal breaking contact and increasing distance, we additionally incorporated squared infant age into models 8 and 9. Further, we included the two-way interaction between habitat and infant age in model 7, and its interactions with infant age and squared infant age in models 8 and 9. To test the effect of different habitats, we compared the full models with the respective reduced models lacking our test predictor (habitat) and its interactions with infant age and squared infant age, respectively, using LRTs 86. In case of a non-significant interaction, we re-ran the model without the interaction term to facilitate the interpretation of the main effects in the model.
For models 1 to 9, we assessed model stability by excluding the levels of the random effects one at a time and comparing the estimates from the obtained models with the estimates from the model based on all data using a self-written R function provided by Roger Mundry. This did not indicate any obviously influential case. To rule out collinearity, we determined Variance Inflation Factors (VIFs) for respective standard linear models excluding the random effects using the function vif of the package ‘car’ (version 3.0.2 91). According to these, none of the models indicated issues regarding collinearity (all VIFs < 1.16). We confirmed that overdispersion was no issue in all models with a Poisson (models 1–4, 5, 8, 9) or binomial (model 7) error structure. Yet, there was an indication of underdispersion in model 3, probably reflecting the low rates of physical aggression in our dataset. Further, there was indication of slight underdispersion in models 7–9, likely due to limited data points (dispersion parameters = 0.98/ 0.95/ 0.45/ 0.77 for models 1–4, 0.95 for model 5 and 0.65/ 0.67/ 0.52 for models 7–9). Visual inspections of a QQ-plot 92 of the residuals and a scatterplot of the residuals plotted against the fitted values 93 did not reveal obvious deviations from the assumptions of normally distributed and homogeneous residuals in the Gaussian model (model 6).
Ethical note
We obtained permits to study Macaca nemestrina from the Department of Wildlife and National Parks Peninsular Malaysia (permit holder: Dean of School of Biological Sciences, Universiti Sains Malaysia). We obtained permits to enter the forest reserve bordering the oil palm plantation from the Forestry Department Peninsular Malaysia (permit holder: Asyraf Mansor, School of Biological Sciences, Universiti Sains Malaysia). No written permit was needed to enter the plantations, but we informed the local management about the study. This non-invasive study was conducted in line with Universiti Sains Malaysia’s animal ethics requirements.