Kin bias and male pair-bond status shape male-male relationships in a multilevel primate society

Male-male social relationships in group-living mammals vary from fierce competition to the formation of opportunistic coalitions or the development of long-lasting bonds. We investigated male-male relationships in Guinea baboons (Papio papio), a species characterized by male-male tolerance and affiliation. Guinea baboons live in a multi-level society, with units of one reproductively active “primary” male, 1–6 females, and offspring at the core level. Together with “bachelor” males, several units form a party, and 2–3 parties constitute a gang. We aimed to clarify to which degree male relationship patterns varied with relatedness and pair-bond status, i.e., whether males had primary or bachelor status. Data were collected from 24 males in two parties of Guinea baboons near Simenti in the Niokolo-Koba National Park in Senegal. Males maintained differentiated and equitable affiliative relationships (“strong bonds”) with other males that were stable over a 4-year period, irrespective of their pair-bond status. Remarkably, most bachelor males maintained strong bonds with multiple primary males, indicating that bachelor males play an important role in the cohesion of the parties. A clear male dominance hierarchy could not be established due to the high degree of uncertainty in individual rank scores, yet bachelor males were more likely to be found at the low end of the dominance hierarchy. Average relatedness was significantly higher between strongly bonded males, suggesting that kin biases contribute to the social preferences of males. Long-term data will be needed to test how male bonds affect male tenure and ultimately reproductive success. Males living in social groups may employ different strategies to increase their reproductive success, from fierce fighting to opportunistic alliance formation or the development of long-term bonds. To shed light on the factors that shape male strategies, we investigated male-male social relationships in the multilevel society of Guinea baboons (Papio papio) where “primary” males are associated with a small number of females and their offspring in “units” while other males are “bachelors.” Strong bonds occurred among and between primary and bachelor males and strongly bonded males were, on average, more closely related. Bachelor males typically had multiple bond partners and thus play an important role in the fabric of Guinea baboon societies. Across primate species, neither dispersal patterns nor social organization clearly map onto the presence of strong bonds in males, suggesting multiple routes to the evolution of male bonds.


