Decoupling social status and status certainty effects on health in macaques: a network approach

Subjects and housing

The subjects of this study were 252 rhesus macaques (71 males, 181 females; age range: 3–29 years, mean = 7.7 years) from three outdoor social groups (126–185 individuals per group) that were part of the breeding colony at the California National Primate Research Center (CNPRC) in Davis, CA. Animals were socially housed in 0.2 ha outdoor enclosures containing multiple A-frame structures, suspended barrels, swings, and perches and were free to interact as they chose. Animals were fed a standard diet of monkey chow twice per day at approximately 0700 h and between 1430 and 1530 h. Fresh fruit or vegetables were provided one time per week and seed mixture provided daily. Water was available ad libitum. Animals housed in outdoor enclosures were managed with a minimum level of disturbance. At the end of the study all subjects were relinquished to the CNPRC breeding colony. This research was approved by the University of California Davis Institutional Animal Care and Use Committee.

Rank and dominance certainty data collection

Each group was observed for 5–7 weeks, one group in spring (Group A: March–April 2013) and two in fall (Group B: September–October 2013, Group C: September–October 2014) as part of a larger study on social networks and health. Unexpected management events resulted in lengthening (e.g. Group A) or shortening (e.g. Group B) the six-week study period by a few days in certain groups. We used an event sampling design to collect all instances of aggressive and submissive interactions. Two observers collected data for 6 h on 4 days per week from 0900–1200 and 1300–1600 h. Aggressive and submissive events were recorded as a series of dyadic interactions. A total of 13,567 events were recorded during 444 h of observation. Aggression was categorized according to severity and included threat (open mouth stare, brow flash, ear flap), mild aggression (threat and follow, lunge, push, slap, chase < 6 meters), moderate aggression (grapple, wrestle, chase > 6 meters), and intense aggression (pin or bite). Submission categories included freeze/turn away, move away, run away < 6 meters, run away > 6 meters, and crouch. Data on all dyadic aggressive interactions were used for constructing aggression networks to calculate dominance rank and certainty.

Dominance rank and certainty were calculated using the Perc package in R (Fujii et al., 2015) which uses a new network-based approach that combines information from direct dominance interactions with information from multiple indirect dominance pathways (via common third parties) to quantify dyadic dominance relationships (Fushing et al., 2011; Fujii et al., 2015). Essentially, individuals with dominance pathways that run counter to the primary direction of aggression (or submission) in the hierarchy have more ambiguous status than individuals who do not. The method begins by using a modified random walk algorithm to exhaustively identify all directed pathways in the network (e.g., A→B←C is not a directed pathway, but A→B→C is). To determine how much to weight the imputed (i.e. indirect) ‘wins’ from these pathways, the transitivity of the network is calculated as the proportion of transitive triangles (as opposed to cyclic or nontransitive triangles). In networks with high transitivity such as ours (> 95% of triangles are transitive, such that A > B > C and A > C), a ‘win’ that is imputed from an indirect pathway is more likely to reflect the true direction of dominance, and thus given higher weight, than in a network with lower transitivity. Regardless, wins from direct interactions are always weighted more than imputed ‘wins’ from network pathways. Finally, given the number of direct wins and imputed ‘wins’ from pathways, we calculate the probability that the row animal is dominant to the column animal using a Beta distribution to incorporate a source of statistical probability that reflects the level of uncertainty expected given the current data on each pair (for a more detailed explanation, see: Fushing, McAssey & McCowan, 2011). The matrix of these dyadic dominance probability values (range: 0–1) thus represent the cumulative information from all network pathways between each pair of animals. A dyadic dominance probability of one reflects the highest possible certainty that the row animal outranks the column animal, whereas 0.5 means the dominance relationship is perfectly ambiguous. The matrix of dyadic dominance probabilities was then used to generate the lowest cost linear rank order (see Fushing et al., 2011 for details).

