Social responses to the natural loss of individuals in Barbary macaques

In recent years, there has been considerable interest in investigating how animal social structure is affected by the loss of individuals. This is often achieved using simulations that generate predictions regarding how the removal of ‘key’ individuals from a group affects network structure. However, little is known about the effects of such removals in wild and free-ranging populations, particularly the extent to which naturally occurring mortality events and the loss of a large proportion of individuals from a social group affects the overall structure of a social network. Here, we used data from a population of wild Barbary macaques (Macaca sylvanus) that was exposed to an exceptionally harsh winter, culminating in the death of 64% of the adults from two groups. We analysed how social interaction patterns among surviving individuals were affected by the natural loss of group members using social networks based on affiliative (i.e., grooming) and aggressive social interactions. We show that only the structure of the pre-decline grooming networks was conserved in the post-decline networks, suggesting that grooming, but not aggression networks are resilient against the loss of group members. Surviving group members were not significantly different from the non-survivors in terms of their affiliative and agonistic relationships, and did not form assorted communities in the pre-decline networks. Overall, our results suggest that in primates, patterns of affiliative interactions are more resilient to changes in group composition than aggressive interaction patterns, which tend to be used more flexibly in new conditions.


Introduction
Social relationships, described through the content, quality and patterning of social interactions among group members (Hinde 1976) can have a strong impact on individual fitness, especially among group living animals, such as primates (Ostner and Schuelke 2018). Usually, individuals more socially integrated into their group-either via strong or a broad network of social bonds-tend to experience improved physiological wellbeing and health (Smith and Christakis 2008;Brent et al. 2011;McFarland et al. 2015), higher rates of reproductive success (Cameron et al. 2009;McFarland et al. 2017;Page et al. 2017;Ostner and Schuelke 2018) and improved survival Lehmann et al 2016;Ostner and Schuelke 2018). Consequently, the loss of closely connected individuals can have a negative impact on individual health (Christakis and Fowler 2009;Holt-Lunstad et al. 2010), influencing, for example, cardiovascular and pulmonary functions (Uchino 2004). Such losses, however, can be mitigated by an increase in affiliative association with other individuals (Engh et al. 2006; 1 3 Firth et al. 2017;Hobbs and Burke 2017). Social relationships have historically been described using composite measures of dyadic social interaction (e.g., Hinde 1976; Silk et al. 2003). More recently, however, social network analysis (SNA) has allowed researchers to move beyond the level of dyads to look at the patterning of social interactions across the entire group of individuals, examining not only the distribution of social relationships in group living animals (Lehmann and Ross 2011;Stanton and Mann 2012;Sosa 2016) but also the mechanisms governing animal social structure (Barrett et al. 2012;Farine et al. 2012). Animal social networks are often depicted as static entities. In reality, however, they are dynamic and constantly 'evolving' (Farine 2018). Individuals respond flexibly not only to changes in group composition but also to changes in the environment, which in turn affects social network structure (e.g., Wey and Blumstein 2012). For example, seasonal changes in resource availability have been shown to significantly affect the structure of social networks in Cape Giraffes (Giraffa camelopardalis giraffe: Wolf et al. 2018), Tasmanian devils (Sarcophilus harrisii: Hamede et al. 2009), breeding degus (Octodon degus: Wey et al. 2013), Asian elephants (Elephas maximus: da Silva et al. 2011), chacma baboons (Papio ursinus: Henzi et al. 2009) and vervet monkeys (Chlorocebus pygerythrus: Henzi et al. 2017). Indeed, even to simply maintain the existing social network structure individuals would need to adapt their social behaviour when faced with changes in group size or composition. Furthermore, the degree to which social networks are flexible in response to sudden shifts in social settings has also been demonstrated in animals living in groups that are stable over long periods of times, such as primates and some birds. For instance, in chacma baboons, affiliative networks undergo structural changes during periods of social instability, such as changes in the alpha male  or the death of a group member , which can lead to a reduction in individuals' stress levels (Engh et al. 2006). In captive rooks (Corvus frugilegus) artificially introduced changes in group composition were found to lead to adjustments of individual social relationships that served to 'recover' the original network structure (Borgeaud et al. 2016). Similar adjustments and flexibility in network structure following the loss of an individual have also been observed in great tits (Parus major: Firth et al. 2017).
A number of studies (e.g., Godfrey et al. 2013;Goldenberg et al. 2016) have demonstrated that animal social networks are flexible (i.e., networks can be structurally reshaped to accommodate new social or ecological conditions) and resilient (i.e., network retains their structure following a major perturbation). Note, while a network may be resilient to change, this does not necessarily mean dyadic relationships remain the same. To maintain social network structure when group size changes, social interactions are likely to have to change as well (Flack et al. 2006; Barrett et al. 2012; see Shizuka and Johnson 2020 for overview). While overall network flexibility in response to external conditions has been previously documented (Godfrey et al. 2013;Goldenberg et al. 2016;Puehringer-Sturmayr et al. 2021) the response of individuals to the (permanent) loss of a large number of group members, remains poorly understood (Farine 2021). Since social network structure is based on individual interactions, it is expected that the extent of structural changes following the loss of individuals will depend on the characteristics of the lost individuals, such as their personality, dominance status, age, sex and social position. For example, the removal of individuals that occupy a particularly central position in a network can lead to the fragmentation of existing social networks (e.g., in captive chimpanzees, Pan troglotydes: Kanngiesser et al. 2011; killer whales, Orcinus orca: Williams and Lusseau 2006;Columbian ground squirrel, Spermophilus columbianos: Manno 2008). In addition, the removal of certain age or sex classes may disproportionally affect the existing network structure (Williams and Lusseau 2006), because the position of an individual in a network is often related to its age and sex as, for example, found in Barbary macaques (Macaca sylvanus: Sosa 2016), olive baboons (Papio anubis: Fedurek and Lehmann 2017) and African lions (Panthera leo: Abell et al. 2013). Thus, targeted removals of certain age-sex classes through, for example, culling or trophy hunting, may reflect individual differences in boldness or visibility, which may bias the impact that such removals have on a social network. However, exactly what these impacts are is difficult to study, as it is very often either unethical or methodologically difficult to experimentally carry out such removals (Shizuka and Johnson 2020). Because of this, previous studies assessing the effects of individual loss on network structure tend to be based on either simulated removals (Raborn et al. 2002;Lusseau 2003;Williams and Lusseau 2006;Manno 2008;Kanngiesser et al. 2011;Levé et al. 2016;Fedurek and Lehmann 2017;Mourier et al. 2017;) or experimental short-term removal of individuals (Flack et al. 2006;Beisner et al. 2015;but see Cowl et al. 2020), with only a few studies investigating such effects in wild populations (Barrett et al. 2012;Franz et al. 2015). Studies using simulated removals usually compare the pre-removal structure of the network with that of the post-removal network (Shizuka and Johnson 2020). This allows researchers to assess to what extent targeted removals affect the network structure. However, simulated removals are unable to account for social adjustments and network rewiring (Wey et al. 2008;Blonder et al. 2012;Farine 2021) or for indirect effects deaths or dispersal of individuals may have on a network, such as changes in competitive regimes (Blonder et al. 2012;Shizuka and Johnson 2020). Neither simulations 1 3 nor short-term experimental approaches can reliably predict how social networks respond to the loss of individuals in the long-term (Barrett et al. 2012;Franz et al. 2015;Annagiri et al. 2017). In addition, little is known about the effects of removing a larger or random proportion of individuals (Naug 2009;Firth et al. 2017), especially in wild populations (Carter et al. 2009;Elliser and Herzing 2011), even though naturally occurring fluctuations in wild animal populations, due to the death or emigration, are not uncommon (Boesch and Boesch-Achermann 2000;Lehmann and Boesch 2004;Baigger et al. 2013;McFarland and Majolo 2013a). A recent study of Bechstein's bats (Myotis bechsteinii) showed that after a significant population decline, social relationships observed during the pre-decline period were still detectable for up to 1 year after the decline (Baigger et al. 2013), suggesting that although network rewiring occurred, there was a high resilience of bat networks to disturbance. Similarly, social network structure in Indian queenless ants (Diacamma indicum) appeared remarkably resilient to the experimental removal of leaders (Annagiri et al. 2017), and in African elephants (Loxodonta africanus) the hierarchical topology of social networks stayed stable over time despite changes in group composition due to, for example, poaching (Goldenberg et al. 2016). These studies suggest that, when possible (i.e., dyads are still present), social structure in many species are preserved following social and environmental disturbances. A recent study on the Cayo Santiago populations of rhesus monkeys that experienced a natural disaster which destroyed most of the island's vegetation, found that individuals invested more time building new relationships rather than strengthening existing ones (Testard et al. 2021). However, such studies on real-world network flexibility and resilience are still relatively rare (Elliser and Herzing 2011;Shizuka and Johnson 2020).
Here, we used data from two groups of wild Barbary macaques that were exposed to an exceptionally harsh winter that culminated in the death of 64% of the adults (30 of 47 individuals) across the two study groups  to assess how the loss of a large proportion of the group affects social network structure. Barbary macaques live in stable social groups with highly differentiated social relationships, as reflected by their social interaction patterns (McFarland and Majolo 2011) and are thus useful model candidates for investigating the effect of disturbances on network structure. Taking advantage of a rare natural experiment, our over-arching objective was to assess how resilient Barbary macaque social networks are following the natural removal of a large proportion of animals from the population. We compared macaque social networks before and after the significant population decline based on two different behaviours, grooming and aggression. We examined to what extent the macaques' pre-decline social network structure is preserved, and what network properties are more prone to change in response to the major shift in size and composition of the group. While the value attributed to affiliative relationships is largely considered to lie in the maintenance of a social bond over long periods (McDonald 2007;Cameron et al. 2009;Ostner and Schuelke 2018), agonistic relationships are often less stable and more responsive to changes in competitive regime (Sterck et al. 1997) or changes in social hierarchy (Schino et al. 1990). We, therefore, hypothesize that the structure of affiliative networks would more likely be preserved following perturbations than that of the aggression networks. To assess the extent to which network changes affected individual relationships, we also examined whether affiliative and aggressive relationships between survivors differed from those of nonsurvivors in the pre-decline period and whether these potential differences were preserved in the post-decline networks.