Introduction
Individuals living in social groups typically develop differentiated relationships. Recurring affiliation and mutual support can give rise to the development of "strong bonds," which have been defined as relationships that are stable over time, equitable, and stronger (defined by more frequent or longer lasting affiliative interactions) compared to other relationships within the same group (Silk 2002). Long-term investments in social relationships are believed to mitigate the costs of group-Communicated by K. Langergraber living (reviewed in Ostner and Schülke 2018;Thompson 2019) and have been shown to provide adaptive benefits, such as increased longevity and offspring survival for females (e.g., baboons, Papio spp.: Silk et al. 2003Silk et al. , 2009Silk et al. , 2010Archie et al. 2014;feral horses, Equus caballus: Cameron et al. 2009; bottlenose dolphins, Tursiops sp.: Frère et al. 2010;reviewed in Snyder-Mackler et al. 2020) and increased reproductive success in males (e.g., Assamese macaques, Macaca assamensis: Schülke et al. 2010; chimpanzees, Pan troglodytes: Gilby et al. 2013).
Due to male competition over un-shareable fertilization, it was initially argued that male bonds should occur only rarely and, if so, exclusively among kin (van Hooff and van Schaik 1994). Thus, strong bonds between males were thought to emerge due to indirect fitness benefits and be most likely in three scenarios: (1) male philopatry, (2) joint dispersal with kin (natal or secondary), or (3) preferential dispersal into groups with resident kin (van Hooff and van Schaik 1994;van Hooff 2000). However, bond formation among non-kin males has by now been reported for several mammalian species (e.g., African lions, Panthera leo: Grinnell et al. 1995;chimpanzees: Langergraber et al. 2007; bottlenose dolphins : Gerber et al. 2019;and Assamese macaques: De Moor et al. 2020). These findings indicate that such relationships can also emerge through direct mutualistic benefits, where all contributors gain more than if they acted alone (Clutton-Brock 2002Ostner and Schülke 2014).
While dispersal patterns have been recognized as factors in shaping male-male relationships, fewer studies have considered the social organization of the species. Multi-level societies have been receiving increasing attention because of the shift in the cost-benefit ratio of group living compared to unilevel groups (Grueter et al. 2020). A multi-level social organization provides the benefits of larger groupings, such as joint resource defense and minimizing predation risk, while potentially decreasing the costs of living in large groups, such as higher food competition, through breaking up into smaller subgroups and using different areas during foraging (Grueter et al. 2020). Multi-level societies occur across a broad range of taxa, including elephants, zebras, and cetaceans (Grueter et al. 2020). "Core units," also known as "harems" or "one-maleunits," typically comprise one reproductively active male and several females with their young. Several core units form higher levels of aggregations (Grueter et al. 2020). Males that are not reproductively active ("bachelor males") may either be integrated at different levels of the society or form separate allmale groups (Grueter and Zinner 2004).
Prominent examples for a multi-level social organization in the primate order are geladas (Theropithecus gelada), hamadryas baboons (Papio hamadryas), Guinea baboons (Papio papio), and golden snub-nosed monkeys (Rhinopithecus roxellana) (reviewed in Grueter et al. 2017). In geladas and hamadryas baboons, "leader" males benefit from "leader-follower" associations in the form of longer tenure, higher number of females, and more offspring born within the unit (Snyder-Mackler et al. 2012;Chowdhury et al. 2015). In golden snub-nosed monkeys, males of different core units engage in cooperative defense of females against "satellite" (bachelor) males to increase paternity certainty . Kin-based alliances of bachelor males, in turn, jointly attack unit males to obtain reproductive opportunities . Despite their similar social organization, hamadryas and Guinea baboons differ from geladas and golden snubnosed monkeys in terms of their dispersal patterns: Guinea and hamadryas baboons show female-biased dispersal (Swedell et al. 2011;Kopp et al. 2014Kopp et al. , 2015Städele et al. 2015), while geladas are female-philopatric (Snyder-Mackler et al. 2014) and golden snub-nosed monkeys display malebiased dispersal (Chang et al. 2014;Huang et al. 2017).
Here, we investigate male-male relationship patterns in Guinea baboons. Guinea baboons live in nested multilevel societies where "units," composed of one "primary" male, one to six associated females, and immatures, comprise the core level (Goffe et al. 2016). "Primary" males maintain largely exclusive sexual relationships with the associated females. Several units, together with "bachelor" males that are typically not reproductively active, make up a "party." Two to three parties regularly aggregate into "gangs" and these gangs have overlapping home ranges (Patzelt et al. 2014). Genetic analyses provided evidence for female-biased dispersal (Kopp et al. 2015) with females transferring between units, parties, and gangs (Goffe et al. 2016), while males are predominantly philopatric. The social organization of Guinea baboons is similar to that of hamadryas baboons; in terms of their multi-level social organization, both species show intriguing parallels with that of early human societies (Fischer et al. 2017;Swedell and Plummer 2019).
Previous investigations of male-male relationships in Guinea baboons showed that adult males maintain high levels of spatial tolerance, engage in affiliative interactions, and support each other in coalitions (Patzelt et al. 2014;reviewed in Fischer et al. 2019). No correlation between the degree of affiliation and relatedness between adult males was observed, but the sample size in this initial study was small (Patzelt et al. 2014). More importantly, at the time of the initial study by Patzelt and colleagues, nothing was known about male-female association patterns and the existence of primary and bachelor males. We therefore specifically considered the distinct relationships primary and bachelor males may have, with the aim to clarify how their relationship patterns compared to those of hamadryas baboons and geladas. Our dataset comprised all large juvenile, adolescent, and adult males (N = 24) living in two parties.
To characterize the relationships among males, we first assessed whether male affiliative relationships fulfill the three criteria used for identifying "strong bonds" (Silk 2002), namely whether affiliative relationships among males of the same party are differentiated, equitable, and stable over time. We paid particular attention to the relationships among primary and bachelor males in order to clarify whether their relationship patterns conform to a "leader-follower" model or whether bachelor males rather formed sub-groups at the periphery of the parties.
Second, we re-evaluated the dominance hierarchy among males. In many species, rank relationships are central to predict interactions between subjects. In previous studies, we failed to identify a clear rank hierarchy between males (e.g., Kalbitzer et al. 2015). There could be two reasons for this absence of a linear rank hierarchy. One is that there is indeed no clear hierarchy; the other is that the sampling was insufficient to detect it. We therefore used a larger dataset and applied a novel method (Sánchez-Tójar et al. 2018) that allows for estimating the certainty with which a dominance hierarchy can be assessed.
Third, we assessed the average genetic relatedness of males with different bond strength and type. Kin selection theory (Hamilton 1964) predicts that strongly bonded males, including primary males and associated bachelor males, should be more closely related compared to other male-male dyads. Yet, while such patterns were reported in the similar societies of hamadryas baboons (Städele et al. 2016), previous investigations failed to find such a link in our study population. We here therefore set out to re-examine this question with a larger data set.

Field site and study subjects
Data for this study were collected over 19 months, from April 2014 to October 2015, at the Centre de Recherche de Primatologie (CRP) Simenti field station in the Niokolo-Koba National Park (PNNK), Senegal (described in Maciej et al. 2013). The Guinea baboon community in the area comprises over 400 individuals, including five habituated parties in two gangs, the Mare and the Simenti gang. We selected two parties with the highest numbers of large juvenile, adolescent, and adult males as our study groups (party 9 from the Mare gang and party 6 from the Simenti gang). Party was used as the level of analysis because males spend most of their time in spatial proximity, and show higher rates of affiliation, ritualized greeting behaviors, and coalition formation at the level of the party (Patzelt et al. 2014;Fischer et al. 2017;Dal Pesco and Fischer 2018). Party size and composition varied during the study period due to maturation, mortality, and immigration/emigration, with an average of 44.5 individuals in party 6 (range 2014-2015 = 40 to 48, average adult sex ratio: 0.89 males per female) and 44.5 individuals in party 9 (range 2014-2015 = 38 to 51, average adult sex ratio: 0.48 males per female; see supplementary Table S1.1a; all supplementary files are available in Online Resource 1). We performed continuous behavioral observations of all large juvenile males, small and large adolescent males, and adult males. Each month, two independent observers determined developmental stages and assessed age categories using physical markers (see supplementary Tables S1.2a, b). Males were introduced as focal subjects when they reached the status "large juvenile male." Three males disappeared during the study period (likely due to predation), while three males that transitioned to large juvenile status were included later in the study. We observed a total of 24 individuals (N = 14 in party 6, and N = 10 in party 9). The average observation time per subject was 17 months (range = 3 to 19). Additional data from two subsequent years (2016 and 2017) were added to investigate male relationship stability and equitability (see supplementary Tables S1.1a, b, S1.3 for detailed information about data collection between 2014 and 2017).