Sparse and missing data (i.e. pairs of animals that are either infrequently or never observed to interact) are a common problem in animal behavior, and most ranking methods (such as the IS&I method (de Vries, 1998), Elo-rating (Neumann et al., 2011; Viswanath, 2016) and the Bradley-Terry model: Boyd & Silk (1983) and Albers & de Vries (2001)) are vulnerable to sparse data. The Percolation and Conductance method addresses the problem of sparse and missing data by gathering dominance information from network pathways. For example, in our study group A only 42.3% of all possible dyads had at least one agonistic interaction in the network. The remaining dyads were never observed to interact, making estimation of their pairwise dominance relationships from direct interactions prone to error. Furthermore, of the 42.3% that did interact, they averaged less than two interactions per dyad (x = 1.89). Adding dominance pathways dramatically increased the information per dyad to an average of 207.4 dominance pathways across all pairs (using up to 3-step pathways: A→B→C→D). Given that transitivity is very high, the dominance information from pathways has a high probability of agreeing with the true relationship. Furthermore, it is likely that macaques have the cognitive capacity to use this dominance information because many social vertebrates are capable of transitive inference such as deducing A > C from A > B and B > C (McGonigle & Chalmers, 1977; Davis, 1992; Bond, Kamil & Balda, 2003; Grosenick, Clement & Fernald, 2007).

From these dyadic dominance probabilities (which contain both rank direction and certainty information), an average dominance probability was calculated for each subject to provide a metric of the overall certainty of each animal’s rank. Prior to averaging, we transformed the dyadic dominance probabilities (initially bounded by 0–1) to range between 0.5 (ambiguous) and 1.0 (certain), thereby focusing on the information about certainty and ignoring rank direction. We also transformed ordinal dominance ranks for each group (derived from permuting the rows and columns of the dominance probability matrix) into the proportion of animals outranked within their respective groups (i.e. 0 is the lowest ranked animal and 1 is the highest ranked animal) to account for differences in group size. For the purposes of graphical representation, individual-level average dominance certainty was categorized into high, moderate, and low certainty categories for each group. We examined empirical distributions of average dominance certainty values (performed on each group separately) using the segmented package in R to identify logical break-points for these categories (Muggeo, 2008). Groups differed in their break-points due to differences in distribution, and these group-specific categorizations of low vs. moderate vs. high dominance certainty are thus reflected in our plots.

Social network measures and independence

Statistical analysis of social network metrics requires some attention to the lack of independence in network data. The calculation of a network metric for one node (e.g., an individual) typically involves some of the same edges as the same calculation for other nodes (Croft, James & Krause, 2008). As a result, a change in the value of one node could potentially affect the values of other nodes. We experimentally verified that this lack of independence was not an issue in this study by recalculating network measures after removing single individuals from the network. We selected three individuals with differing betweenness centrality (i.e. highest, median, and lowest) and experimentally removed each individual from the aggression network, one-at-a-time, recalculating the values of the remaining individuals each time, for a total of three recalculations (one for each individual). We chose to examine betweenness centrality because it is defined by shortest paths between all pairs in the network, similar to the counting of all paths between pairs in our Percolation and Conductance method. Bivariate and Spearman rank order correlations between the original dominance certainty and the new values indicate very little change (bivariate pairwise correlations for (a), (b) and (c): 0.999 ≥ r ≥ 0.932; Spearman pairwise correlations: 0.998 ≥ r ≥ 0.941). Furthermore, our Percolation and Conductance method incorporates statistical probability/uncertainty into the calculation of each dyadic dominance probability value, which serves to dampen any effect of dependence in the data.

Blood sample collection

Blood samples were collected on a single day (between 0800–1200 h) during the fifth week of the 5–7 week observation period for each group using the CNPRC’s standard method for biannual physical examinations. On sampling days, all animals in the group were immobilized (10 mg/kg of ketamine) and given standard physical examinations by veterinary staff (e.g., checked for injuries, weighed, assessed for pregnancy). Blood samples (5 mL) were also obtained from the femoral vein and serum was aliquoted and stored at −80 °C for later assay. Due to the large number of animals, blood was collected in batches of ∼15 samples. Batches were labeled with collection start and end times and animal identity to track the order of sample processing and control for effects of the capture and sedation procedure for animals processed later in the morning.

General health assessments