Data collection
We collected data on two groups (groups F and L) of wild Barbary macaques living in the Middle Atlas Mountains of Morocco. These groups were fully habituated and fed on a completely natural diet ). An adult female from group F died at the beginning of the study and was, therefore, excluded from the current analyses. Thirty of our 47 study animals died or disappeared during the exceptionally cold and snowy winter between December 2008and January 2009. We examined the social networks of the two groups both pre-and post-population decline (group F: predecline = Oct-Jan 2008, post-decline = Feb-Jun 2009; group L: pre-decline = Oct-Jan 2008, post-decline = Mar-Jun 2009). During the pre-decline period, groups F and L comprised 18 (7 females, 11 males) and 29 (10 females, 19 males) adult individuals (> 4 years), respectively. In the post-decline period, groups F and L comprised of six (4 females, 2 males) and six (3 females, 3 males) adult individuals, respectively.
Data were collected using continuous 20-min focal observations (Altmann 1974). The order of focal observations was randomized, and each subject was only sampled once per day. During focal observations, data were collected on all occurrences of grooming and aggressive behaviors (i.e., bite, charge, chase, displace, grab, lunge, or slap) exchanged between our focal animal and all other group members. During the pre-decline data collection period, 664 h of focal data were collected from our study animals (mean = 14.1 ± 9.1 SD h/subject). During the post-decline data collection period, 231 h of focal data were collected from our study animals (mean = 19.3 ± 4.9 SD h/subject). It was not possible 1 3 to record data blind, because our study involved focal animals in the field.