Data collection
During 2014 and 2015, we conducted behavioral observations during morning (6:30-12:30) and afternoon (15:00-18:00) sessions for a total of 410 observation days (884 contact hours for party 9 and 941 contact hours for party 6). All data were collected on Samsung Note 2 handhelds using electronic forms developed for our long-term data collection using the Pendragon 7.2 software (Pendragon Software Corporation, USA). It was not possible to record data blind because our study involved focal animals in the field. We recorded demographic changes (presence, birth, absence, or death), health status, and female reproductive state on a daily basis (Goffe et al. 2016). Aggression, displace/supplant, avoidance, unprovoked submission, coalitionary support, copulation and grooming were recorded ad libitum. We conducted focal follows (Altmann 1974) of 20 min that were balanced between subjects and between morning and afternoon sessions, with an average of seven follows per individual per month (2014-2015 total: 956 h and 2961 focal follows for both parties). The behavioral data collection protocol included recording continuous focal animal activity (moving, feeding, resting, and socializing) and all occurrences of social behaviors such as approaches (to within 1 m), grooming, contact-sitting, and ritualized greeting involving the focal animal (Dal Pesco and Fischer 2018). For all grooming and contact-sit bouts, durations were additionally recorded to the nearest second. Scan sampling (Altmann 1974) was used before and after each focal follow to record all neighbors of the focal male within 10 cm, 1 m, 5 m, and 10 m (total number of proximity scans 2014-2015=5911).

Data analyses and modeling
All data and statistical analyses (including figure preparation) were conducted using the R environment (version 3.4.4; R Development Core Team 2018) in the RStudio interface (version 1.1.456; RStudio Team 2018). We ran general and generalized linear mixed models using the R package "lme4" (version 1.1.14; Bates et al. 2014). Detailed information about sample size, data standardization/transformation, model structure, diagnostics (assumptions, collinearity, overdispersion, model stability), and effect sizes can be found in the supplementary section S2 and in each model table (supplementary  Tables S6.1a-S7.3b in supplementary sections S6 and S7). p values were obtained from the likelihood ratio test performed with the R function "drop1" with the argument "test" set to "chisq" (Barr et al. 2013). Effect sizes were calculated with the "r.squaredGLMM" function of the "MuMIn" R package (version 1.42.1; Burnham and Anderson 2002). Additional specific functions and packages are mentioned in each sub-section.

Genetic sampling, genotyping, and relatedness analysis
We collected fecal samples of all large juvenile, adolescent, and adult males (N = 24) belonging to our study parties and characterized genetic variation by analyzing 24 polymorphic autosomal microsatellite markers. Genetic sampling, storage, DNA-extraction, and genotyping methodologies are described in detail in Dal Pesco et al. (2020). When DNA extracts from tissue samples were available from previous studies (see Patzelt et al. 2014), at least one additional fecal sample was genotyped in order to cross-check individual identity. For all other individuals, analyses were conducted on three independent fecal samples to rule out identification errors during sample collection.
We calculated descriptive statistics for all 24 markers, estimated F IS , expected, and observed heterozygosity, and tested for Hardy-Weinberg equilibrium (HWE) and null alleles presence for all loci following the methodologies in Dal Pesco et al. (2020). All loci were polymorphic with allele numbers averaging 3.5 (range = 2 to 6). One locus (D1s548) showed signs of null alleles and significant deviations from HWE and was therefore excluded. Thus, a total of 23 loci were included in the following relatedness estimation. All details regarding the descriptive statistics of genetic markers can be found in the supplementary Table S3.1.
We estimated dyadic relatedness for all 24 males belonging to our study parties for a total of 134 within-party dyads using the R package "related" (Pew et al. 2015; also see Wang 2011; see methodological details in Dal Pesco et al. 2020). The Wang estimator (Wang 2002) showed the highest correlation coefficient between simulated relatives and genetic relatedness estimates (0.83), and was therefore selected as the relatedness estimator. Pedigree reconstruction analysis was not performed, as our genetic database is too recent to provide information on mother-offspring pairs for large juvenile, adolescent, and adult individuals.