Health indicators were recorded twice weekly by one observer between 0900–1200 and 1300–1600 h. The observer located and visually inspected each study subject in the group and scored presence/absence of liquid stool (i.e. observing defecation of liquid stool or observing fresh liquid fecal matter on the tail or rump). From these data, we counted the total frequency of bouts of diarrhea per subject across the 5–7 week observation period. A bout of diarrhea was defined as either a single observation day with liquid stool or multiple consecutive observation days with liquid stool with the end of a bout marked by at least one observation day with no evidence of liquid stool. Occasionally animals were temporarily removed from the cage for veterinary care and were not available for health observations—these absences were recorded to control for total observation days.

Pro-inflammatory proteins

We measured three pro-inflammatory proteins (interleukin-6 (IL-6), tumor necrosis factor-alpha (TNF-α), and C-reactive protein (CRP)) from frozen serum. These proteins were chosen because they are markers of general systemic inflammation and are demonstrated risk factors for multiple diseases (e.g. type 2 diabetes and atherosclerosis, Pradhan et al., 2001; Libby, Ridker & Maseri, 2002).

Cytokine assay

Serum levels of IL-6 and TNF-α were measured using commercially available, species specific Milliplex multi-analyte profiling (MAP) reagents purchased from EMD/Millipore (Billerica, MA, USA), and utilizing Luminex Xmap technology (Luminex, Austin, TX, USA). Color coded polystyrene microbeads coated with specific antibodies for IL-6 and TNF-α were incubated with the serum samples, washed, and then were further reacted with biotinylated detector antibodies followed by Streptavidin-PE to label the immune complexes on the beads. After a final washing to remove all unbound material, the beads were interrogated in a BioPlex dual laser (BioRad, Hercules, CA, USA). The median fluorescent index for each sample was compared to a standard curve to calculate the concentration. Samples were tested in duplicate and had an intra-assay coefficient of variability of 15.3%. Samples falling below the threshold sensitivity of the assay (1.6 pg/mL) were assigned a value of one.

C-reactive protein assay

Concentrations of CRP were determined using a latex particle immunoturbidmetric method on the Beckman Coulter AU480 clinical chemistry analyzer.

Data analysis

Data analysis proceeded through a two-step process. First, hypothesized models (see below and Table 1) were tested for four dependent variables including three pro-inflammatory proteins (CRP, TNF-α, and IL-6) and one general health outcome (frequency of diarrhea bouts). The sets of hypotheses that guided model-fitting were explored for the following reasons. Age and sex were included as main effects in all models because both have been previously found to influence health (Klein, 2000; Sansoni et al., 2008). The impact of rank on health is also known to vary by sex in some populations (Creel, MarushaCreel & Monfort, 1996; Abbott et al., 2003), and it is reasonable to expect that the impact of rank may manifest differently at different ages. Therefore, sex by rank and age by rank interaction terms were also explored. We also examined a sex by dominance certainty interaction because the inherent structural differences in male vs. female rank acquisition (i.e. individual-level features such as age and body size drive male rank: (Dittus, 1975; Sprague, 1992; Sprague, 1998; Sebastian, 2015); family-level features such as agonistic support drive female rank: (Sade, 1972; Datta, 1986)) suggest that ambiguity may arise more readily amongst males than females and the presence of ambiguous relationships may have greater costs for females than males.

All data were analyzed using generalized linear mixed-effects regression models (McCullagh & Nelder, 1989; Raudenbush & Bryk, 2002) with group as a random effect. Models were run using a negative binomial distribution for IL-6, TNF-α and diarrhea bouts, and a gamma distribution for CRP. Appropriate distribution(s) for each outcome were chosen based on descriptive statistics, histograms, and Q-Q plots. We also evaluated the goodness of fit of these distributions using the Pearson chi-square statistic (SAS Institute Inc., 2009). We also included variables to control for known confounds: sample collection order for blood-based measures, and total health observation days for diarrhea. Due to the unavailability of blood samples for some animals, the cytokine and CRP analyses were run on 234 of the 252 subjects. Finally, ten of the study subjects showed CRP levels above 10 (a sign of possible infection), and these animals were included in analyses because they represent an important part of the health continuum. We note, however, that excluding these subjects from analyses did not change the magnitude or direction of the effects.