Social network analyses
We created social networks containing all group members over 4 years of age independently for grooming and aggressive interactions. Both types of interactions have previously been shown to be important predictors of survival in the current study population (Lehmann et al. 2016). For each group, we created separate social networks for the observation period before and after the population decline. In the grooming network, ties represent the rate of grooming exchanged within a given dyad (i.e., minutes per hour). In the aggression network, ties represent the rate of aggression exchanged within a given dyad (i.e., frequency of interaction per hour). Only focal sessions that were conducted when both members of the dyad were present were included in the calculation. All networks were directed, meaning that both in-coming (i.e., receiving grooming or aggression) and outgoing (i.e., giving grooming or aggression) ties were considered while constructing the networks. We used five network metrics: degree, strength and clustering coefficient to assess individual network integration and density, and average path length as measures of overall network connectivity and cohesion. These network metrics have been commonly used to describe the overall structure of a network and have been previously used in simulated knockout studies (Croft et al. 2008;Manno 2008;Fedurek and Lehmann 2017). Degree (derived from asymmetric matrices) indicates an individual's number of grooming or aggression partners, while strength measures the overall weight of all grooming or aggression ties (i.e., the total amount of time spent in grooming or aggressive interactions). Both variables have been shown to affect health and reproductive success in primates (Silk et al. 2003(Silk et al. , 2009Brent et al. 2011). It is worth noting that when looking at mean values of 'in' degree or strength, these have the same values as the corresponding means of 'out' degree or strength and thus we only present one value for each metric. Individual (local) clustering coefficient measures the degree to which the social partners of an individual are directly connected themselves and is calculated as the proportion of the existing ties to all possible ties between the individual's partners (Watts and Strogatz 1998; Barrat et al. 2004;Opsahl and Panzarasa 2009). In other words, the clustering coefficient shows the extent to which social interactions occur primarily within triads. Clustering coefficient in aggression networks has been shown to be a predictor of survival in this population (Lehmann et al. 2016). The calculation of the clustering coefficient was based on arithmetic mean of triplets (Opsahl 2009), while isolated nodes and nodes with only one connection were excluded from the analyses as the resulting clustering values of such nodes are, by definition, undefined (Opsahl 2009). Because discarding disconnected nodes from our sub-sampled networks might inflate clustering coefficient values, we repeated the analysis assigning '0' values to isolated individuals (results are in supplementary material, Table S1). Density (i.e., the extent to which a network is saturated), indicates the proportion of existing ties to all possible ties in the network for the predecline and post-decline networks. Density is an important variable to control for in analyses like this, as it has been shown to affect other network metrics, such as clustering coefficient (Croft et al. 2008). Average path length is the average number of steps along the shortest paths for all possible pairs of network's nodes (Barabási and Albert 1999). Average path length was calculated based on the lengths of the existing paths, discarding disconnected nodes (West 2001). Average path length has been argued to be an important indicative of parasite transmission dynamics in group living animals (Grear et al. 2013).
All measures were calculated based on weighted data except density, degree and average path length, which are, by default, calculated from binary measures. We used Kendal correlation to assess if our individual based network metrics (degree, strength and clustering coefficient) were significantly correlated with each other within networks.