Male-male social bonds: strength, equitability, and stability
We measured the strength of dyadic affiliative relationships within parties by computing the Dyadic Composite Sociality Index (hereafter DSI; Silk et al. 2006b). This index ranges from 0 to infinity and, for a given dyad, measures the deviation in affiliative behavior compared to all other dyads in the same party. The mean DSI value is 1, lower values represent affiliative relationships weaker than average, and higher values indicate stronger than average relationships. For each male-male dyad, we calculated the DSI on a yearly basis (January to December). The only exception was the analysis of average dyadic relatedness and male sociality, for which the DSI was computed for the whole study duration (i.e., 2014 and 2015 combined). The following behaviors were included in our index calculation (see supplementary Table S4.1): grooming frequency (Gf) and duration (Gd), contact-sit frequency (Cf) and duration (Cd), and frequency of approaches to within 1 m (Af). Note that only approaches that were not followed by social behavior (positive or negative) within 10 s were considered in the DSI calculation to avoid redundancies with the other behaviors. The DSI was calculated using the following formula: DSI = [(Gf ij /Gf meanX ) + (Gd ij /Gd meanX ) + (Cf ij /Cf meanX ) + (Cd ij /Cd meanX ) + (Af ij /Af meanX )]/5; where terms containing ij always represent the dyadic behavioral frequency/duration for a certain dyad ij and terms containing meanX represent the mean frequency/duration for a certain behavior in party X . One requirement of the DSI is that the components are correlating (Silk et al. 2006a, b). We therefore assessed the correlation between all behavioral components in our index performing a Kendall's tau correlation test using the "cor.test" function in the "stats" package (R Development Core Team 2018). All behavioral components included in the composite index were positively correlated (see supplementary Table S4.2). Based on the methodology by McFarland et al. (2017), the number of strong bonds per male was based on the number of higher-than-average DSI values.
We calculated male-male relationship equitability within parties by computing the Grooming Symmetry Index (GSI; Silk et al. 2006a) over a 4-year period on a yearly basis (January to December) for each male-male dyad. All males present for at least 1 year were included in this analysis (N = 24). No gaps occurred between observation periods of individual males. The GSI measures how evenly grooming is balanced within a dyad by considering the duration of grooming given and received. The GSI ranges from 0 to 1, where 1 represents a completely equitable relationship and 0 corresponds to a completely unequitable relationship where one individual received all the grooming and gave none. The GSI was calculated using the following formula: GSI = 1 − abs[(Gd i→j /Gd i↔j ) − (Gd j→i /Gd j↔i )]; where Gd indicates grooming duration, → indicates grooming directionality and ↔ represents the total grooming exchanged by a certain dyad ij (i.e., given and received). Only dyads that exchanged at least four grooming bouts per year were included in this calculation and the following analysis. To investigate the link between male-male relationship strength and equitability we tested whether grooming equitability (GSI) was significantly higher for dyads with greater dyadic relationship strength (DSI). We ran a general linear mixed model with a Gaussian error structure (Baayen 2008) with GSI as the response and DSI as the main predictor of interest. We included year and party membership as fixed control factors, and male identities and dyad identity as random intercepts. The following random slope components were included in the model: DSI within ID1, ID2, and dyad identity. To control for multiple membership of individuals (IDs) within dyads, we ran the model 10,000 times, randomly assigning dyad members to either ID1 or ID2.
We calculated male-male relationship stability using the Partner Stability Index over a 4-year period (PSI; Silk et al. 2006aSilk et al. , 2013. The PSI measures variation in individual partner preference based on each individual's top partners across several time periods. The PSI ranges from 0 to 1, where 1 is the highest stability value, for individuals with the same top partners over all periods, and 0 is the lowest stability value, for individuals that changed top partners in every period. All males present for at least two consecutive years (N = 23) were included in this analysis, and no gaps occurred between observation periods of individual males. Because for some males, no single third top partner could be identified, we calculated all PSI values using the top two partners. The PSI was calculated using the following formula: PSI = [(ns − u)/(ns − s)]; where n is the number of periods the individual was present in the study group, s is the number of top partners considered in the analysis (i.e., two), and u is the number of unique top partners the individual had across periods. To test if the observed preference patterns were different than those expected by chance, we compared the observed PSI values for each individual male against mean expected PSI values for random partner choice calculated from 10,000 permutations (Silk et al. 2012;Kalbitz et al. 2016) using an exact Wilcoxon signed-rank test ("wilcoxon.test" function in the "stats" package R Development Core Team 2018; two-sided, paired, confidence level = 0.95). Furthermore, we calculated the percentage of males with an observed PSI value greater than 95% of permuted values by computing the proportion of permuted values lower than the observed value for each male and by assessing how many males had a proportion greater than 95. To investigate the link between male-male relationship strength and stability (Kalbitz et al. 2016), we additionally tested if the PSI for each male's top two partners was different than the PSI for his third and fourth strongest partners. Only the males that had unique third and fourth partners for at least 2 years were included in this analysis (N = 21). For each individual male, the PSI calculated based on the top two partners was compared with the PSI calculated based on the third and fourth partners using an exact Wilcoxon signed-rank test ("wilcoxon.test" function in the "stats" package R Development Core Team 2018; two-sided, paired, confidence level = 0.95).

Male dominance hierarchy assessment
Ad libitum and focal data on male-male aggression (i.e., threats, lunges, chases, physical fights), displace/supplant, avoidance, and unprovoked submission were used to assess the dominance hierarchy among large juvenile, adolescent, and adult males of the same party. Decided dyadic interactions were used to compile a winner/loser matrix to determine dominance relationships for each studied party. We excluded any aggressive interaction that followed one or more polyadic interactions within the same aggressive event (polyadic event). When several interactions per dyad occurred within the same aggression event, only the highest intensity interaction was considered. We used the "aniDom" package (version 0.1.3; Farine and Sánchez-Tójar 2017) to calculate the randomized Elo-rating scores for males of each party (Sánchez-Tójar et al. 2018), as this method performed best for both intermediate and low hierarchy steepness in a data simulation study (Sánchez-Tójar et al. 2018). This method randomizes the order in which interactions occur to avoid temporal biases and produces several individual rank scores, allowing for the calculation of mean individual ranks and the respective 95% range of individual scores (Sánchez-Tójar et al. 2018). By providing information about the scores' distributions, this method allows the estimation of the individual ranks' uncertainty. The use of this methodological approach also allowed us to assess the hierarchy steepness and uncertainty independent of group size and sampling effort. We assessed the degree of hierarchy orderliness using the triangular transitivity index (Shizuka and McDonald 2012). Compared to others, this index is not influenced by dataset sparseness or variation in group size (Shizuka and McDonald 2012). For methodological details and functions used in this analysis, see supplementary section S5a.

Male status and relationships among primary and bachelor males
Data on female-male associations and male status (primary, bachelor) were recorded daily. Female unit transfers were recorded and verified using copulations, grooming bouts, contact-sit bouts, greetings, and aggression events from focal and ad libitum data (see Dal Pesco and Fischer 2018). We used male-male spatial proximity to identify association patterns among males (following Chowdhury et al. 2015), so that we could directly compare our results to those for hamadryas baboons (Chowdhury et al. 2015). We first used a change point analysis to detect discontinuities in the distribution of male-male approaches (within 1 m) including all large juvenile, adolescent, and adult males within-parties (R package "changepoint"; version 2.2.2, Killick et al. 2016). Males were classified as "associated males" when the approach rate was above the change point threshold (see details in supplementary section S6a and supplementary Fig. S6.1). For each primary male, we determined the number of associated bachelor males, i.e., associated males that were not primary males. To control for demographic changes and the variation in male status, we calculated the number of associated bachelor males per primary male as a yearly average weighted by the duration of the association in days.

The role of kinship in male sociality
We tested whether the dyadic relatedness was significantly higher on average for strongly bonded dyads versus nonstrongly bonded dyads, and for primary and their associated bachelor males versus all other dyads. Male-male dyads were characterized as strongly bonded if their DSI (calculated over both years, 2014 and 2015) was above 1 and as primary/ associated bachelor males if they fell into this category at least once during this period. We ran two general linear mixed models with a Gaussian error structure (Baayen 2008) with dyadic relatedness estimates (Wang estimator, Wang 2002) as the response and relationship type 1 (strongly bonded versus non-strongly bonded) and type 2 (primary and their associated bachelor males versus other male-male dyads) as the main predictor, respectively. In both models, we included party membership as a fixed control factor and male identities as random intercepts. To control for multiple membership of individuals (IDs) within dyads, we ran both models 10,000 times, randomly assigning dyad members to either ID1 or ID2. Note that in this analysis, the two study years were considered as a single study period because the relatedness estimate is constant over time.

Male-male social bonds: strength, equitability, and stability
Male-male affiliative relationships were differentiated, equitable, and stable over time. The distribution of the Dyadic Composite Sociality Index (DSI) was highly skewed, indicating that relationships were differentiated (see Fig. 1, supplementary Fig. S4.1). The DSI for the entire study period (see supplementary Tables S4.3, S4.4 for yearly values) ranged from 0.00 to 11.90 with a mean of 1 and a median of 0.10. Only 22.4% of dyads (30/134 dyads) had a DSI above group average, and the top 10% of relationships had a DSI above 3.27. To illustrate the variation in behaviors in relation to DSI, average approach rates per hour were 0.11 for a DSI of 0.1, 0.49 for a DSI of 1, and 1.42 for a DSI of 5. The average number of strong bonds per male (calculated as the number of higher-than-average DSI values) was 2.50 (SD = 1.93; range = 0 to 6): 3.29 for party 6 (SD = 2.13; range = 0 to 6; N = 14) and 1.40 for party 9 (SD = 0.84; range = 0 to 3; N = 10). Average DSI across all strongly bonded male dyads was 3.93, indicating that these dyads affiliated almost four times as often/long compared to the average of the party.
Across the 4 years, a total of 699 grooming bouts were observed between males belonging to the same party. Of a total of 133 male-male within party dyads included in this analysis, 48 dyads exchanged at least one grooming bout, indicating that males were selective in their choice of grooming partners. The yearly Grooming Symmetry Index (GSI) considering only dyads which exchanged at least four grooming bouts during a given year (N = 30) was 0.59 (SD = 0.27; range = 0.00 to 0.98). We found some (weak) evidence that dyads with higher DSI values showed significantly higher GSI values (mean lmer estimate ± SE = 0.08 ± 0.04, mean p = 0.050 [range: < 0.001 to 0.883] based on 10,000 permutations with 71% of models having a p value ≤ 0.050; see supplementary Tables S7.1a, b).
Partner choice was not random and relationships were stable over time. The average observed Partner Stability Index (PSI) was 0.72 (SD = 0.22; range = 0.25 to 1.00) with 69.6% of males having observed PSI values greater than 95% of the permuted values. The observed PSI values were significantly higher than mean permuted PSI values (Wilcoxon signed-rank test: V = 271, N = 23, p < 0.001; Fig. 2 These results indicate that our observed preference patterns were different from those expected by chance, i.e., by random partner choice (see Fig.  2a, b). Stronger bonds were also more stable: the PSI values based on the top two partners were higher compared to the PSI values based on the 3rd and 4th ranked partners (Wilcoxon signed-rank test: V = 85.5, N = 21, p = 0.040).

Male dominance hierarchy assessment
Out of 1026 within-party male-male aggressive interactions recorded during 2014-2015, only 19.1% (196 interactions) could be used to assess dominance because the others were non-decided (42.7% of 1026), non-dyadic (59.7% of 1026), or comprised repeated interactions within the same bout of aggression. These 196 interactions were used in combination with 209 displace/supplant, avoidance, and unprovoked submission interactions, for a total of 405 interactions, to determine the rank hierarchy. Both parties showed a hierarchy of intermediate/low steepness with high variation in randomized Elo-rating scores and great overlap between the 95% score range for most males (Fig. 3). Those few individuals whose scores did not overlap with those of the highest-ranking individuals were all large juvenile, adolescent, or old males for the entire duration of the study (see Fig.3; also see supplementary online material in Dal Pesco and Fischer 2018). Hierarchies showed high transitivity, indicating high levels of orderliness. Yet, due to the extremely low rate of aggression (Patzelt et al. 2014;Kalbitzer et al. 2015) and usable proportion of interactions, our attempts to assess dominance yielded uncertain estimates. Therefore, dominance could not be included as a predictor in our analyses (see supplementary section S5b, supplementary Fig. S5.1, supplementary Tables S5.1a, b for detailed explanations).

Relationships among primary and bachelor males
The change point analysis determined that a total of 34 malemale dyads should be classified as "associated male-male dyads." Associated dyads exchanged grooming and contactsit bouts at a mean (± SD) rate of 0.39 bouts/h (± 0.36); nonassociated dyads did so at a mean (± SD) rate of 0.01 (± 0.02) bouts/h (see supplementary section S6a). Associated males interacted at significantly higher rates with the unit females (e.g., interaction rate per hour of focal observation; associated males: mean ± SD = 0.21 ± 0.31; non-associated males: mean ± SD = 0.01 ± 0.05; see supplementary section S6b, supplementary Fig. S6.2, and modeling supplementary Tables S6.1a, b, c). While four of these male-male dyads included two males with bachelor-status for the entire study duration, the remaining 30 dyads included at least one male that had primary status for part or all of the study. As male status changed throughout the study period, dyad type varied accordingly. Nevertheless, no associated dyads were composed of two males that both had primary status for the entire duration of the study and only nine dyads were composed of two primary males for a limited duration (duration in days, mean ± SD: 96 ± 111, range = 2 to 283 on 579 total days). During the study period, 13 of the 17 primary males (76.5%) were associated with at least one bachelor male. More than half of these primary males (8 of 13 males; 61.5%) were associated with more than one bachelor male. The weighted average number of associated bachelor males per primary male was 1.65 (SD = 1 .47; range = 0.00 to 4.08): 2.56 (SD = 1.44; range = 0.00 to 4.08; N = 14) for party 6 and 0.62 (SD = 0.55; range = 0.00 to 1.43; N = 10) for party 9. Most bachelor males were classified as either large juvenile/adolescent or late prime/old (12 of 15 males: 80%). All bachelor males were associated with at least one primary male and interacted with that primary male's females (also see supplementary section S6). Remarkably, 66.7% of bachelor males (10 of 15 males) were associated with multiple primary males (average number of units 2.53 ± 1.36 SD; range = 1 to 5; median = 3.00. Strong affiliative relationships between males were not restricted to dyads consisting of one primary and one bachelor male, however. Keeping into account demographic and status changes, an average of 27.9% of strong bonds occurred between two primary males or two bachelor males (see Fig. 4). Bachelor males did not associate in all-male groups nor form a sub-group within the party. Instead, bachelor males were an integral part of the party and were very well socially integrated, maintaining strong bonds with one or more primary males as well as with other bachelor males (Fig. 4). We did not observe bachelor males that were not affiliated with any primary male.

The role of kinship in male sociality
Relatedness estimates for male-male within-party dyads averaged 0.06 (SD ± 0.24) and ranged from − 0.56 to 0.69 (median = 0.05). Strongly bonded males had significantly higher average relatedness estimates than non-strongly bonded males  Fig. 3 Randomized Elo-rating scores for male Guinea baboons of our two study parties (party 9, N = 10, and party 6, N = 14, in the upper and lower plot, respectively). Each point represents the average randomized Elo-rating score per male with the respective 95% score range (i.e., within 2.5% and 97.5% quantiles). Note that rank scores decrease from left to right, with the highest individual rank score at the upper left and the lowest individual rank score at the lower right. Males that were large juvenile/adolescent or old for the entire duration of the study (2014)(2015) are labeled as * and **, respectively (mean lmer estimate ± SE =− 0.21 ± 0.05, mean p < 0.001 [range: < 0.001 to 0.003] based on 10,000 permutations, Fig. 5a, see supplementary Tables S7.2a, b). Note however that the range of dyadic relatedness overlapped to a great extent (strongly bonded: range = − 0.31 to 0.69; non-strongly bonded: range = − 0.56 to 0.57). The same pattern was found for primary males and their associated bachelor males, where dyads composed of primary/associated bachelor males had significantly higher average relatedness than all other malemale dyads (mean lmer estimate ± SE = − 0.20 ± 0.05, mean p < 0.001 [range: < 0.001 to 0.003] based on 10,000 permutations, Fig. 5b, primary/associated bachelor males: range = − 0.22 to 0.69; all other male-male dyads: range = − 0.56 to 0.57; see supplementary Tables S7.3a, b).

Discussion
Guinea baboon males form highly differentiated and equitable affiliative relationships with partner preferences being stable over a 4-year period (and most likely longer). These findings justify the claim that male Guinea baboons form strong social bonds. Strong bonds have previously been reported between females in female-philopatric baboons (Silk et al., 2006a(Silk et al., , b, 2009) and between males in male-philopatric chimpanzees (Mitani 2009) and in dispersing male Assamese macaques (Kalbitz et al. 2016). Strong male-male bonds were found both in uni-level groups (e.g., Barbary macaques, Macaca sylvanus: Young et al. 2014; Assamese macaques: Kalbitz et al. 2016), as well as in multi-level societies (e.g., Guinea baboons: Patzelt et al. 2014;bottlenose dolphins: Connor and Krützen 2015;Gerber et al. 2019). Given the variety of social systems where strong bonds between males occur, neither social organization (multi-level or uni-level) nor dispersal patterns (female-or male-biased) are obvious predictors of the occurrence of male-male bonds.
In Guinea baboons, primary males maintained strong bonds both with other primary males and with bachelor males, and strong bonds also occurred between bachelor males, indicating that these bonds are not restricted to the unit level. With few exceptions, most bachelor males were either large juvenile/adolescents or late prime/old males, while males in their prime generally had primary status. In our population, similar to the findings for hamadryas baboons (reviewed in Grueter et al. 2017), we did not see any all-male groups. In contrast, in both snub-nosed monkeys and geladas, all-male groups-also called "bachelor groups"-can be observed (reviewed in Grueter et al. 2017). Furthermore, in our study population, we did not observe bachelor males that were not affiliated to any primary male (termed "solitary" in hamadryas  (2014 and 2015) for all male-male within-party dyads. Party 6 (N = 14) outlined in violet (upper left) and 9 (N = 10) outlined in green (lower right). For illustration purposes, different node colors reflect average male status: dark gray nodes for primary males (with primary status for ≥ 80% of study time); light gray nodes for transition males (with primary status for ≥ 20% and < 80% of study time); white nodes for bachelor males (with primary status for ≥ 0% and < 20% of study time) baboons). Instead, bachelor males could simultaneously be associated with multiple primary males. Our findings show that in Guinea baboons, bachelor males are an integral part of the group ("party") and can be very well socially integrated. This sets male Guinea baboons apart from male hamadryas baboons and geladas. Inter-unit interactions are infrequent in hamadryas baboons, with occasional affiliation between leader males and follower/solitary males, but with interactions between leaders limited to threats, avoidance, and stereotyped greeting behavior ("notifying";Kummer 1968;Schreier and Swedell 2009). In geladas, leaders usually ignore each other (Dunbar 1983) and male-male affiliation is restricted to allmale ("bachelor") groups (Pappano 2013) or rare interactions between leaders and their followers, which most often stay at the periphery of the unit (Dunbar and Dunbar 1975).
The occurrence of follower/bachelor males differs between the three species. Most Guinea baboon primary males had at least one associated bachelor male (76.5% versus 55.4% in hamadryas baboons and 33.3% in geladas; Chowdhury et al. 2015;Snyder-Mackler et al. 2012, respectively), and most of these males had more than one associated bachelor male (61.5% versus 38.9% in hamadryas baboons; mean number of bachelor/follower males per unit: 1.65 versus 0.80 in hamadryas baboons; see Chowdhury et al. 2015). Note that we differentiated between associated and non-associated males using proximity at a much closer distance (1 m) than that used in hamadryas baboon studies (5 m) (Chowdhury et al. 2015). This was necessary as Guinea baboons are highly gregarious, and preferred associations are difficult to discern using larger distances (Goffe et al. 2016). Contrary to hamadryas baboons, where follower males tend to associate with a single leader male and further associations are shortterm and usually do not extend to females (Chowdhury et al. 2015); in Guinea baboons, bachelor males have multiple associations that may last several years and involve social interactions with both primary males and their associated females.
In contrast to an earlier study with a smaller sample size (Patzelt et al. 2014), we found here that the average relatedness was significantly higher between strongly bonded males and between primary males and their associated bachelor males compared to all other dyads. These results suggest that male-male relatedness plays a role in shaping relationship patterns between males. Given the high mate fidelity of females once they joined a unit (Goffe et al. 2016), infants born into the same unit will most likely be paternal siblings. Thus, some of the early peers are related. Observations from captive Guinea baboons in the Chicago Zoo (Boese 1975) indicated that males establish bonds already as juveniles; whether this Primary maleassociated bachelor male dyads All others male-male dyads Strongly bonded dyads Non-strongly bonded dyads Dyadic estimated relatedness (Wang) Dyadic estimated relatedness (Wang) Fig. 5 Average dyadic male-male relatedness estimates (Wang 2002) by a relationship type 1 (strongly bonded dyads versus non-strongly bonded dyads); b relationship type 2 (primary and their associated bachelor males versus other male-male dyads). Points depict within-party dyads (134 dyads, 24 males). Boxplots depict the median, the lower, and upper quartiles (25% and 75%), and the range excluding outliers (vertical line). Short horizontal whiskers depict the bootstrapped 95% confidence intervals (party manually dummy coded and centered) pattern holds in the wild and whether males preferentially recruit partners from their natal units can only be clarified when long-term demographic data are available. While individuals should bias their behavior toward closely related partners due to indirect fitness benefits (Hamilton 1964), it is also possible that males establish bonds with non-kin, if such bonds provide sufficient direct benefits (Mitani et al. 2002;Krützen et al. 2003;Langergraber et al. 2007Langergraber et al. , 2009Gerber et al. 2019;Sandel et al. 2019;De Moor et al. 2020). Non-kin ties have been deemed crucial in vampire bats, where they can result in food-sharing network expansion and aid in dealing with the disappearance of key food-sharing partners (Desmodus rotundus: Carter and Wilkinson 2015; Carter et al. 2017). In the absence of pedigree information, we are not able to assign the males in our study to specific kin classes, nor are we able to identify male-male dyads that are certain to be unrelated. Therefore, it is not yet possible to decide to which degree kinship and/or familiarity via shared unit membership, for instance, play a role in bond formation.
Corroborating and expanding earlier findings (Patzelt et al. 2014;Kalbitzer et al. 2015), we found that male hierarchies in Guinea baboons were of intermediate/low steepness, individual rank scores were highly variable with great overlap between inter-individual score ranges and, most importantly, estimates remained uncertain. Combined with past unsuccessful attempts to establish a significant linear hierarchy (Patzelt et al. 2014;Kalbitzer et al. 2015; also see supplementary online material in Dal Pesco and Fischer 2018), our results reinforce the view that steepness and linearity are not a key element of male dominance in this species. The low levels of aggression and the absence of an agonistically enforced hierarchy show that rank is not a useful construct to predict interactions between males in this population (see also supplementary section S5c).
Although the lack of a dominance hierarchy is striking, these patterns are not uncommon in multilevel societies, where male status ("leader"/" non-leader" or "follower") is often used to define rank, with leaders being considered dominant over non-leaders, the latter of which tend to be younger males or older ex-leaders (e.g., Colmenares 1990;Bergman et al. 2009; but see Zhang et al. 2008 for a linear hierarchy between one-male-units in snub-nosed monkeys). We found a similar effect of age in Guinea baboons, where the individuals with the lowest rank scores were all large juvenile, adolescent, and old males (also see supplementary online material in Dal Pesco and Fischer 2018). Moreover, except in very rare cases, all bachelor males were large juvenile, adolescent. or late prime/old adult males, similar to observations in other multilevel societies (e.g., Colmenares 1990;Bergman et al. 2009). Hence, bachelor males can be found at the low end of the dominance hierarchy.
The finding that male Guinea baboons form strong, equitable and stable bonds corroborates the view that male coresidence and relatively low contest potential promote bond formation. Following predictions from kin selection theory, we found that strongly bonded males and primary and their bachelor males were on average more closely related. Future investigations should focus on the factors that promote malemale bond formation and the adaptive benefits of strong bonds between males, such as enhanced coalitionary support and/or increased male reproductive benefits. of the park Colonel Ousmane Kane for his cooperation and logistical support; all the staff and field assistants of the CRP Simenti, in particular Mustapha Faye, Armél Louis Nyafouna, Elhadji Dansokho, Lamine Diedhiou, Moustapha Dieng and Touradou Sonko for their support in the field. We thank Matthis Drolet, Lauriane Faraut, Mathieu Cantat, Fransziska Wegdell and Eréna Dupont for their assistance in the field. We also thank Roger Mundry for his precious statistical support and Matthis Drolet for invaluable comments on the manuscript. We would like to thank the editor and the two anonymous reviewers for all their useful comments, which greatly improved the quality of the manuscript.
Funding Open Access funding enabled and organized by Projekt DEAL. This research was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)-Project number 254142454/GRK 2070.
Data availability The data that support the findings of this study are openly available on OSF as "Data and scripts for: Kin bias and male pair-bond status shape male-male relationships in a multilevel primate society" via this link: https://osf.io/cfeqv/?view_only=908f656c01a741d1bf96ef55c0e4c5d9.

Compliance with ethical standards
Conflict of interest/competing interests The authors declare that they have no conflict of interest.
Ethical approval All applicable international, national, and/or institutional guidelines for the care and use of animals were followed. All procedures performed in studies involving animals were following the ethical standards of the institution or practice at which the studies were conducted. Approval and research permission were granted by the DPN and the MEPN de la République du Sénégal (research permit numbers: 0383/24/03/2009; 0373/10/3/2012). Research was conducted within the regulations set by Senegalese agencies as well as by the Animal Care Committee at the German Primate Center.

Consent to participate Not applicable.
Consent for publication All authors agreed to the publication of this manuscript.
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/.