Sulawesi Crested Macaque (Macaca nigra) Grooming Networks Are Robust to Perturbation While Individual Associations Are More Labile

Social and environmental disturbance occurs naturally, and species in bonded social groups should be resilient to it. Empirical evidence of social responses to disturbance in primates, however, remains limited. We constructed social networks using group-level scan samples (N = 299) to test the robustness of grooming networks in a captive group of 20 Sulawesi crested macaques (Macaca nigra) to two management interventions involving environmental and social disturbance. During the first, the institution removed six castrated males and one female, contracepted six of the nine remaining females, and moved the group to a new enclosure. The second involved the introduction of a novel, reproductive male five weeks later. Networks remained stable following the first intervention. However, after introduction of the male, the number of grooming partners and the frequency of grooming with non-maternal kin increased in female-only networks. We observed less marked increases in the grooming frequency and number of grooming partners in whole group networks. Ten weeks later, network structure was more similar to that of pre-intervention networks than post-intervention networks. Our results suggest that reproductive males play a more important role in structuring Sulawesi crested macaque social networks than castrated males, as networks expanded and relationships between non-maternal kin occurred more frequently after introduction of the reproductive male. However, network responses to interventions appeared to be temporary as networks following a period of acclimation more closely resembled pre-intervention networks than post-intervention networks. Our study demonstrates the utility of social network analysis for understanding the impact of disturbance on stable social groups.