Sub-sampling procedures
To assess to what extent the smaller post-decline networks differed from pre-decline networks, we used a sub-sampling technique (following Fedurek and Lehmann 2017), whereby we pseudo-randomly (i.e., conserving group structure as explained below) removed individuals from the pre-decline grooming and aggression networks and then calculated the network metrics described above. This is needed to differentiate the effects of a reduction in group size per se (simulated by the sub-sampling) from those of actual changes in network rewiring. Although it has been shown that survival is not random and at least partly linked to social interactions Lehmann et al. 2016), we have chosen to create null models based on random removals for two reasons: first, the links between survival and social structure have not been sufficiently studied and we do not yet know how universal these effects are; second, as other modelling studies have used random null models (Manno 2008;Fedurek and Lehmann 2017;Mourier et al. 2017) we wanted to maintain comparability when assessing how social interactions are re-shaped after the loss of group members. If observed post-decline networks resemble those obtained from sub-sampled networks, this would suggest high preservation of pre-decline social structure, while significant deviations from sub-sampled networks would suggest a significant level of rewiring.
Sub-sampling (n = 1000) was performed in such a way that the number and the sex ratio of individuals of the subsampled networks matched the natural 'removal' condition (note, that 1000 sub-sampled networks are well below the possible number of unique combinations). After each subsampling the mean of the individual network metrics of interest was calculated (following Fedurek and Lehmann 2017). The distribution of these mean network metrics for the sub-sampled and post-decline networks were then compared statistically (Bejder et al. 1998). The sub-sampling derived p values were calculated as the number of times the observed mean network metrics were higher than mean values for the same metric from sub-sampled networks divided by 1000 (i.e., the number of sub-sampling events). p values were considered significant when the observed network value was smaller (p < 0.025) or greater (p > 0.975) than 97.5% of the corresponding values resulting from the described above sub-sampling procedure (two-tailed test; following the procedure used by Bejder et al. 1998; see Farine 2017 for overview). For example, if the observed network metric value derived from the grooming or aggression network was greater than 97.5% of the corresponding values resulting from the sub-sampled network, the resulting p value would be p > 0.975 and would be considered significant.
Pre-decline networks were generally sparser compared to the corresponding post-decline networks (see Table 1), which might result in higher network metric values of the latter. To control for such possible effects of density, we ran an additional analysis whenever post-decline network values were significantly different from the sub-sampled pre-decline values. In this additional analysis, we filtered the post-decline networks by progressively increasing the cut-off threshold of the weakest connections until the density of these networks was within the same range as the sub-sampled pre-decline networks (i.e., density values differed by no more than 0.04). It has been argued that network filtering might potentially remove vital network information, as weak ties can play an important role in human and non-human primates social networks (Granovetter 1973;McFarland et al. 2017). Nevertheless, filtering procedures have been widely applied in SNA studies to mitigate noise resulting, for instance, from infrequent and random interactions, while preserving stronger and socially more relevant ties (Croft et al. 2008), especially in saturated networks (Fedurek et al. 2021). Network metrics resulting from filtered post-decline networks were bootstrapped (i.e., resampled with replacement, n = 100) and compared the resulting mean values with corresponding mean values (n = 100) resulting from the subsampled pre-decline networks using Wilcoxon sum rank test. This was done for strength and clustering coefficient but not degree as networks of the same size and density will always exhibit the same mean degree values.
The post-decline structure of social networks might have been, to some extent, determined by social interactions between survivors in the pre-decline period, if pre-decline grooming/aggression patterns observed between survivors were retained in the post-decline networks. To assess if this was indeed the case, we quantified the extent to which predecline social interaction patterns between survivors (n = 6) differed from those of randomly chosen group members Table 1 Observed and sub-sampled network metrics (with standard deviation) for the grooming relationships in two groups of Barbary macaque Observed metrics are based on pre-/post-decline observations, while sub-sampled metrics are based on pseudo-random removals of individuals, simulating the disappearance of the same proportion of randomly chosen group members. p values correspond to the comparison of observed post-decline and sub-sampled same-sized pre-decline network metrics and are considered significant when p < 0.025 or p > 0.975 (corresponding to a two-tailed test; see text for explanation). Significant results are indicated in bold. Italics indicate results close to significance. n = 6 'survivors' for each group using the same procedure as described above. That is, we created n = 1000 observed pre-decline networks, where the number of individuals and the sex ratio of the surviving individuals were preserved. The distribution of the network metrics derived from those networks composed of randomly selected group members (n = 6 individuals per network) was then compared against the observed measure of the pre-decline networks comprised of survivors only. p values were calculated as the number of times the observed mean value of the individual network metric in the pre-decline network between survivors was larger than the corresponding mean network metric value of randomly chosen individuals (i.e., survivors and non-survivors) from the pre-decline network divided by 1000 (i.e., the number of sub-sampling events; see Bejder et al. 1998;Whitehead 2008;Lusseau et al. 2009). Because we had no directional prediction regarding the social relationships amongst survivors compared to the rest of the group in the pre-decline networks, we used two-tailed statistics. Results were considered significant when the observed value was greater (p > 0.975) or smaller (p < 0.025) than 97.5% of the sub-sampled (average) network values (Bejder et al. 1998).
We calculated the assortativity coefficient (Newman 2003) for weighted networks, for our pre-decline networks, to determine whether survivors and non-survivors were assorted in separate communities in the pre-decline networks. The assortativity coefficient approaches − 1 in highly disassortative networks (e.g., survivors interact only with non-survivors), 0 when individuals interact randomly compared to an a priori specified assortment category, and 1 in highly assortative networks (survivors and non-survivors interact only among each other; Newman 2003). To calculate to what extent the assortativity coefficient of the pre-decline networks differs from what would be expected by chance, we resampled (n = 1000, without replacement) the assortment category (i.e., randomly assigning the survivor/non-survivor category to all group members) of the pre-decline networks and at each step calculated the assortativity coefficient. The distribution of the simulated values of the coefficient was then compared with the corresponding observed value. Subsampling-derived p values were calculated as the number of times the observed values of assortativity coefficient were higher or lower than the corresponding simulated values for the metric divided by 1000 (i.e., the number of simulations).
Finally, while our sub-sampling procedure assesses the preservation of the number of pre-decline social connections, they do not take into account the identity of the individuals. Thus, we used a Mantel test (Mantel 1967) to assess the preservation of individual dyadic relationships (i.e., the identity of grooming/aggression actors and recipients is taken into account). The Mantel test determines the extent to which the distribution of pre-decline relationships amongst the survivors are preserved in the post-decline networks i.e., by correlating pre-decline interaction matrices comprised of survivors only (each network contains n = 6 individuals, the pre-decline matrices of survivors were extracted from the full pre-decline matrices) with the corresponding postdecline matrices. p values were considered significant if p < 0.05.
Sub-sampling procedures and all statistical analyses were conducted in R (R Development Core Team 2018). We used the Mantel test for asymmetrical matrices (9999 iterations) using the 'ape' package (Paradis and Schliep 2019) in R. While conducting the Mantel test, the 'ape' package randomly permutates rows and columns of the second matrix to calculate z-statistics based on sum of rows multiplied by columns of lower triangle of the data matrices (Paradis and Schliep 2019). To calculate assortativity coefficient we used the 'assortment' package (Farine 2014) in R.

