The relationship between individual phenotype and the division of labour in naked mole-rats: it’s complicated

Background The naked mole-rat (Heterocephalus glaber) is among the most social mammals on the planet, living in eusocial groups of up to 300 individuals that contain a single reproductive female and up to three reproductive males. A critical aspect of their complex social system is the division of labour that allows non-breeders to form an effective workforce. Age- or weight-based polyethisms are widely cited as explanations for how labour is divided, but evidence in support of these hypotheses has been equivocal. Methods To assess the extent to which individual working behaviour is determined by sex, age, weight and social rank, we studied the behaviours of 103 animals from eight captive colonies. We performed focal sampling and ran mixed-effects models to assess which factors explained variation in working behaviour during six ten-minute observation periods per individual. Results Contrary to widely-held beliefs, we found that working behaviour did not decrease linearly with weight, although polynomial regressions indicated younger and medium-sized individuals worked most frequently, while high-ranking individuals worked for the shortest periods of time. Working behaviour and its relationship with individual characteristics also varied between colonies. Conclusions While age- or size-based polyethisms may have some influence on working behaviour, we argue that other characteristics of the individual and colony are also important. In particular, the interactions of individual, social and environmental factors must be considered in order to understand the emergence and effectiveness of the division of labour that is so critical to many social organisms.


INTRODUCTION
The evolution of cooperation is a major driver of sociality and can provide the opportunity for labour to be divided so that individual contributions vary by task or by total effort (Maynard Smith & Szathmáry, 1999). The benefits of an effective division of labour have been studied in depth in social insects, in which a combination of behavioural flexibility and physical adaptations meet colony needs (e.g., Wilson, 1987;Gordon, 1989;Jeanson & Weidenmüller, 2014). In contrast, where tasks are divided in social mammal groups, the underlying causes of this division are poorly understood. After almost two decades without further research, a recent analysis has raised further questions regarding earlier findings. Mooney et al. (2015) showed that digging behaviour was not associated with body mass or age, although aggression towards foreign mole-rats was positively correlated with body fat. In light of these results, we performed a largescale study of the determinants of working behaviour in naked mole-rats. By observing individuals across eight captive colonies, we examined the effects of weight and age on working behaviour. Additionally, we tested for the effect of rank, a characteristic that could influence how much an individual contributes to collective activity but which has not been thoroughly tested (but see Clarke & Faulkes, 2001).