Introduction
Social organization results from the interactions between group members (Hinde 1976) and has important consequences for the emergence of cooperative behavior, disease transmission, and reproductive success (Dunbar and Shultz 2010;Leu et al. 2016;McCowan et al. 2008;Shultz et al. 2011;Weber et al. 2013). In most social groups, interactions between conspecifics are not random; individuals can choose with whom to interact with and for how long, choices that vary depending on social system, hierarchy, and kinship (Formica et al. 2012;Silk et al. 2009). Social relationships, particularly among kin, provide support during periods of social stress (Abbott et al. 2003;Young et al. 2014), and strong relationships among females encourage coalitionary support during conflict and cooperative antipredator behavior (Micheletta et al. 2012;Sterck et al. 1997). Moreover, intersex social relationships can shield offspring from male harassment (Smuts 1985). Ultimately, social integration has fitness benefits through its association with increased infant survival and lifespan, an association that persists independently of maternal rank (Silk et al. 2003(Silk et al. , 2009. As the structure of a social group affects its function in a given environment (Strogatz 2001), and gregarious, long-lived species rely on social connections and information transfer to survive and exploit their habitat (Lusseau 2003), social groups should develop mechanisms to persist through changing landscapes to prevent network disruption or fragmentation.
Social stability provides physiological and psychological benefits to individuals in a social group (Sapolsky 1983(Sapolsky , 2005, and the disruption of social associations may affect group cohesion and network function (Flack et al. 2005(Flack et al. , 2006Sapolsky 1983). Moreover, there are fitness costs of instability as disruptions to a social network can affect reproductive success, communication, or resource acquisition (Silk et al. 2003;Wey et al. 2008). Existing literature suggests that social networks can remain stable following the removal of individuals through two processes: network plasticity or network robustness. Plasticity describes the reconfiguration of associations in a network (Barrett et al. 2012;Stroeymeyt et al. 2018), while robustness suggests that the removed links were redundant to the network, and that remaining links can maintain the structure and function of a network despite the removal of individuals (Albert et al. 2000;Lusseau 2003). In most social networks, only a few key individuals have a higher number of connections than the mean for that network (Barabási and Albert 1999). The removal of these individuals can result in the loss of network structure, and networks could destabilize or fragment (Flack et al. 2006). In contrast, the addition of novel individuals to the network may result in unpredictable consequences for network structure (Barabási et al. 2002).
Social network analysis describes the strength and nature of social relationships by characterizing a population as a sum of direct and indirect associations between individuals in the group and by describing group cohesion and information flow (Brent et al. 2011). In primates, previous research highlights the variation of social network responses to social and environmental disturbance. In wild savannah baboons (Papio cynocephalus), for example, the natural death or dispersal of high-ranking males had minimal short-term and no long-term effects on female social networks , while social networks in wild female chacma baboons (Papio hamadryas ursinus) expanded for 3 months after the loss of a close relative (Engh et al. 2006a). In captive primates, the social networks of pigtailed macaques (Macaca nemestrina) became smaller, less cohesive, and less diverse after the temporary removal of high-ranking males (Flack et al. 2006), while social networks in captive squirrel monkeys (Saimiri sciureus) and capuchins (Cebus apella) were more cohesive in the 2 days following an enclosure move (Dufour et al. 2011). Different responses to disturbance may result from the type of disturbance that occurred, as well as the timescale of observations . Captive primates may also have more extreme responses to disturbance than those in the wild, as the diversity of interaction partners is limited and social hierarchies are more easily enforced (Stevens et al. 2008). Moreover, as captive breeding programs transfer individuals primarily based on their genetic attributes, it is possible that animal managers remove individuals that function as integral components for the maintenance of social cohesion, further augmenting responses to social disturbance in captive primates when compared to primates in the wild (Flack et al. 2006;McCowan et al. 2008).
Sulawesi crested macaques (Macaca nigra) reside in large multimale, multifemale groups and are a highly socially tolerant macaque species characterised by low levels of intense aggression and a high tendency to reconcile (Duboscq et al. 2013;Riley 2010;Thierry et al. 2000). In this opportunistic study, we examined the impact of large-scale social and environmental management interventions on the social networks of captive Sulawesi crested macaques housed at Chester Zoo. We collected data from July to December 2015, during which time the group was subject to two management interventions: during the first, the zoo removed seven adult group members, implanted six of the nine remaining females with hormonal contraception, and moved the remaining group to a novel enclosure. During the second intervention, five weeks later, the zoo introduced a novel male to the group. We focussed our analyses on the stability of grooming associations, as grooming has a predominantly social rather than hygienic function in primates (Dunbar 1991;Henzi and Barrett 1999) and is used to reinforce social relationships, maintain social tolerance, and reconcile following conflict (Borgeaud and Bshary 2015;Henzi and Barrett 1999).
As Sulawesi crested macaques are highly socially tolerant, individuals in stable social groups should interact with all group members equally (Duboscq et al. 2013;Thierry 2000). As such, we predicted that baseline social networks were well connected. If social relationships help buffer responses to social stressors (Abbott et al. 2003;Swedell 2012;Wittig et al. 2008) and help reinforce social relationships (Henzi and Barrett 1999), we predicted that the strength of grooming associations and number of grooming partners would increase following each management intervention. Moreover, if social stability and cohesion convey fitness benefits (Silk et al. 2003(Silk et al. , 2009Wey et al. 2008), we predicted that following a period of acclimation to the novel social and environmental conditions, social networks would regain their baseline structure. Lastly, we examined the influence of disturbance on kin relationships. If kin relationships help mitigate stress responses to social stressors (Abbott et al. 2003) and increase in populations with a higher proportion of females , we predicted that grooming associations would correlate more with relatedness following both management interventions than during the baseline observation period.

Study Subjects and Observation Periods
We collected data at Chester Zoo between the hours of 09:00 and 15:00 from July to December 2015. We divided the study into four observation periods, each lasting five weeks: a baseline observation period prior to management interventions; an observation period following the removal of seven adults (six males and one female), contraception of six of the nine remaining females, and the transfer of the remaining group to a novel enclosure (post-intervention 1); an observation period following the introduction of a novel, intact male to the group (post-intervention 2); and a terminal period 10-15 weeks after the introduction of the novel male, by which time the macaques also had access to an outdoor enclosure. We collected data for the terminal period 10 weeks after the male introduction to determine whether networks revert to their baseline state.
The initial group consisted of 18 adults (8 castrated males and 10 females) and 2 juvenile females, all of which were natal to the group. Group size varied from 20 individuals during the baseline, to 13 individuals following the first intervention, and 14 individuals following the second intervention and during the terminal observation period. The number of females decreased by only one individual throughout the study (baseline: N = 12; post-intervention 1, post-intervention 2, and terminal: N = 11).
During the enclosure move, the area of the enclosure decreased from 191.5 m 2 to 186.5 m 2 , and enclosure height increased from 5.7 m to 6.3 m-8.5 m. Group density decreased by only 0.02-0.03 individuals/m 2 between observation periods (baseline: 0.10 individuals/m 2 , post-intervention 1: 0.07 individuals/m 2 , post-intervention 2: 0.08 individuals/m 2 ). Although the enclosure sizes in our study are smaller than those in comparable studies (754.76 m 2 , Flack et al. 2006;2023.43

Data Collection
We collected data using group instantaneous scan sampling (Altmann 1974) every 30 min throughout the observation period. We collected 80 baseline (mean ± standard deviation (SD) = 16 ± 6 samples per week), 80 post-intervention 1 (16 ± 4), 71 postintervention 2 (12 ± 2), and 68 terminal scan samples (10 ± 2). We assumed that scan samples were independent, as for most behaviors (excluding resting) scan samples at an interval of 15 min maintain independence (Poulsen et al. 2001). Our analyses focus on the effects of management interventions on allogrooming (henceforth grooming) networks. We examined changes to the structure of whole group and female-only networks. Whole group networks include all individuals present in each observation period, while female-only networks focus on the females that were present throughout the study (N = 11; the female that was removed during the first intervention was omitted from baseline matrices). Female-only networks included the two juvenile females, as excluding juveniles from social networks can significantly alter social network structure (Fedurek and Lehmann 2017).

Social Network Analysis
We constructed binary association matrices for each observation period using presence (1) or absence (0) scores for interactions during each set of scan samples. We used these association matrices to construct all social networks. All networks were undirected (symmetric), weighted, and did not include self-loops (an interaction with oneself, e.g., autogrooming). By using undirected networks, we maximized network densities and minimized the number of network metrics analyzed (removing in-and out-measures) (Lehmann et al. 2016). We weighted networks using the simple ratio index (SRI) (Cairns and Schwager 1987). We calculated the SRI by dividing the total number of times we observed individuals a and b interacting by the total number of times we observed individuals a and b interacting plus the total number of times we observed individuals a and b separately.
We evaluated differences in overall network structure following social and environmental management interventions using group-level measures of graph density, network diameter, and assortativity. Graph density describes how well connected a network is and calculates the proportion of potential connections between individuals that are actually connected. Network diameter is an indicator of how large a social network is and gives the value of the shortest distance between the two most distant individuals. Lastly, network assortment (assortativity) describes the tendency of similar individuals to be connected. For our analyses, we investigated assortment by network degree. High assortativity indicates that individuals prefer to interact with individuals of a similar degree and with a similar access to social resources, decreasing the diversity of interaction partners in a group (Flack et al. 2006).
Furthermore, we calculated several individual-level network metrics: binary degree, graph strength (weighted degree), harmonic centrality, and betweenness centrality. As direct measures of an individual's involvement in the network, the number and weight of connections that an individual influences an individual's binary degree, graph strength, and harmonic centrality, whereas betweenness centrality, an indirect measure of network involvement, relies on the interactions between an individual's neighbors.
Binary degree describes the number of individuals the focal individual interacts with. Graph strength, also known as the weighted degree, is the sum of an individual's edge weights and reflects the number of times an individual interacted with others. Closeness centrality measures the closeness of one individual in the network to all others. This measure, however, considers edge weights as distances, rather than the strength of the association. We therefore used a modified closeness centrality, the harmonic centrality, which calculates centrality based on the sum of reciprocal distances between individuals rather than the reciprocal of the sum of distances. Given that self-loops are not considered in these networks, scores are then normalized by group size (N − 1), bounding the score between 0 (isolated individual) and 1 (completely connected to every other individual in the network) (Opsahl et al. 2010;Rochat 2009).
Betweenness centrality is the number of shortest paths between two individuals going through the focal individual and assesses the importance of individuals in overall network cohesion. As the networks are weighted, we use a weighted variant of traditional betweenness centrality in which we consider the weight of the ties between individuals as well as the number of shortest paths positively, and in equal proportions (tnet α = 0.5) (Opsahl et al. 2010). Individuals with high betweenness values are important for network connectivity and have the potential to disconnect the network if removed (Bret et al. 2013;Kanngiesser et al. 2011).
We used a Quadratic Assignment Procedure (QAP) test to determine the extent to which grooming associations correlate with relatedness. The QAP test computes the Pearson's correlation coefficient (C) between two matrices before randomly rearranging the columns and rows of the matrices to recalculate the correlation coefficients. These permutations are used to estimate the significance of the correlation by calculating the proportion of permutations (P r ) in which the correlation coefficient of random matrices is larger than or equal to the observed correlation (P r ≥ C) (Krackhardt 1987). Statistical significance at P = 0.05 is implied by a low proportion (P r ≤ 0.05) (Imbert et al. 2016). We compiled data on relatedness using pedigree information from the Zoological Information Management System (ZIMS) database. We had complete information for matrilines but did not have full paternity data for the entire group; therefore, we could only evaluate matrilineal relatedness (Duboscq et al. 2017;Whitehouse et al. 2017). As such, the relatedness coefficient between a mother and her offspring is 0.5, between siblings is 0.25, and between an aunt and niece is 0.125. Cercopithecine primates are female bonded and the strongest, and most beneficial relationships in groups are normally between close female kin (Cords 2012;Silk et al. 2009). Excluding paternal relatedness from our analyses may mask some additional kin effects on the social responses of the Sulawesi crested macaques, as paternal relatives in many species interact more frequently than non-kin (Schülke et al. 2013;Schülke and Ostner 2008;Smith et al. 2003;Widdig et al. 2002).
To assess the persistence of associations between individuals within the network, we used QAP tests to assess the correlation of associations between weighted association matrices. Initial analyses examined the correlation between the baseline, post-intervention 1, and post-intervention 2 matrices as well as the correlation between the terminal, baseline, and post-intervention 2 matrices. As matrices used in QAP tests must contain the same individuals, we could only determine correlations between observation periods for networks with the same group composition. We therefore excluded the comparison between terminal and post-intervention 2 matrices.
Lastly, we assessed whether changes to network structure following the first intervention resulted from the removal of individuals or from changes to individual associations. To do this, we created knockout models by recalculating baseline networks after the exclusion of the removed individuals, following methods presented by Flack et al. (2006). As female networks remained consistent throughout the observation period, we created knockout models only for whole group networks.

Statistical Analysis
Network Permutations and Linear Models Social network metrics vary with observation effort and sampling frequency, e.g., if the observer is able to observe certain individuals more easily and more frequently than others (Farine and Whitehead 2015). Moreover, individual-level network data are inherently non-independent, as social network analysis counts all associations for a focal individual for all individuals connected to that individual too. Social network data therefore violate the assumptions of traditional statistical analyses (Croft et al. 2011;Farine and Whitehead 2015). These biases are controlled for by using null models constructed using permutations of the underlying network data (Croft et al. 2011;Farine 2017;Lusseau et al. 2008).
In this study, network permutations involve sequentially swapping associations between individuals in the binary matrix for each sampling period (observation period). Following each swap, new edge weights for the network are calculated. To assess whether a given network metric for a network is statistically different from a network in which individuals are randomly associating, we calculated the proportion of permuted network values (P r ) that are larger or smaller than the empirical network value. As there was no a priori prediction for the direction of change between the metrics of the observed and permuted networks, all analyses evaluate two-tailed significance, with significant differences found when 2.5% of the permuted values are more extreme than the empirical value (P r < 0.025 or P r > 0.975) (Blaszczyk 2018;Farine 2017). As such, we calculated significance at P ≤ 0.05 as twice the proportion of permuted values that are larger or smaller than the empirical value. To ensure that the permutation-based Pvalue had stabilized, we permuted networks 2,000 times.
We assessed the effect of observation period on each of our individual-level network metrics using linear mixed models (LMMs) in which the network metric was the dependent variable, observation period the independent variable, and the identity of the individual was set as the random factor. As network variables did not correlate strongly with one another, we included all individual level network metrics in analyses (Electronic Supplementary Material [ESM] Table SII). As we included multiple observation periods in each LMM, we calculated differences in the mean network metric across observation periods using methods analogous to the Tukey test, modified to incorporate network permutations. We calculated the slope between observation periods from empirical networks and compared this value to those calculated from network permutations to assess whether or not differences between the networks are due to chance (Farine 2013;Farine and Whitehead 2015). Initial analyses that examined the effect of management interventions included the baseline, post-intervention 1, and post-intervention 2 networks. To control for multiple comparisons, we adjusted P-values using Bonferroni correction, resulting in an α-significance of 0.02.
Following initial analyses, we used the modified Tukey test to compare differences between terminal networks and the baseline and post-intervention 2 networks to examine network structure following a period of acclimation. These models did not include post-intervention 2 networks, as group composition differed between networks and there was no period of acclimation between the first and second management interventions. To control for multiple comparisons, we adjusted P-values using Bonferroni correction, resulting in an α-significance of 0.03. Lastly, we used the modified Tukey test to compare knockout networks to baseline and post-intervention 1 networks. Again, we did not compare networks following the second intervention to the terminal period from these analyses due to differences in group composition between networks, and the α-significance for these comparisons was 0.03.
To account for differences in group size, we scaled several of our network metrics. As we calculated network diameter as a function of the distance between individuals, which positively correlates with group size (Nunn et al. 2015), we divided network diameter by group size. Moreover, as network assortment reflects differences in group size (Barrat et al. 2008;Hetman et al. 2008), we also divided network assortment by group size. As the number of individuals the focal individual is interacting with influences the binary degree and graph strength, we divided binary degree and graph strength by the total number of potential individuals the focal individual could interact with (group size -1) (Freeman 1979;Kasper and Voelkl 2009). Lastly, we bound each individual's betweenness and harmonic centrality between 0 and 1 by dividing each value by the highest betweenness score in that observation period so that values were relative to one another (Borgatti et al. 2018;Wasserman and Faust 1994). We ensured that each model met the assumptions of linear models by examining the distribution of residuals as well as the residual plots.
To determine whether changes to the social group affect all individuals uniformly or whether certain individuals are disproportionately affected, we compared the relative rank of individuals in subsequent observation periods (baselinepost-intervention 1, post-intervention 1post-intervention 2, postintervention 2terminal) (Farine 2018) using the network permutation method described in the foregoing. We determined an individual's rank for each individual-level network metric and normalized this value by dividing each rank by the highest rank in that observation period to ensure that we could compare rank across observation periods. We included individuals in relative rank analyses only if they were present in all observation periods. We adjusted P-values using Bonferroni correction, resulting in an α-significance of 0.02.
The removal of individuals from a group can have significant effects on network structure, depending on their role in structuring social associations. To assess whether the network metrics of removed individuals differed significantly from individuals that remained in the group, we compared the network metrics of the removed individuals to those of the remaining individuals using Welch's two-sample t-tests. We calculated statistical significance by comparing the proportion of permuted t-values that were more extreme than the observed tvalue. As the male introduction involved only one individual, a similar comparison was not possible for the second management intervention. Instead, we present the percentile in which the network value for the introduced male lies within values for the whole group.

Statistical Software
We ran all analyses using R (R Development Core Team 2013). We used asnipe to construct and permute social networks (Farine 2013); igraph to calculate graph density, graph strength, and harmonic centrality (Csardi and Nepusz 2006); tnet to calculate betweenness centrality and network diameter (Opsahl et al. 2010); lme4 to run LMMs; and multcomp to calculate conditional R 2 values (R 2 c ), which incorporate random effects. We used assortnet to calculate assortativity (Farine 2016) and sna to carry out QAP tests (Butts 2016).

Ethical Note
The interventions described in this research were part of scheduled, management interventions for the Sulawesi crested macaques at Chester Zoo. Animal care staff trained all animals to enter transport crates before moving them from one enclosure to the other. The removal of individuals occurred for welfare reasons, and the animal care staff followed appropriate protocols to minimize stress during the procedure. The authors declare that they have no conflict of interest.

Data Availability
The datasets compiled and/or analyzed during the current study are available from the corresponding author on reasonable request.

Social Network Analysis
Approximately 25% of all possible grooming associations between individuals occurred during the baseline and following the first intervention in both whole group and female-only networks (Table I). After the second intervention, density doubled in whole group and female-only grooming networks to ca. 54%, highlighting that new grooming relationships developed (Fig. 1, Table I). Following an initial increase in diameter following the first intervention, the diameter of whole group networks was lowest after the introduction of the novel male (post-intervention 2; Table I). Network diameter in female-only networks consistently decreased following each management intervention (Table I). Associations across networks and observation periods did not assort by degree (Table I), such that all individual were likely to associate, regardless of their degree.
Welch's t-tests indicate that no network values for removed individuals differed significantly from remaining individuals, suggesting that removed individuals likely did not play a disproportionate role in structuring the social network (Table SIII). However, both graph strength and the betweenness centrality of the novel male were larger than the group mean and were in the upper 30th percentile for all values in this period (Table SIII).
Graph strength and binary degree were significantly higher in empirical than in permuted networks in all observation periods (Table SIV). Harmonic centrality did not differ significantly from random in the majority of observation periods, although it was significantly higher in whole group networks during the baseline, post-intervention 1, and terminal observation periods than in random networks. In both baseline networks, betweenness centrality was significantly higher than in random networks (Table SIV).
The QAP analysis indicated that grooming associations in all observation periods correlated significantly with matrilineal relatedness, particularly in female-only networks. Following the first intervention, the correlation of grooming associations with matrilineal relatedness increased in whole group networks from 0.14 (P = 0.04) to 0.21 (P = 0.04), but decreased in female-only networks from 0.36 (P = 0.02) to 0.32 (P = 0.02). The second intervention increased associations between maternal kin in whole group networks (C = 0.27, P < 0.01), but not in female-only networks (C = 0.27, P = 0.04), in which grooming associations were the least correlated with matrilineal relatedness.

Baseline -Post-intervention 1
The relative rank of individual network measures did not differ significantly from those in permuted networks, highlighting that the first intervention did not disproportionately affect remaining group members (Table II, Fig. 2). Binary degree and graph strength remained stable in both networks following the first intervention (Table III, Fig. 2). We observed inverse associations between networks for the harmonic centrality, which increased in female-only networks and decreased in whole group networks. Betweenness centrality did not increase in either network (Table III, Fig. 2). These changes did not differ significantly from permuted networks (Table III). Lastly, females maintained their grooming relationships following the first intervention, as association matrices were moderately but significantly correlated with the baseline matrices (C = 0.42, P < 0.01).  Bold P-values indicate factors that differed significantly from permuted networks using a two-tailed approach.

Post-intervention 1 -Post-intervention 2
Again, the relative rank of individual network metrics did not differ from those in permuted networks, indicating that all changes to network metrics following the second intervention were homogeneous across the group (Table II). In whole group networks, graph strength remained stable despite an increase in binary degree, while the harmonic and betweenness centrality increased. Again, these changes did not differ from permuted values (Table III, Fig. 2). In femaleonly networks, binary degree, graph strength, and betweenness centrality increased following the introduction of the novel male, while there was no change in harmonic centrality. No contrasts between post-intervention periods varied significantly from permuted networks (Table III). Associations in femaleonly networks were moderately but significantly correlated between post-intervention periods (C = 0.36, P < 0.01).

Baseline -Post-intervention 2
In whole group networks, graph strength increased from 0.03 to 0.05, although this change did not remain significant following Bonferroni correction (Table III, Fig. 2). In female-only networks, binary degree increased more than that expected from random networks (Table III), and the number of interaction partners in female networks was significantly higher following the introduction of the novel compared to the baseline network. Observation periods from left to right: baseline, period following the removal of seven adults (six males and one female), contraception of six of the nine remaining females, and the transfer of the remaining group to a novel enclosure (post-intervention 1), period following the introduction of a novel, intact male to the group (post-intervention 2), and a terminal period 10-15 weeks after the introduction of the novel male (terminal). Line thickness denotes edge weight. Blue squares, red circles, and yellow triangles depict males, females, and juveniles, respectively. Letters denote individual identity for males (M), females (F), and juveniles (J) in age order. An asterisk (*) denotes the male that was introduced in the second intervention.  Observation periods include a baseline period; a period following the removal of seven adults (six males and one female), contraception of six of the nine remaining females, and the transfer of the remaining group to a novel enclosure (post-intervention 1); and an observation period following the introduction of a novel, intact male to the group (Post-intervention 2). Bold P-values indicate contrasts that differed significantly from permuted networks following Bonferroni correction using a two-tailed approach.
Interestingly, female-only grooming networks correlated most with the baseline following the second intervention (C = 0.44, P < 0.01).

Terminal Period
Whole group networks during the terminal period more closely resembled baseline networks in terms of density and diameter, while female-only networks were denser and had a smaller diameter than in the baseline network (Table I). Neither female-only nor whole group networks showed any assortativity (Table I). Individual changes to network metrics were relatively homogeneous across the group, as the relative rank of individuals did not differ significantly than that expected by chance (Table II). Graph strength and binary degree did not differ from the baseline in terminal female-only networks (Table IV, Fig. 2). In whole group networks, binary degree and graph strength were significantly higher than in baseline networks, but lower than immediately following the male introduction (Table IV). In whole group networks, the harmonic centrality was lower than during the baseline or following the introduction of the novel male, while betweenness was equal to that following the introduction of the novel male; neither metric differed significantly from the baseline or from networks following the introduction of the novel male (Table IV, Fig. 2). In female-only networks, the harmonic centrality was higher and betweenness centrality lower during the Observation periods include a baseline period (B); a period following the removal of seven adults (six males and one female), contraception of six of the nine remaining females, and the transfer of the remaining group to a novel enclosure (Post-intervention 1; PI1); an observation period following the introduction of a novel, intact male to the group (Post-intervention 2; PI2); and a terminal period 10-15 weeks after the introduction of the novel male (terminal; T). Squares with dashed lines depict males, circles with full lines depict females, triangles with dotted lines depict juveniles, and red stars depict mean values.

Table III
Coefficients of modified Tukey tests and network permutations for models comparing social network metrics in whole group and female-only grooming social networks for Observation periods include a baseline period; a period following the removal of seven adults (six males and one female), contraception of six of the nine remaining females, and the transfer of the remaining group to a novel enclosure (Post-intervention 1); and an observation period following the introduction of a novel, intact male to the group (Post-intervention 2). Bold P-values indicate contrasts that differed significantly from permuted networks following Bonferroni correction using a two-tailed approach. We highlight values that were no longer significant following the Bonferroni correction with a superscript a ( a ). As there was minimal variation in the betweenness centrality of whole group networks ( b ), mixed models incorporating individual as a random factor were overfit. We therefore used a linear model to evaluate differences in betweenness in whole group networks terminal period than during any other observation period (Table IV, Fig. 2). Neither measure differed significantly from the baseline observation period or the period following the introduction of the novel male. During the terminal period, grooming associations in both networks correlated more strongly with matrilineal relatedness than in other periods (whole group C = 0.44, P < 0.001; female-only C = 0.47, P < 0.01). Grooming associations during the terminal period correlated most with the grooming associations following the introduction of the novel male (terminalmale introduction: C = 0.42, P < 0.01; terminalbaseline: C = 0.33, P = 0.03).
Taken together, these results suggest that networks during the terminal period regained some of their overarching baseline structure (density, diameter), while the strength and number of associations remained elevated in whole group networks.

Knockout Networks
The knockout network model did not tend to mirror empirical changes in networks between the baseline and the period following the first intervention, suggesting that observed changes in networks were the result of behavioral responses to social and environmental disturbance, rather than simply the removal of several individuals. Knockout networks mirrored the increase in graph strength following the first intervention, although this did not differ  Bold P-values indicate contrasts that differed significantly from permuted networks following Bonferroni correction using a two-tailed approach.
significantly from random networks, as well as the increase in diameter (Tables I and V, Fig. 3c), highlighting that grooming associations responded flexibly to the removal of several individuals, contraception, and an enclosure move. Associations in the knockout model were more strongly correlated with matrilineal relatedness (C = 0.25, P = 0.03) than the empirical network following the first intervention was (C = 0.21). The knockout model mirrored the direction but not the magnitude of change in several network metrics including betweenness centrality and the harmonic centrality (Table V, Fig. 3).

Discussion
Our results suggest that, overall, grooming networks are relatively robust to social and environmental disturbance in Sulawesi crested macaques, even if significant changes in the number of grooming partners and the strength of associations occur. A key to the stability of grooming networks may lie in the social structure of Sulawesi crested macaques; social associations are weakly constrained by kinship (Baker and Estep 1985;Duboscq et al. 2013Duboscq et al. , 2017, allowing individual associations to respond flexibly to disturbance. Whole group and female-only networks were robust to the removal of a third of the group, an enclosure move, and the application of hormonal contraception. These results support findings from female baboons and suggest that social networks may be reasonably robust to the removal of individuals . Social network responses to the removal of individuals should, however, depend on the role of the removed individuals in the group (Flack et al. 2006). That the networks in this study were robust to the removal of several castrated males and one adult female, yet plastic to the introduction of a single intact male, is likely influenced by the reproductive status of the males present in that Bold P-values indicate contrasts that differed significantly from permuted networks following Bonferroni correction using a two-tailed approach. We highlight values that were no longer significant from permuted networks following the Bonferroni correction with a superscript a ( a ).
observation period. Castration decreases testosterone production, social dominance, and aggression, particularly aggression directed toward females (Michael and Zumpe 1993;Richards et al. 2009). The introduced intact male is therefore more likely to be dominant and more embedded in the group than the removed males and should play a key role in maintaining group cohesion (Flack et al. 2005;McCowan et al. 2011). Fig. 3 (a) Social network diagrams and (b-e) predictors ± SE from linear mixed models comparing whole group network metrics for a baseline network; a network simulating the removal of individuals (knockout); and a network following the removal of several adults (six males and one female), female contraception, and an enclosure move (post-intervention 1) for the Sulawesi crested macaques housed at Chester Zoo in July-December 2015. In the social networks, blue squares, red circles, and yellow triangles depict males, females, and juveniles, respectively.
Large changes in female grooming networks following the introduction of a novel and intact male suggest that social networks may be less robust to the addition of individuals than to their removal. Immigrating individuals can disrupt existing social hierarchies and increase the risk of infanticide, and as such present considerable stress to residents owing to potential fitness costs (Alberts et al. 1992;Engh et al. 2006a, b;Parmigiani and Vom Saal 1994;Pereira and Weiss 1991;Pusey and Schroepfer-Walker 2013;Wittig et al. 2008). Here, we observed signs that grooming networks expanded following the introduction of the new male as the number of grooming partners increased and remained elevated during the terminal period. As grooming in Sulawesi crested macaques is associated with a decrease in behaviors normally associated with anxiety such as aggression and self-directed behaviors (Aureli and Yates 2010), increases in the number of grooming partners may help buffer the stress of novel social conditions (Brent et al. 2014;Young et al. 2014). However, as immigration events in wild Sulawesi crested macaques appear only mildly stressful (Marty et al. 2017), new social associations in the presence of an intact male may, alternatively, highlight the importance of reproductive males in structuring social networks.
Male-female associations facilitate social relationships among unrelated and previously unaffiliated females (Silk et al. 2017). Associations between non-kin can reduce levels of nepotism (Silk et al. 2017), augmenting social cohesion and tolerance in the group (McCowan et al. 2008;Tiddi et al. 2011). As such, the increased affiliation between female non-maternal kin following the introduction of the male may be an additional factor contributing to the expansion of the grooming networks. In contrast, decreased affiliation between maternal kin, coupled with an increase in the number of grooming partners, may suggest that associations between paternal kin increased following the introduction of the novel male. Associations between paternal kin are posited to underlie high degrees of social tolerance in macaques (Schülke and Ostner 2008), and many weak paternal kin-associations can potentially mask effects from few strong associations between maternal relatives (Duboscq et al. 2017). However, as paternity data were not available for the group, we could not assess how grooming associations correlated with paternal relatedness.
Differences in group composition and cohesion may result in intraspecific variation in patterns of social affiliation (Duboscq et al. 2013). However, despite differences in group composition, networks during the terminal period somewhat regained their baseline structure, suggesting that social network responses to social disturbance may be temporary responses to novel social conditions rather than a result of group composition. Although social associations between individuals change over time (Blonder et al. 2012), social disturbance can "reset" social associations to a previous baseline state (Formica et al. 2017). Ultimately, this short-term social plasticity is adaptive, allowing a social group to maintain functionality under varying internal and external pressures (Wolf and Linden 2012). While plasticity in networks is driven primarily by homogeneous changes in individual-level network metrics across the group , biologically, social plasticity is driven by plasticity in the neural networks, hormones, and genes involved in different social contexts (Teles et al. 2016;Wolf and Linden 2012).
Several types of intervention occurred simultaneously during the first management intervention; as a result, we are unable to fully disentangle the effects of each intervention type. The use of contraception may have affected grooming relationships. Female Sulawesi crested macaques have exaggerated sexual swellings that reliably indicate when they are most fertile (Higham et al. 2012); hormonal contraception in Sulawesi crested macaques would have severely reduced, if not eliminated, sexual swellings in six of the nine females in the observed group (Portugal and Asa 1995). Female macaques are groomed less frequently by other females when maximally swollen (Aujard et al. 1998;Topping and Melfi 2006); therefore eliminating sexual swellings in the majority of females could result in increased grooming effort. Although we observed that grooming efforts increased following the management interventions, contracepted females were neither the least central nor the least groomed individuals in either network, suggesting that increases in grooming are more likely to be attributable to social changes than to the effects of contraception.
Maintaining social stability has important fitness consequences, and our results present insight into how social groups change their structure in response to group composition among other variables. While social associations are non-random (Formica et al. 2012;Silk et al. 2009), behavior itself is innately plastic, and, as we illustrate, individual responses to social and environmental disturbance are flexible. Nevertheless, homogeneous changes in associations result in robustness in the group: while grooming networks are largely structurally robust to social and environmental disturbance, associations in the network are less so, and changes in individual association are plastic across individuals. Importantly for captive management, we illustrate that castration may limit the role of male macaques in structuring social associations. We also demonstrate the utility of social network analysis as a practical tool with which to monitor the social implications of management recommendations. This work contributes to literature highlighting the utility of social network analysis for understanding the impacts of social and environmental disturbance on stable social groups.
Acknowledgments VBC was supported by the Chester Zoo Conservation Scholar Scheme. JMDL was awarded a University of Manchester Faculty of Life Sciences Research Scholarship. SS is supported by a Royal Society University Research Fellowship, Royal Society grant number UF110641. We sincerely thank the Chester Zoo Primate and Applied Science teams for their support with this project, Damien Farine for his helpful discussions regarding data analysis, as well as the editors and reviewers of this and previous versions of this article.
Author Contributions VBC, KJ, SLW, and SS conceived the idea, VBC collected the data, VBC analyzed the data with support from JMDL and SS, VBC wrote the manuscript; all authors provided editorial advice.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.