Results
Figures 1 and 2 illustrate the pre-and post-decline grooming and aggression networks of the two study groups.
Both, grooming and aggressive interactions were quite common in the observed groups, with grooming interactions being twice as common as aggressive interactions. Thus, grooming networks were generally denser than aggression networks (except post-decline networks in group L, Fig. 2, Table 1), although networks were relatively sparse in the pre-decline period and denser in the post-decline period. However, even though post-decline groups consisted of 6 individuals only, not all individuals interacted with every other individual (density < 1; Figs. 1 and 2, Table 1).
Some of our network metrics were significantly correlated with each other in some networks (supplementary material, Tables S2 and S3); however, no one metric was consistently correlated with another metric across all networks, and we, therefore, maintained all five network metrics in the analysis.

Comparison of pre-and post-decline network structures
Grooming. Although observed post-decline network metrics differed from pre-decline values, when controlling for differences in network size via sub-sampling, we found that despite the change in group composition, the general network structure of the pre-decline grooming network was preserved in the post-decline period: with the exception of density (both groups) and degree (group F), all other post-decline grooming network metrics were not significantly different from same-sized (sub-sampled) predecline network metrics (Table 1, Fig. 3). Because the change in density may have affected other network metrics, we controlled for this (see supplementary material) and found that the results remained the same even after controlling for differences in density (supplementary material, Table S4).
Therefore, our results suggest that the general structure of the grooming networks was mostly resilient to change (Table 1), although in some cases (group F) individuals on average increased the number of grooming partners they had (as indicated by mean grooming degree values).
Aggression. In aggression networks, the calculated values for pre-and post-decline network metrics differed (Table 2); however, contrary to grooming networks, for aggression networks all but one metric, average path length, also differed significantly between sub-sampled (i.e., controlling for network size) and observed post-decline networks (Table 2, Fig. 1 Grooming (a, c, e) and aggression (b, d, f) networks of Group F before the population decline (a-d) and after (e and f). Pre-decline grooming and aggression relationships between survivors are depicted separately in graph c and d, respectively. Grey nodes represent survivors, white nodes represent nonsurvivors. The thickness of a tie represents its relative strength compared to other ties within the same network, while the arrows represent the direction of a tie; rectangles and circles represent females and males, respectively Fig. 4). Results remained the same after controlling for density (supplementary material, Table S4).
This suggests that the aggression network structure is not preserved following a major population decline. Postdecline networks were denser, and individuals had (on average) significantly more aggression partners (higher mean degree), and more frequent (higher mean strength) and clustered (higher clustering coefficient) aggressive interactions in the post-decline period than expected from the subsampling of aggression networks. Fig. 2 Grooming (a, c, e) and aggression (b, d, f) networks of Group L before the population decline (b-d) and after (e-f). Pre-decline grooming and aggression relationships between survivors are depicted in graph c and d, respectively. Grey nodes represent survivors, white nodes represent nonsurvivors. The thickness of a tie represents its relative strength compared to other ties within the same network, while the arrows represent the direction of a tie; rectangles and circles represent females and males, respectively  Table 2 Observed and sub-sampled network metrics (with standard deviation) for the aggression relationships in two groups of Barbary macaque Observed metrics are based on pre-/post-decline observations, while sub-sampled metrics are based on pseudo-random removals of individuals, simulating the disappearance of the same proportion of randomly chosen group members. p values correspond to the comparison of observed post-decline and pre-decline sub-sampled same-size network metrics and are considered significant when p < 0.025 or p > 0.975 (corresponding to a two-tailed test; see text for explanation). Significant results are indicated in bold. n = 6 'survivors' for each group