We used an Information Theoretic approach to determine which variables best explained each of our health outcomes. First, we ran all mixed-effects regression models to address our complete set of hypotheses outlined in Table 1. For each model, we present AICc scores, dAICc, model likelihoods (L = exp(−(1/2 * dAICc)), Akaike weights, and evidence ratios (ratio = weight of best model/weight of comparison model) as outlined in Burnham, Anderson & Huyvaert (2011). We then used model weights and dAICc to select a set of candidate models for each outcome variable. When the weight of the best model was < 0.90 (Burnham & Anderson, 2002; Symonds & Moussalli, 2011), we included in the candidate set those models with dAICc ≤ 5 and discuss the inferences based on all models in this candidate set. Analyses were performed in SAS 9.4 and R 3.3.1 (R Core Team, 2013). Plots were produced in Stata 14.1. Although plots of model results are typically generated from marginal effects tables, this was not appropriate for our data due to the presence of discontinuous regions. We have used the alternative method of generating model specific plots from predicted values (Hardin & Hilbe, 2007; Rabe-Hesketh & Skrondal, 2012).

Rank and dominance certainty

Rank and dominance certainty were associated in a nonlinear manner (Fig. 2). Specifically, animals of high- and low-rank exhibited higher dominance certainty than those in the middle of the rank distribution, a pattern that held true for all three study groups. In addition, there was variability in dominance certainty at each level of rank. Specifically, dominance certainty ranged from 0.70–0.97 and 0.72–0.93 among the highest and lowest ranked tertiles of the sample, respectively, which is very similar to the range for the sample as a whole (0.70–0.97).

Pro-inflammatory measures

Results for CRP indicated a single model with the highest weight (w = 0.93, H8; Table 2) in which the relationship between rank and CRP was dependent on dominance certainty and sex. High-ranking males with low dominance certainty had higher CRP, whereas little to no effect of rank was found in males with high dominance certainty (Table 3; Fig. 3A). For males with highly certain dominance relationships, increasing rank by 0.25 (moving up a quartile in rank) was associated with a reduction in CRP levels by 1.17 times. In contrast, for males with low dominance certainty, increasing in rank by 0.25 was associated with a 1.80 times increase in CRP. In females there was no change in CRP with increasing rank at any level of dominance certainty (Table 3; Fig. 3B). Consistent with predictions, older animals also had higher levels of CRP than adults, subadults, or juveniles (Table 3).

For IL-6 two models with similar weight accounted for nearly 80% of total weight (Table 4), and both models provided evidence that the effect of rank on IL-6 levels was dependent on dominance certainty (w = 0.45, H5; Fig. 4) and also potentially on sex (w = 0.33, H8; Fig. 5). Specifically, for animals with more certain dominance relationships, an increase in rank of 0.25 was associated with 1.41 times lower levels of IL-6 (H5), an effect that may be stronger in males than in females (H8; effect size: 2.16 vs. 1.14, respectively; Table 3; Fig. 5). Among animals with less certain relationships, an increase in rank of 0.25 was associated with a 1.50 times (H5) increase in IL-6 with effects possibly being stronger in males than in females (H8; effect size = 2.11 vs. 1.24, respectively). Taken together these results also indicate that high-rank animals with more ambiguous dominance relationships (i.e. low or moderate dominance certainty) exhibited higher levels of IL-6 than high-rank animals with more certain dominance relationships and that this effect may be specific to males.

Two simpler nested models were also part of the candidate set for IL-6 (Table 4). Model H3 (w = 0.08) included a main effect for dominance certainty plus age and sex terms, and indicated that an increase in dominance certainty of 0.10 was associated with 1.28 times lower levels of IL-6 (Table 3). Model H1 (w = 0.06) included only age and sex terms. For all four candidate models young animals (juveniles and subadults) had lower levels of IL-6 than adults (Table 3).

For the pro-inflammatory cytokine TNF-α, the model for H5 had the highest weight (w = 0.91; Table 5) indicating that the effect of rank on cytokine levels was dependent on dominance certainty for both males and females (Table 3; Fig. 6). For animals with more ambiguous dominance relationships (i.e. low dominance certainty) increasing in rank by 0.25 was associated with 2.07 times higher levels of TNF-α whereas increasing rank by 0.25 was associated with 1.78 times decrease in TNF-α for animals with highly certain dominance relationships (Fig. 6). Among low-ranked individuals, those with more certain dominance relationships had slightly higher levels of TNF-α compared to those with lower dominance certainty. Finally, juveniles and subadults had lower levels of TNF-α than adults (Table 3).

Plots of the raw data for CRP, IL-6, and TNF-α relative to rank and dominance certainty are presented in Figs. S1–S3.

Diarrhea bouts

Results for diarrhea bouts showed no single model with highest weight (Table 6). The set of candidate models included H4, H3, H7b, H2, H5 and H7a. Dominance certainty was included as a main effect only (i.e. not part of interaction terms) in models H4, H3 and H7a and was predictive of the incidence of diarrhea in each of these models (Table 7). According to these models, reducing the certainty of one’s dominance relationships by 0.10 (e.g., from 0.85 to 0.75) was associated with between 1.95 and 2.44 times greater incidence of diarrhea (Table 7), regardless of rank or sex. Models H5 and H7b were more complex, incorporating dominance certainty into interaction terms to determine whether its impact on health differed by rank or sex. The presence of these models in the candidate set suggest that the effect of dominance certainty might be dependent on rank or the effect of dominance certainty might be more pronounced in males (Table 7). However, addition of these interaction terms did not improve model fit, as AIC increased by 1.6 and 2.0 in H7b and H5, respectively (Table 6).

Dominance rank was included as a main effect only in models H2, H4, and H7b (Table 6). According to these models, reducing dominance rank by 0.25 (e.g. from outranking 75% to outranking 50%) was associated with between 1.25 and 1.42 times greater incidence of diarrhea (Table 7). Models H5 and H7a incorporated dominance rank into interaction terms with dominance certainty (described above) and with sex. Addition of the dominance rank × sex interaction did not improve model fit and AIC increased by 2.09 (Table 6). As expected, older animals showed a higher incidence of diarrhea bouts than adults and juveniles across all models (Table 7).

Rank and dominance certainty

Variation in social status, and thus its potential impact on health, includes not only whether one has high or low status but also the relative predictability or certainty of one’s status relationships. Dominance rank and dominance certainty showed a nonlinear (U-shaped) relationship where individuals of both high and low dominance rank had relatively greater dominance certainty than individuals in the middle of the hierarchy. However, variability in dominance certainty at all levels of dominance rank suggests that these two metrics may describe two complementary aspects of social status that can have independent and interacting effects on an animal’s health. The idea that dominance rank stability can have a critical impact on the status-health relationship has been suggested previously (Sapolsky, 1992; Sapolsky, 2005; Marmot & Sapolsky, 2014). The current paper, however, expands on previous research by using a more general measurement of each individual’s “fit” within the group hierarchy, as opposed to examining the status-health relationship in stable vs. unstable groups, or only in individuals that experienced dominance rank reversals.

Status and measures of inflammation

Our research suggests that whether individuals have high or low dominance rank in the hierarchy does not fully represent the complexity with which social status affects health. Instead, dominance certainty modified the impact of rank on biomarkers of inflammation. Under conditions of low dominance certainty, we found that high-ranking animals had higher levels of pro-inflammatory cytokines and CRP. This is consistent with previous findings in African wild dogs, baboons, chimpanzees, dwarf mongooses, and ringtail lemurs (Creel, MarushaCreel & Monfort, 1996; Cavigelli, 1999; Muller & Wrangham, 2004; Gesquiere et al., 2011), where high dominance rank individuals (or at least the alpha individual, as in Gesquiere et al. (2011)) had elevated GCs. Because dominance certainty is quantified as the consistency in the flow of dominance between pairs in the aggression network, having low certainty may be evidence of an individual-level tendency to protest or challenge others’ dominance rank or precursors of a rank change (within a hierarchy that is stable overall). Thus, high-ranking animals with ambiguous dominance relationships likely have reduced predictability that others will submit to them and greater potential risk of losing status (e.g., Crockford et al., 2008), and this type of uncertainty could represent a major psychosocial stressor. Our findings are consistent with previous research indicating that instability, whether due to dominance style (e.g. chimpanzees, Muller & Wrangham, 2004) or due to current social factors (e.g. baboon dominance rank reversals, Sapolsky, 1992), is particularly bad for high-ranking individuals, possibly due to potential for loss of status. For example, in an examination of rank reversal in male baboons, Sapolsky (1992) reported that despite similar levels of participation in aggressive interactions, only males about to lose their rank showed elevated GCs, not those about to rise in rank.

In contrast, for animals with a high degree of certainty in their dominance relationships, we found a small reversal of the effect with high-ranking animals having slightly lower levels of inflammatory markers than low-ranking animals. This effect is weak, but notably, more closely matches results from humans, baboons, macaques, meerkats, and spotted hyenas in which low status individuals exhibited risky health profiles (Sapolsky & Share, 1994; Shively & Clarkson, 1994; Goymann et al., 2001; Young et al., 2006; Ostner, Heistermann & Schülke, 2008; Gesquiere et al., 2011; Marmot & Sapolsky, 2014). In contrast to high status individuals, where uncertainty may have been associated with a potential loss in status, dominance uncertainty may not be a stressor for low-status individuals because they have little to lose; it may even be the case that low-ranking individuals with low dominance certainty may have an opportunity to increase their status—a possibility that needs to be explored further in future studies. Among low-ranking animals, being certain of one’s low dominance rank is likely to be stressful because such individuals can reliably expect to receive aggression, harassment, and/or intimidation from dominants, have little control over the occurrence of such interactions, and have fewer social outlets to cope with this harassment (Schino, 2001; Sapolsky, 2005). This is largely consistent with Shively & Clarkson’s findings that female macaques that were experimentally arranged to lose status (via group membership manipulations) increased atherosclerosis by 500%, whereas those that gained status increased atherosclerosis by a far smaller amount—44% (Shively & Clarkson, 1994). Thus, dominance rank changes may be stressful for both parties, but losing status is likely worse than gaining status.

Sex differences

Although sex differences in the status-health relationship are not often studied, sex differences have been noted previously (Kaplan & Manuck, 1999). Our results for the pro-inflammatory proteins, CRP and IL-6, suggest that there may be sex differences in the impact of rank on inflammation. Elevated levels of CRP in high-ranking individuals with lower dominance certainty was found only in males, not in females. The lack of an effect for females of this species is not surprising given the behavioral biology of rhesus macaques, specifically due to sex differences in how dominance rank is gained and maintained. Rhesus macaque males emigrate, and once they are established in a new group, they can increase their rank through both alliances and direct competition. Female macaques, on the other hand, remain in their natal group and inherit their dominance rank from their mothers with all the females of one family outranking all the females from another family (Lindburg, 1971; Sade, 1972; Melnick, Pearl & Richard, 1984). This process results in a male hierarchy that is more labile and changeable, but changes in dominance rank are rarer for females (Sade, 1972; Berard, 1990). For females, due to inherited rank and lifetime tenure in social groups, we would expect to see little dominance uncertainty, except during periods of larger scale group instability, such as when one matriline threatens to overthrow another matriline (Ehardt & Bernstein, 1986; Beisner et al., 2011).

Diarrhea

Dominance uncertainty was a better predictor of diarrhea than dominance rank with lower certainty being associated with greater risk of diarrhea. While animals with lower dominance certainty may be more vulnerable to diseases causing diarrhea, stress is also a known contributor to diarrhea. It may be the case that the increased stress of uncertainty in social relationships is contributing to the incidence of non-pathogenic diarrhea (e.g., Stasi et al., 2012; Buckley, O’Mahony & O’Malley, 2014). The small effect of dominance rank for diarrhea is consistent with other studies in baboons, macaques, meerkats, and spotted hyenas that demonstrate a wide range of negative biomarkers of health, as well as poor health outcomes in low-ranking individuals (Sapolsky & Share, 1994; Shively & Clarkson, 1994; Goymann et al., 2001; Young et al., 2006). Notably, unlike the biomarkers of inflammation discussed above, there was not strong evidence for an interaction between dominance rank and dominance certainty. This result suggests that the impact of social status on health is specific to the type of status measure being examined (e.g. certainty vs. rank) as well as the particular health or fitness outcome of interest (e.g. biomarkers of inflammation).