MATERIALS & METHODS
We selected individuals from eight colonies of captive-bred naked mole-rats at Queen Mary University of London and the study was carried out in accordance with institutional guidelines. They were all captive born descendants of animals that were originally captured in Kenya during the 1980s. Our youngest subjects were four months old during observations, as naked mole-rats do not exhibit the full behavioural repertoire until around three months (Jarvis, 1981;Jarvis, 1991). Details of the individuals and colonies used in the study are given in Table 1 (Results, below). Each colony was housed separately in a single room, in a network of acrylic tunnels and boxes. The temperature of the room was maintained between 26-32 • C. Ambient light and noise were constantly recorded and changed unpredictably as people moved around the room and surrounding area. The animals were fed ad libitum every day. The diet generally consisted of carrot, sweet potato and butternut squash, although occasionally other fruits and vegetables were provided.
We observed all individuals in colonies that contained fewer than 20 in total, and randomly selected 20 individuals from each of the remaining colonies. The observation sequences of colonies and individuals within colonies were also randomly generated. If a colony's breeding female was not selected randomly, she was added to the colony's observation list manually (colonies 11A, 11C and FK100) as a means of comparison in our focal studies. Each individual was observed three times in the morning and three times in the afternoon to avoid confounding variation from daily behaviour patterns. Animals were weighed in each of the three weeks they were observed, and means were calculated and used in the analyses. Animals were identifiable through RFID microchips and were marked twice a week with black marker pens so we could distinguish between them during observations.
Females were classified as breeders if they had a perforated vagina during the observation period, they were seen mating, or became pregnant, and we classified all other females as non-breeders. External genitalia are generally monomorphic in non-breeding naked molerats (Jarvis, 1991;Peroulakis, Goldman & Forger, 2002). However, the dark-red vaginal membrane can become prominent in some individuals, which is thought to reflect a partial release from reproductive suppression and may be a sign of incipient reproductive activation (Jarvis, 1991). We recorded which females had prominent vaginal membranes and excluded them from the analyses. Males seen mating were recorded as breeders and the rest were classified as non-breeders (Jarvis, 1991). While we recorded ano-genital nuzzling when possible, we did not classify animals as breeders if they were observed in these interactions as some ano-genital nuzzling involves non-breeders (Jarvis, 1991). Historical records of breeding status based on this classification were also used.
On entering the animal room, we allowed the colony to settle down from any increased behavioural activity for a period of up to 10 min and we minimised movement and noise while observing to avoid disturbing the animals. Behaviours were recorded between 08.00 and 18.00. Individuals were observed sequentially in their burrow system for ten minutes each, with each individual observed six times for a total of 60 min per individual.
The broad behavioural categories we recorded were as described in the ethogram of naked mole-rat behaviour: transport of food and nest material, digging and offspring tending . As it was not possible to distinguish between digging and transporting behaviours, these were classified together as working behaviour. Offspring were only present in one colony during observations and the few instances of tending offspring were not included in the analyses.
We tested for an effect of dominance rank on working behaviour and aggression by establishing the dominance hierarchy within each colony using passing behaviour, which is a reliable indicator of dominance hierarchies in naked mole-rats (Clarke & Faulkes, 1997). We recorded which individuals passed over the top of other individuals during face-to-face encounters in tunnels. Interactions not thought to indicate rank include tail-to-face encounters, passing in chambers or corners of tunnels, when one individual digs throughout the encounter, and when individuals do not pass directly over the top of one another.
We used the Elo rating system to calculate individual dominance ranks within colonies (R package EloRating, Neumann & Kulik, 2020). The Elo rating method has several advantages over matrix-based methods such as its ability to calculate ranks within small groups and account for the loss of individuals during observations. After calculating each individual rank, we scaled each rank to account for variation in group size by dividing by the number of individuals observed in the respective colony. Lower rank values (those towards zero) represent more dominant animals. We calculated the steepness of each colony's dominance hierarchy using the steepness function from the EloRating package, which is based on David's Scores (De Vries, Stevens & Vervaecke, 2006).
Ten-minute observation periods that contained no working behaviour were assigned a value of zero. To investigate which factors predicted whether an individual showed any working behaviour, we conducted logistic mixed-effect models, using the R package lme4 (Bates et al., 2020, p. 4). For individuals that showed working behaviour in a given session, we determined the factors that predicted variation in the duration of working behaviour using linear mixed-effect models, also in lme4. For a discussion of ''two-part'' models, see (Duan et al., 1983;Min & Agresti, 2002). We checked the residuals of the logistic models for uniformity, dispersion, zero-inflation and the presence of outliers using the DHARMa R package (Hartig, 2020). For the linear models, we confirmed the residuals were normally distributed and had similar variances.
In both our logistic and linear models, we first constructed null models in which we included individual and colony fitted as random effects, with individual nested withincolony. The response variable was working behaviour per observation session. For each null model, we then constructed a set of separate models, each containing one of the individual characteristics as a fixed effect: sex, age, weight and rank. Due to the correlations between age, weight and rank (older individuals tend to be bigger and of higher rank, Schieffelin & Sherman, 1995;Clarke & Faulkes, 1997;O'Riain & Jarvis, 1998); we did not create models that contained more than one fixed effect. As relationships between work and age, weight and rank could be non-linear, we also created polynomial linear regressions that included the individual characteristic-squared and -cubed. We compared the performance of the models as described below and report the best performing for each characteristic.
Akaike Information Criterion (AIC) values estimate how well a model approximates the unknown reality relative to other models; smaller AIC values indicate better models (Burnham & Anderson, 2002). Second-order AICs (AICcs) were generated using the aictab function from the AICcmodavg package (Mazerolle & Linden, 2020) for each model and compared to see whether the inclusion of an individual characteristic reduced the AICc. The relative likelihood of a model given the data is calculated by exp(-( 1 2 ) AICc), where AICc is the difference in AICc between two models (Burnham, Anderson & Huyvaert, 2011). The ratios of model likelihoods can be used to calculate an evidence ratio, which indicates the extent to which the data support one model over another (Burnham, Anderson & Huyvaert, 2011). We report the AICc values, AICcs and likelihoods relative to the null model for each alternative model. Along with AIC values as indicators of relative model performance, we report Nakagawa's R 2 for each model to estimate how much variation in our working behaviour data is explained by the independent variables (Nakagawa & Schielzeth, 2013;Nakagawa, Johnson & Schielzeth, 2017), calculated using the performance R package (Lüdecke et al., 2019).
All statistical analyses were carried out in R Studio Version 3.6.0 (R Core Team, 2014) and figures were made using the ggplot2 package (Wickham, 2016). To assess which characteristics predicted whether working behaviour was observed or not, we compared each of our logistic regression models with a single fixed effect to the null model and assessed model fit. AICcs showed that sex, weight and rank did not improve model fit ( AICcs, versus null model: Sex = +1.48, Weight = +0.79, Rank = −0.06; model outputs in Table 3, SI5). Effect sizes and standard errors support the inference that these characteristics do not explain variation in whether an individual was observed working. In contrast, adding age was associated with a reduction in AICc value of 7.24 and the likelihood of this model was approximately 37 times higher than that of the null model. The model coefficient suggests individuals were less likely to be observed working as age increased and this result was significant at alpha = 0.05. While including weight did not improve the performance of the model (SI5), including weight-cubed reduced the AICc by 3.89 in a model that had a likelihood 7 times higher than the null model. Nakagawa's R 2 estimates suggest all models had limited explanatory power (conditional R 2 between 0.202 and 0.230), including those that performed better according to AICc, and more variation was explained by the colony and individual random effects than by any of the fixed effects (Table 3).
To assess which characteristics predicted the duration of work observed, we compared each of our linear models with a single fixed effect to the null model and assessed model fit. Linear models excluded observation periods during which no working behaviour was recorded. AICc values suggest adding sex, age, weight or rank as predictors did not improve the null model, which included only colony as a random effect ( AICcs, versus null model: Sex = +2.04, Age = +2.00, Weight = +1.28, Rank = +0.56; model outputs in Table 4, SI5). Effect sizes and standard errors support the inference that these characteristics do not predict variation in the amount of time individuals were observed Table 3 Results of logistic mixed-effects models predicting whether working behaviour was observed or not. The format of the models was: Presence or absence of working behaviour ∼ Colony|Individual + Fixed Effect where Presence of working behaviour is whether working behaviour was observed during the ten-minute observation period, Colony and Individual are random effects with Individual nested within Colony, and Fixed Effect was omitted in the null model and one of sex, age (months), weight (grams) and rank (0-1, scaled) was included in the corresponding models. The scripts used to run the models are available in File S3. Relative likelihoods are calculated as exp(-( 1 2 ) AICc). Nakagawa's R 2 estimates the variance explained by the fixed effects (marginal variance) and both fixed and random effects (conditional variance) (Nakagawa & Schielzeth, 2013  The plots are fitted with locally weighted (loess) regression lines that display localised trends in the data. These are different from the models described below which include colony and individual as random effects. The size of the points reflects the number of data points at each location (geom_count, gg-plot2). Rank is scaled to account for group size and the most dominant individuals have ranks closer to zero. Age and weight-cubed were associated with the probability that an individual was observed working; younger and mid-sized individuals worked most frequently. The plots are fitted with locally weighted (loess) regression lines that display localised trends in the data. These are different from the models described below which include colony and individual as random effects. Rank is scaled to account for group size and the most dominant individuals have ranks closer to zero. Rank-squared was associated with the duration an individual was observed working; more dominant animals worked for shorter periods of time.
Full-size DOI: 10.7717/peerj.9891/ fig-2 working. While including rank did not improve the performance of the model (SI5), including rank-squared reduced the AICc by 7.21 in a model that had a likelihood 37 times higher than the null model. Figure 2C suggests higher ranking individuals may work for shorter periods, although this effect appears to plateau and may even reverse among mid-and low-ranking animals. Nakagawa's R 2 estimates suggest all models had limited explanatory power (conditional R 2 between 0.169 and 0.182), including those that performed better according to AICc, and more variation was explained by the colony and individual random effects than by any of the fixed effects (Table 4). In summary, age did predict whether an individual worked or not, but was not associated with the duration over which an individual worked. None of the other variables (sex, weight and rank) predicted either whether an individual was recorded working or not, or the duration of working behaviour. Weight-cubed also predicted whether an individual was observed working, while rank-squared predicted the duration of working behaviour. None of the models had high explanatory power and the majority of the variation in working behaviour was unexplained.
Plots of the duration of working behaviour and the relationships between working behaviour and individual characteristics within-colonies demonstrate the between-colony variation in our data . Each colony's dominance hierarchy was assessed for how steep it was and assigned a value between 0 and 1, with 1 being the steepest. All steepness estimates were below 0.2.

DISCUSSION
Previous studies indicate that some naked mole-rats spend more time working than others (Jarvis, 1981;  and that differences between individuals may be stable for long periods of time (Mooney et al., 2015). Our analyses of 618 observation periods, encompassing 103 non-breeders from eight colonies, indicated that the probability of being observed working was explained by age and weight-cubed but not by sex, weight or rank. Initially, work frequency seems to increase with age and weight, until intermediate weights and around the age of two, after which the frequency of work plateaus and may even decrease (Fig. 1). The assignment of different roles to different age or weight classes distributes work across a colony without the need for costly cognitive evaluations of colony needs and worker availability based on patchy information (Robinson, 1992). Indeed, the presence of overlapping generations of offspring is considered a key criterion for eusociality (Crespi & Yanega, 1995) and different age and weight cohorts will almost always be present within eusocial groups. Despite the effect of age on probability of working, when we excluded periods in which no working behaviour was recorded, we found that only rank-squared predicted the amount of time spent working during the observation period. Figure 2C suggests low-and mid-ranking individuals worked for longer periods. Work duration was shorter among high-ranking, more dominant individuals.
In the absence of clear linear predictors, we suggest that age-or size-based polyethisms may depend on other conditions, such that colony members alter their working behaviour in response to the complex interaction of individual-and group-level pressures. Indeed, Table 4 Results of generalised linear mixed-effects models predicting the duration of observed working behaviour. The format of the models was: Duration of working behaviour ∼ Colony|Individual + Fixed Effect Where Duration of working behaviour is the duration of working behaviour observed during the ten-minute observation period, Colony and Individual are random effects with Individual nested within Colony, and Fixed Effect was omitted in the null model and one of sex, age (months), weight (grams) and rank (0-1, scaled) was included in the corresponding models. The scripts used to run the models are available in File S3. Relative likelihoods are calculated with exp(-( 1 2 ) AICc). Nakagawa's R2 estimates the variance explained by the fixed effects (marginal variance) and both fixed and random effects (conditional variance) (Nakagawa & Schielzeth, 2013  Mooney et al. (2015) found that individuals changed their behaviours to compensate for the removal of their colony mates, indicating a level of flexibility in response to the needs of the group. Moreover, a relationship between age and behaviour was only evident when frequent-workers were removed, for which younger individuals compensated by increasing their work rate (Mooney et al., 2015). On the other hand, in the eusocial Damaraland molerat (Fukomys damarensis), age does seem to play a key role in cooperative behaviour (Zöttl et al., 2016;Thorley et al., 2018). Both studies of Damaraland mole-rats found helping behaviours increased until the age of one, after which point there was either a plateau (Zöttl et al., 2016) or reduction (Thorley et al., 2018) in helping behaviour. The similarity of helping behaviour in Damaraland and naked mole-rats is interesting, and implies convergent evolution of these patterns given that sociality is thought to have evolved separately in the two species (Faulkes & Bennett, 2013). Taken together, our findings imply that older individuals perform fewer bouts of working behaviour, although these bouts do not differ in duration compared to those of younger individuals. In mammals, older individuals tend to be less active (Ingram, 2000;Marck et al., 2017) and our results may reflect this general mammalian trend. There is also evidence activity decreases after the of age three in the Ansell's mole-rat (Schielke, Begall & Burda, 2012). The relatively small number of working bouts seen in older mole-rats might also relate to them having other roles, which were not recorded here. For example, while we did not analyse aggressive or defensive behaviours, it is possible that older individuals, who also tend to be larger (O'Riain & Jarvis, 1998), spend more time defending the colony (Mooney et al., 2015). Older individuals may also invest less in working behaviour in order to maintain energy reserves for future challenges for dominance within the colony, particularly given that these animals are more likely to succeed in the event that a breeder dies or is removed (Reeve, 1992;Clarke & Faulkes, 1997;Clarke & Faulkes, 1998).
We found no evidence that work decreases linearly with body size, but work frequency increases with size until intermediate body weights are reached. After this, work frequency plateaus and may decrease. However, we are cautious in our interpretation of this result due to the relatively sparse data from animals of small and large body weight, and the tendency of polynomial models to overfit data (Lever, Krzywinski & Altman, 2016).
Across the African mole-rats, the link between age, size and working behaviour is still far from clear. For example, in naked mole-rats, Mooney et al. (2015) also reported no relationship between body size and working behaviour, although larger individuals were more aggressive towards conspecifics (O'Riain & Jarvis, 1997). On the other hand, body size does seem to be a good predictor of space-use in the Ansell's mole-rat (Šklíba et al., 2016). In contrast, a reported association between size and working behaviour in the social Micklem's mole-rat (Fukomys micklemi) was attributed to an underlying age-based polyethism (Van Daele et al., 2019), while a similar association in the Damaraland mole-rat may relate to growth rates (Zöttl et al., 2016). An early report on Damaraland mole-rats found that heavier individuals worked more (Gaylard, Harrison & Bennett, 1998). In this respect, it is interesting to note that naked mole-rats exhibit unusual growth patterns; growth rates vary between litters and individuals can gain size rapidly in response to social changes, even after extended periods of constant size (O'Riain & Jarvis, 1998). This decoupling of size and age may further complicate the relationship with behaviour.
We found no linear effect of social rank on working behaviour, but did find an effect of rank-squared. Naked mole-rat groups have strict dominance hierarchies (Clarke & Faulkes, 1997) and periods with high rates of aggression and fighting during which individuals are often killed (Clarke & Faulkes, 1997;Clarke & Faulkes, 2001;Medger et al., 2019). Highranking naked mole-rats are more likely to become breeders if the previous breeders are removed (Clarke & Faulkes, 1997;Clarke & Faulkes, 1998), and may reduce energy spent on working behaviour in order to prepare for future dominance challenges. Alternatively, high-ranking individuals tend to be older and heavier (Clarke & Faulkes, 1997;Clarke & Faulkes, 1998) and may contribute less to general working behaviour but more to colony defence, as has been reported elsewhere (O'Riain & Jarvis, 1997;Mooney et al., 2015). More data are need before we can be certain about the distribution of work among heavier, older and higher-ranking individuals. In our study, we recorded an average of four interactions per individual, whereas the Elo-rating method for assigning rank performs better when there are around ten interactions per individual (Sánchez-Tójar, Schroeder & Farine, 2018). Thus, given the large confidence intervals around estimates of rank, we would like further evidence to confirm the impact of rank, which might be more evident over longer observational periods.
Our observations reveal considerable variation in the total amount of working behaviour and the relationships between work and individual characteristics within each colony . With data from eight colonies, we did not have the power to test colony-level hypotheses, however, individual behaviours may be affected by a number of colony-level factors, such as food availability (Reeve, 1992), worker availability (Mooney et al., 2015) or the age-structure of the colony (Gaylard, Harrison & Bennett, 1998;Šklíba et al., 2016). Many published results from naked mole-rat behaviours have come from studies of one or a few colonies and were thus unable to account for colony-level variation. In this study, more variation in working behaviour was explained by individual and colony variables than by any individual characteristics. Our findings indicate that future research should increase the number of colonies used and statistically control for colony-level effects wherever possible.
This study looked at the factors that could predict whether an individual works and the variation in the amount of time individuals spend working. Labour could also be divided by task. Task specialisation has been observed in some social insect species (Charbonneau & Dornhaus, 2015) and, although it could be an important benefit to sociality in mammals, it has rarely been recorded (Stander, 1992;Gazda et al., 2005;Hurtado, Fénéron & Gouat, 2013;Gazda, 2016). Mooney et al. (2015) found that the amount of time spent on different tasks was stable in naked mole-rats, although Thorley et al. (2018) found no evidence for task specialisation in the eusocial Damaraland mole-rat. We could not test worker specialisation due to the lack of pup tending observed. Although some authors have shown that repeated exposure to a task increases efficiency in insects (Langridge, Sendova-Franks & Franks, 2008), this is not always the case (Dornhaus, 2008;Santoro, Hartley & Lester, 2019). Future research could use experimental tasks to explore the presence of specialisation in naked mole-rats and determine whether it increases individual or colony efficiency. Assessing worker, pup tending and defence behaviours like Mooney et al. (2015) would be the best way to establish whether task specialisation exists, as has been suggested in studies of other mole-rats (Van Daele et al., 2019). Gordon (2016) suggested that the term 'division of labour', implying internally-driven choices, is too rigid, and instead terms relating to individual behaviours should emphasise the influences of external conditions and social interactions. Further, we agree that workers should not be classified according to discrete 'castes' in mole-rats unless clear supporting evidence is reported (Šklíba et al., 2016). Polyethisms based on individual characteristics could be the start of an effective division of labour but the relationship appears to be complicated and variable. The outstanding questions are the extent to which polyethisms can be fine-tuned and the mechanisms that facilitate these changes. To date, the causes of within-individual variation are still largely unstudied in social mammals and insects (Jeanson, 2019). Given our results, combined with those of Mooney et al. (2015) and the mixed findings of earlier studies (Faulkes et al., 1991;Jarvis, 1991;, future research should focus on the interaction between internal, social and environmental