Social interactions amongst survivors in the pre-decline period
When comparing the pre-decline interaction patterns between survivors to those of rest of the group, we found that there were no significant differences between grooming network metrics of survivors compared to randomly chosen group members (Table 3; Fig. 5). Similarly, pre-decline aggression networks did not significantly differ between survivors and randomly selected individuals although p values for degree and strength were close to significance (Table 4 and Fig. 6), tentatively suggesting that on average survivors might have had more aggression partners and a higher aggression frequency amongst themselves compared to the rest of the group. Furthermore, we did not find any evidence that individuals in the pre-decline grooming and aggression networks were assorted (using survival as category for assortativity), i.e., assortativity coefficients were not different from what would be expected in random networks (group F grooming: r = − 0.033, p = 0.599; group F aggression: r = − 0.028, p = 0.802; group L grooming: r = 0.149, p = 0.936; group L aggression: r = − 0.064, p = 0.354). This further strengthens our finding that survivors did not form separate communities in the pre-decline period.

Preservation of individual social interaction patterns following perturbation
Grooming and aggression interactions amongst survivors in the pre-decline networks significantly predicted interactions in the post-decline period in the larger social group L, but not in the smaller group F (Table 5).

Discussion
Comparing network matrices across networks of different sizes and densities is challenging, as many metrics are affected by network size and density. However, understanding how individual social interactions and thus network structure respond to (major) changes in group size is of crucial importance for our understanding of population dynamics and the importance of key-players within social groups. The aim of our study was to identify how a naturally occurring loss of a large number of individuals affects primate social networks and to assess if/how social interactions between survivors might have contributed to their communal survival. Overall, the results of our study were consistent across the two social groups, which differed in both size and demographic composition. We found that the properties of the grooming networks were generally better preserved in the post-decline period than those of the aggression networks. Even network metrics that are derived from both direct and indirect connections (such as clustering coefficient and average path length) showed a high level of resistance to change in the grooming network. However, in the aggression network all post-decline network metrics differed significantly from pre-decline values (with the exception of average path length). Interestingly, social interaction patterns between survivors did not differ significantly from those between all other individuals, suggesting that survival  was not linked to special social relationships amongst survivors, but rather a survivor's wider social integration across the group (McFarland and Majolo 2013). Individual social relationships existing before the decline were largely maintained in the larger social group but not in the smaller one, where pre-decline social interaction matrices were not significantly linked to interactions observed in the post-decline period.

Using a sub-sampling procedure to compare networks
Even though naturally occurring deaths are not likely to be entirely random, as certain characteristics of an individual may enhance their survival probability Lehmann et al. 2016), we used a random sub-sampling procedure in this study as our null model, because we feel that currently this is the only valid starting point available. The link between survival and individual sociality has not been sufficiently researched and is thus not fully understood, so that we cannot rely on a strong theoretical underpinning for assuming anything other than randomness for our subsampling procedure. Equally, although social animals are likely to readjust their social relationships after the death of conspecifics (e.g., Testard et al. 2021), when and how this happens is, as of yet poorly understood. Thus, our null model assuming randomness is currently the most parsimonious approach. Our conclusions, therefore, need to be understood as a deviation from the assumption of randomness, which, surprisingly, was not the case for all metrics investigated. Any changes in network structure described are changes with respect to values derived from sub-sampled networks of the same size rather than changes resulting from the comparison between network measures of the full pre-and postdecline networks (thereby avoiding changes simply due to differences in size). Choosing an appropriate time frame can be a problem in network studies based on observational data and may introduce sampling errors (Hoppitt and Farine 2018). In our study we have largely similar observation times for pre-and post-decline periods; in addition, individual observation times across both periods were very comparable. Although it could be argued that larger pre-decline networks might require a longer data collection period to capture all existing social relationships, our focal animal sampling protocol ensures that all individuals and their interactions had equal opportunity to be recorded, despite network size differences between time periods.
Differences in overall densities are an ongoing challenge for network comparisons. Here, we used a filtering procedure to control for this, whereby we progressively removed edges from the post-decline networks, starting with the weakest, until the densities of the two networks we compared were virtually equal (supplementary material 1). We found that results after controlling for density this way did not change our initial results (Table S4). However, filtering networks by edge weight could potentially distort the overall network structure. In our study we were particularly interested in strong enduring bonds, which are more likely to reflect actual dyadic preferences, thereby shaping network structure. Given this scope, we chose to retain stronger bonds instead of removing bonds randomly in the networks (which would have been an alternative procedure) as this ensured that the filtered networks are more representative of the strong social interactions in the real/true networks. In addition, although weak links have been shown to play an important role in primate networks , they are also more prone to observational error, bias and uncertainty (Farine and Strandburg-Peshkin 2015). Thus, the filtering procedure used allowed us to not only control for density differences, but also effectively compare non-random pre-and post-decline network structures.
Finally, using a sub-sampling procedure may mean that individuals could become (artificially) disconnected from their social partners (as no new ties are introduced), thereby introducing a higher than realistic number of zeros (totally disconnected individuals), which may bias the overall group mean. However, as can be seen in Figs. 3 and 4, we do not appear to have a large number of completely disconnected individuals in our sub-sampled networks and group averages do not appear to be overly biased. Individuals are expected to adjust their social relationships after a loss (Farine 2021) but currently the process of this re-wiring is not yet fully understood, so we could not include this in our sub-sampling procedure. However, our study provides a first insight into the flexibility of Barbary macaque social interaction patterns following a large population decline. Despite the difficulties comparing networks and deriving null models for studies like ours, one of the most interesting findings of our study is the difference in how affiliative and aggressive relationships change in response to perturbation.

Social interactions amongst survivors preand post-decline
Overall, the patterning of social interactions among the surviving individuals did not differ significantly from random subsets of six individuals in the group, nor did the assortativity coefficient indicate an assortment into different communities. Thus, their survival as a subgroup was most likely not driven by their social relationships amongst each other but rather by their own individual, social and physical characteristics (e.g., age or personality). There was, however, a tendency towards significance (i.e., p > 0.95; pre-decline degree and strength values between survivors were higher than 95% of corresponding values generated by pseudo-random subsampling; see Table 4 and Fig. 6) for higher aggression degree and strength amongst survivors than between random group members. This may suggest that the survival of individuals as a subgroup might have been influenced not only by their own social and physical characteristics but also, to some extent, by their specific social relationships within the survival cohort. This, to some extent, is in line with the findings of Lehmann et al. (2016) that the position in the agonistic network can be important for survival in Barbary macaques. It is unlikely that this higher frequency of aggression interactions amongst survivors observed in the pre-decline period reflects rank or dominance status, as the rank in this population was not found to be predictive of survival (i.e., pre-decline dominance status was not a significant factor influencing survival; Lehmann et al. 2016), rather, it might reflect an individual's personality or way to cope with pressure and unpredictability, leading to an enhanced likelihood of survival.

Social interactions and network structure after a major decline in group size
When faced with a dramatic change in group size or composition, individuals can either maximise the preservation of previously existing relationships (thereby retaining the pre-decline social network structure between survivors in the post-decline networks) or they can respond by altering their social relationships to adapt to the new situation (thereby departing from the pre-decline network structure between survivors). We found that Barbary macaques did both, depending on the behaviours used to build the social networks.
When looking at affiliative social interactions, such as grooming, the observed post-decline social network structure did not differ significantly from what would be expected based on sub-sampling, suggesting that the structure (but not necessarily the actual relationships) of the network was largely preserved following the loss of a relatively large proportion (64%) of individuals. On the other hand, when looking at the network based on antagonistic interactions, we found that the number of aggression partners, the frequency of aggressive interactions and clustering increased following the decline in group size. Although such an increase in network metrics could be due to the fact that in smaller groups individuals have more opportunities to interact with all available partners, such a theoretical explanation should hold for both affiliative and aggressive networks, and, therefore, does not explain the differences observed here. So, why are aggression networks more prone to structural changes following perturbations? Aggression in Barbary macaques is often used in the context of the acquisition of food, mates or rank (Kuester et al. 1998). It is unlikely that competition over food would have caused this increased level of aggression ; although the post-decline data were mostly collected during a period associated with lower food availability compared to the pre-decline period, the significant reduction in group size should have led to a decrease in food competition; in addition, Barbary macaques shift their diet preferences and foraging strategies during the winter from clumped (and thus potentially highly contested) fruits and seeds to highly dispersed (i.e., nondefendable) fir or cedar foliage rather (Mehlman 1988), thereby further reducing the level of food competition.
Equally, it is unlikely that the observed differences in aggression network structure are due to competition over mating opportunities, as the mating season of Barbary macaques (Oct-Dec) coincided with the predecline period of data collection, not the post-decline period. Therefore, although higher rates of aggression are typically observed in the mating season this would not explain the high rates of aggression observed in the post-decline period (i.e., encapsulating the non-mating season). However, a loss of a large proportion of group members is likely to result in high social uncertainty regarding dominance relationships. Social perturbation due to changes in the dominance status often results in an increase of aggression in primates (Aureli et al. 1992;Sapolsky 1992) a finding reported for both males and females (Schino et al. 1990). It is thus likely that following the loss of individuals hierarchies needed to be re-established, which might explain the observed increase in the rate, degree and clustering of aggressive interactions following the population decline. In addition, while individuals have a certain level of control over their grooming choices, aggressive interactions are less easily controlled, especially during times of rank uncertainty.
Interestingly, the structure of the pre-decline grooming networks was mostly preserved in the post-decline groups despite differences in density. Three of the five network measures used in this study were not statistically different from the corresponding measures obtained from our sub-sampling procedure. A similar resilience of grooming networks to experimental removals of individuals has been recently observed in Sulawesi crested macaque (Macaca nigra; Cowl et al. 2020). Although it has been shown that in many species, including primates, grooming and affiliative interactions can be used flexibly to create new social bonds (Silk et al. 2006), appease others (McFarland and Majolo 2011), enhance social buffering in times of social uncertainty (Engh et al. 2006) and as a currency in biological market situations (Noë and Hammerstein 1994;Barrett et al. 1999;McFarland and Majolo 2011), our results suggest that the grooming network structure per se, or rather, the distribution of grooming relationships on a group level, is relatively well preserved in Barbary macaques, even when group size and composition changes. This suggests that maintaining a certain level of social integration is more important than maintaining particular individual social relationships for Barbary macaques (see also McFarland and Majolo 2013), especially in winter, when social thermoregulation partner choice can be influenced by grooming networks (Campbell et al. 2018). Thus, the preservation of the grooming relationship structure (and the partner identity, at least in one group-see below) is an important asset of an individual's social capital that serves long-term functions and persists even in groups of reduced size.

Preservation of individual relationships
Although our sub-sampling procedure has shown that the strength of social interactions and the number of social partners was (mostly) preserved in the grooming but not aggression networks, the distribution of these social interactions, i.e., the identity of recipients and donors/actors of grooming and aggression interactions, was preserved only in the larger of the two groups of macaques. It has been shown in a number of species that individuals will reshape and even expand their social relationships after a loss of social partners (Engh et al. 2006;Franz et al. 2015;Firth et al. 2017;Shizuka and Johnson 2020), but little is known about the extent and the flexibility of this process, especially following the loss of a large proportion of group members. Long-term preservation of affiliative relationships following serious group perturbations has been shown in Bechstein's bats (Baigger et al. 2013) and captive rooks (Borgeaud et al. 2016), while Testard et al. (2021 have recently reported that rhesus macaques (Macaca mulatta) broadened their social networks by building new relationships rather than strengthening the previous ones following a natural disaster. In our study we saw the preservation of pre-decline relationships only in the originally larger group L, which has experienced a higher loss of individuals compared to group F. Why relationships were preserved in one group but not the other is currently not fully understood, especially as post-decline group size is identical and group composition is similar. Both groups live in ecologically very similar, and in part, overlapping habitats, thus it is most likely that social factors (e.g., kin structure, nature of social relationships) explain this difference. However, our sample size of only two groups is too small to even speculate, as the results may also be based on limited statistical power in case of the smaller group.

Conclusions
Overall, the results of our study show that the structure of the grooming but not the aggression networks was preserved following the loss of a large proportion of individuals, and we provide only partial evidence that the social relationships within those networks were preserved. We also found that in aggression networks the direct measures of social relationships, such as degree and strength, are more prone to change compared to some emergent properties resulting from such direct social connections, such as average path length. As such, the results of our study show that the network response to a loss of a large proportion of individuals has a multidimensional nature and cannot be easily predicted by applying only one method or looking at only one behavioural dimension. The results of our study may have important implications regarding the consequences of removal or loss of individuals from populations via translocation, or death through disease, poaching or hunting. Although behavioural plasticity is a very important characteristic of group living, enabling individuals to swiftly adjust to rapid changes in physical and social environment (DeWitt et al. 1998;Engh et al. 2006;Wright et al. 2010;Kappeler et al. 2013;Strier 2015;Firth et al. 2017), our study also suggests that the extent of such adjustment might be constraint by social profiles of individuals and might be limited by long-term benefits network stability may provide, such as predictability and partner availability for members of social groups.

Appendix: Illustration of affiliative and aggressive behaviours
The following images (Figs. A1, A2 and A3) illustrate affiliative and aggressive behaviours analysed in our study.

Consent for publication Not applicable.
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:// creat iveco mmons. org/ licen ses/ by/4. 0/.