Two simple movement mechanisms for spatial division of labour in social insects

Automated tracking of four social insect species

Fifty queenright colonies were used in the tracking experiments (Table 1). Honeybee colonies (subspecies A. mellifera carnica) were housed in the campus apiary of the University of Lausanne. Colonies of L. niger were raised from single mated queens collected on campus. T. nylanderi colonies were collected from the University of Lausanne campus, and L. acervorum colonies collected from Anzeindaz, Switzerland. These four species were chosen because of their abundance and easy availability in Switzerland, and because they – or closely-related species – have previously been used as model systems for the study of spatial organisation in social insects^{20,21,23,27,44}. The colony sizes used in our experiments (Table 1) fell within the natural range of sizes experienced by these species in nature, either as recently founded colonies (L. niger colonies are founded by a single queen and progressively grow from a few workers to mature sizes of up to 40,000 workers over the course of several years; new honeybee colonies are founded by swarms counting 2400–41,000 bees⁶¹) or as mature colonies (all colonies of T. nylanderi and L. acervorum used in our experiments were mature colonies collected whole from the field).

In all species, a paper tag bearing a unique two-dimensional barcode was glued to the thorax of individuals to allow automated tracking of their movements (Fig. S1). In the ants, tagging of all individuals was performed in a single session two days before the beginning of the experiment, whilst in the bees, newly-emerged workers (one-day-old or less) were tagged every 3 days over the 21 days prior to the beginning of the experiment (Supplementary Note 1).

Tagged colonies were kept in glass observation nests with a single entrance (internal nest dimensions, A. mellifera: 69 × 45 × 4 cm, L. niger: 70 × 40 × 8 mm; L. acervorum: 63 × 42 × 2 mm, T. nylanderi: 63 × 42 × 1.5 mm). The honeybee observation nests also included a 64 × 44 cm wooden frame enclosing a double-sided wax comb containing honey, pollen, and developing brood²⁰. Bees were free to move between both sides of the comb. In all species, individuals were allowed to freely exit and enter the nest. Ants were provided with ad libitum food (Drosophila, sugar solution) and water in a foraging arena, while bees foraged on natural resources outside. Both the ant and honeybee observation nests were exposed to diurnal cycles of temperature and light (Supplementary Note 1).

High resolution digital video cameras operating at two frames per second were used to identify the location and orientation of each tag across successive images²². All colonies were continuously tracked for three days, which corresponded to the inter-cohort time in the honeybee colonies. The trajectories of each worker, and the physical contacts between workers (Fig. S18 and Supplementary Note 14) were extracted using an existing software pipeline⁶².

Building bipartite site-visit networks

To quantify the spatial preferences of individual ants and bees, the interior of the nest was discretised into a regular hexagonal lattice (Fig. 1a, b). Because the worker body lengths of our four study species span an order of magnitude (from ~ 1.5 mm for T. nylanderi to ~ 15 mm for A. mellifera), the width of the hexagonal bins were defined as 1/4 of the mean worker body-length.

To characterise the spatial preferences of different individuals to different parts of the nest, we counted the number of times ({n}_{i}^{s}) that each individual i visited each hexagonal site s. A visit by individual i to site s began when i crossed the border into s, and was terminated when i crossed the border out of s, regardless of the amount of time spent inside. To prevent stationary individuals located on the border between two adjacent sites from rapidly accumulating many single-frame visits to the two sites, successive visits to a same site were only counted when at least 20s elapsed between the end of the previous visit and the start of the next.

The site-visit data were used to construct a bipartite network, in which individuals (layer 1) were connected by undirected edges to the sites (layer 2) they visited (Fig. 1c, d). Because individuals typically made multiple visits to the same sites, each edge i–s was weighted according to the total number of times individual i visited site s, that is, ({n}_{i}^{s}).

Partitioning site-visit networks into modules

The extent to which the site-visit networks were partitioned into discrete ‘modules’ (i.e., set of workers with similar space-use patterns and the set of sites that they exhibit strong ties to) was assessed using the DIRTLPAwb+ algorithm for partitioning weighted bipartite networks³⁹. This algorithm searches for the partition that maximises the number and strength of the links within modules, whilst minimizing connections between modules. The number of modules was not specified a priori by the user, but was identified by the algorithm. All site-visit networks had positive modularity (Fig. S3), indicating that they could be partitioned into a set of well-separated modules (Figs. 1e–h, S2, and S4–S5). The modules in each partition were then assigned functional labels according to the following rules. First, the module whose sites were on average closest to the nest entrance was labelled ‘forager’ module. Second, the module or modules with the greatest spatial overlap with the brood pile in the ant colonies or the broodnest(s) in the honeybee colonies were labelled ‘nurse’ module(s). After defining the forager and nurse modules, the remaining modules (if any) were labelled as follows. If there was only one module remaining after identifying the nurse and forager modules, as was typically the case in honeybee colonies, it was labelled ‘peripheral’. If there were two modules remaining, as was typically the case in ant colonies, then the module whose sites were on average closer to the nest borders (i.e., to the periphery of the nest) was labelled ‘peripheral’, and the remaining module labelled ‘intermediate’. In some cases, the DIRTLPAwb+ algorithm identified five or more modules (9.0% of all iterations across all species and colonies). In these cases, the supernumerary modules never contained more than 1 or 2 individuals, and as they could not be unequivocally assigned using our labelling scheme, they were left unclassified for these iterations.

Validating network modules

As a network constructed by a purely random process could exhibit apparent modular structure by chance, we tested whether the discovered modules represent statistically significant entities. To do so, we produced 1000 null model random networks for each observed network using an established permutation method for bipartite networks⁶³ (Supplementary Note 2). Comparisons between the maximum modularity of the observed networks with that of the corresponding random networks showed that, in all four species, the observed modularity was significantly greater than expected by chance (Fig. S3).

Constructing worker task profiles

A unique labour profile for each ant and each honeybee was constructed by estimating the activity of each worker in the following four tasks:

1. Entrance guarding: workers were classed as guarding when they were (i) within two body lengths of the entrance, (ii) roughly facing the entrance, i.e., with a body alignment diverging from the direct heading to the entrance by no more than π/2 radians, and (iii) ‘on station’ at the entrance, as defined by trajectory coordinates with an associated first passage time (ref. 64; time taken for the individual to pass beyond a circle centred on its current location with a radius of two body-lengths) in excess of 500s.

2. Patrolling: workers were classed as patrolling⁶⁵ when they were (i) active, and (ii) ‘roaming’, as defined by first passage times of <5min for a circle with a radius of four body-lengths.

3. Queen attendance: workers were classed as attending the queen if they were in physical contact and facing towards the queen, as defined by the trapezoid method for identifying contacts described in Supplementary Note 14.

4. Foraging: workers were defined as foraging when they left the nest and entered the foraging arena.

As the total time an individual allocates to a given task, and the number of times it performs that task can vary (nearly) independently, we quantified both the total time spent on each task, and also the number of bouts of each task. As individuals were occasionally lost from view all measures were normalized by the total number of trajectory fixes for a given individual.

Assigning module scores for individuals and sites

As discrete categories are not always suitable to describe continuous biological processes, we developed a methodological extension that allows for overlapping modules. To do so, we exploited the stochastic nature of the DIRTLPAwb+ algorithm and applied it to each network 1000 times, thus producing an ensemble of slightly different partitions. After eliminating duplicate partitions, we defined the score denoting the membership of each individual i (or site s) to each module M as the proportion of partitions in which i (or s) was assigned to M. Thus, each individual and each site was assigned a set of module scores summing to 1 (Supplementary Note 2).

Measuring diversity

We used an information-theoretic approach to obtain three measures of the heterogeneity of the module scores for individuals and spatial locations. These were termed the module score diversity of (i) the visitors to a given site, (ii) the sites an individual visited, and (iii) the nestmates an individual contacted. To do so, we first defined the typical module score profile of (i) the typical visitor to a given site, (ii) the typical site an individual visited, and (iii) the typical nestmate an individual contacted by calculating the following weighted averages:

where [N^*, I^*, P^*, F^*] denotes the module score profile of (i) the typical visitor of a given site s, (ii) the typical site visited by an individual i, or (iii) the typical nestmate contacted by individual i; n_i,s denotes the number of visits by individual i to site s; c_i,j denotes the number of contacts between individuals i and j; [N, I, P, F] denotes the scores of (i) individual i, (ii) site s, or (iii) individual j for the Nurse, Intermediate, Peripheral and Forager modules. Because these are weighted averages, there is a proportionally greater contribution to the typical score profile by (i) individuals that visit site s more frequently, (ii) sites that individual i visits more frequently and (iii) individuals that individual i contacts more frequently. These typical module score profiles were then used to calculate the module score diversity D of each site s and individual i;

where H_visitor,s is the entropy of the typical scores of the typical visitor to site s, ({[{N}^{*},{I}^{*},{P}^{*},{F}^{*}]}_{visitor,s}); H_visited,i the entropy of the typical scores of the typical site visited by individual i, ({[{N}^{*},{I}^{*},{P}^{*},{F}^{*}]}_{visited,i}); and H_contacted,i the entropy of the typical scores of the typical nestmate contacted by individual i, ({[{N}^{*},{I}^{*},{P}^{*},{F}^{*}]}_{contacted,i}). The diversity index ranges from 0 (e.g., sites that are only visited by specialists for one module – individuals that scored 1 for one module and 0 for all other modules) to 1 (e.g., sites that are visited by module generalists and/or by an equal mix of specialists from all modules).

Module score similarity between workers and sites

The similarity between the module scores of a worker and a site that it visits was quantified using cosine similarity, a standard measure of the distance between two vectors, defined as the cosine of the angle between them. Each node was represented by the vector of its module scores, and the cosine similarity of individual i and site s, was defined as follows:

where M_i and M_s are the module scores of individual i and site s for module M, respectively. Because all module scores were positive, in our study cosine similarity ranged from 0 (orthogonal vectors) to 1 (identical vectors).

Measuring long-range attraction

To assess whether long-range attraction of workers towards their primary module could be a mechanism for producing spatial segregation between workers, we defined a taxis index T_M,i,s, which measures the attraction of an individual i to a focal module M, at each site s that it visited. The taxis index was defined by calculating the projection of the mean resultant vector of all trajectory segments of i starting at site s, (overrightarrow{{v}_{i,s}}), onto a vector pointing directly from site s to the nearest point on the boundary of the focal module M (Fig. S11). Thus, positive values of T_M,i,s occur when an individual exhibits a tendency to head toward M, and negative values indicate a tendency to move away from M.

To test whether individuals display greater attraction towards their primary module than individuals belonging to other modules, we then defined, for each module M and each site s lying outside M, a mean taxis index (overline{{{{{{{{{rm{T}}}}}}}}}_{M,s}}) for either ‘resident’ workers (i.e., individuals whose primary module is M) or ‘non-resident’ workers (i.e., individuals whose primary module is not M). For each species, we then tested whether worker ‘residency’ (resident vs. nonresident) had a significant effect on the mean taxis index (overline{{{{{{{{{rm{T}}}}}}}}}_{M,s}}) at different distances from the module boundary using generalised mixed-effect models with (overline{{{{{{{{{rm{T}}}}}}}}}_{M,s}}) as a dependent variable, worker group and colony size as fixed effects, and colony identity, site identity and focal module identity as random effects. As this analysis involved multiple testing (one for each distance), p-values were adjusted using the Benjamini-Hochberg (BH) correction for multiple testing.

Investigating location-dependent movement

We here describe three fundamental measures of individual locomotion. The first was the site-specific activity probability, which describes the likelihood that a particular individual is in motion at a particular site. The activity probability was obtained by first decomposing each trajectory into an alternating sequence of active and inactive bouts using a combination of change-point analysis⁶ and cluster analysis (Fig. S19, Supplementary Note 15). For each individual i – site s pair, we then calculated the activity probability P_i,s, that is, the proportion of time that i was in the active state when visiting site s (Fig. S20 and Supplementary Note 16).

To further characterise locomotion while individuals were in the active state, we defined two additional measures for every site s, visited by each individual i, namely the mean speed while active v_i,s, and the mean unsigned turn angle while active θ_i,s (Fig. S20). Calculating these measures using only trajectory coordinates associated with active individuals ensured that these measures capture properties of movement, rather than the probability of moving.

Mapping module gradient fields

To estimate the overall module gradient field for a given module M, we performed a spatial multiple regression at each site s. To do so, we first defined a local neighbourhood around each site s, which included all sites whose centres were within one worker body length of the centre of the focal site. The scores for module M were then regressed on the x- and y-coordinates of the centres of all sites in this neighbourhood. This allowed us to extract the equation of best fit for the scores for module M around each site s, that is, M_s = a_M,s*x + b_M,s*y + c_M,s. The coefficients of this equation were then used to define a scalar vector, ({overrightarrow{{{{{{{{rm{g}}}}}}}}}}_{M,s}=left[begin{array}{c}{a}_{M,s} {b}_{M,s}end{array}right]), whose directional component corresponded to the direction of steepest increase in the scores of module M around s, and whose magnitude component ∣g∣_M,s represents the steepness of that increase. Neighbourhoods in which all sites had the exact same score for module M had a magnitude of 0 and an undefined direction. Finally, the local vectors for all sites were combined to produce a two-dimensional gradient field ({overrightarrow{{{{{{{{rm{F}}}}}}}}}}_{M}) for module M across the entire nest. The gradient field typically took values of 0 at the core or far outside the module of interest, where module scores tended to take homogeneous values of 1 (module core) or of 0 (outside of the module). By contrast, the gradient field typically took non-zero values in areas of transition between adjacent modules, with the steepest values coinciding with the module’s borders (Fig. 8b).

Quantifying worker movement in the module gradient field

As worker movement may be affected by spatial heterogeneities and physical features within the nest, such as the presence of physical barriers like nest walls, we first established the typical movement of the average worker within the gradient field of each module. To do so, for each module M and each individual i, each trajectory segment was classified into one of two categories according to whether the individual was heading up-gradient or down-gradient within that module’s gradient field, ({overrightarrow{{{{{{{{rm{F}}}}}}}}}}_{M}). Segments were classified as up-gradient when the absolute angular difference between the trajectory heading and the gradient vector direction at that location was smaller than π/2 radians, and down-gradient when it was greater than π/2 radians.

To quantify the turning behaviour of individuals heading down-gradient in the field of their primary module (i.e., when i is heading out of M_i), we defined a ‘relative down-gradient turn angle’,

where ({bar{theta }}_{i,s,{M}_{i}}^{down}) is the mean turn angle of individual i when heading down-gradient in its primary module field, ({overrightarrow{{{{{{{{rm{F}}}}}}}}}}_{{M}_{i}}), at site s, and the second term, ({bar{theta }}_{s,{M}_{i}}^{down}) is the mean turn angle of all individuals heading down-gradient in ({overrightarrow{{{{{{{{rm{F}}}}}}}}}}_{{M}_{i}}) at site s. Positive values of ({{Delta }}{theta }_{i,s,M}^{{{{{{{{rm{down}}}}}}}}}) occur when a resident worker approaching the border of its primary module from the inside makes bigger turns than the average worker. Such workers tend to turn away from the border, and so are likely stay inside their primary module.

To quantify the behaviour of individuals heading up-gradient in their primary module (i.e., i heading into M_i), we also defined the ‘relative up-gradient turn angle’,

Negative values of ({{Delta }}{theta }_{i,s}^{{{{{{{{rm{up}}}}}}}}}) occur when a resident worker heading into M from the outside turns less than the average worker. Such workers tend not to turn away from the border, and so are likely to enter their primary module.

Statistical analyses

The proportion of specialist workers and of non-overlapping sites were analysed using general linear models (GLM) implemented using the package stats version 3.6.1 for R. The proportion of specialist workers was subjected to a square-root transformation to ensure normality of residuals (Shapiro–Wilk test, n = 50, proportion of specialist workers: W = 0.987, p = 0.84; proportion of non-overlapping sites: W = 0.980, p = 0.53). After fitting the GLM, the significance of the main effects was evaluated using F-tests.

All linear mixed-effects models (LME) and the Poisson generalized linear mixed-effects model (GLMM) were implemented using the package lme4 version 1.1–13⁶⁶ for R. In these models, continuous explanatory variables were scaled where necessary (e.g., in models including both colony size and individual-site similarity as explanatory variables, as these differed in scale by several orders of magnitude), though scaling did not affect the direction or significance of the main effects. To check that the LME model assumptions were not violated, we did not use traditional normality tests because those are not suitable for large sample sizes of more than 300 data points⁶⁷. Instead, we calculated the skewness and kurtosis of each model’s residuals and checked that they were compatible with a normal distribution (i.e., skewness between −2.1 and +2.1 and kurtosis <7.1;⁶⁷). Where necessary, dependent variables were subjected to a square root-, power- or log-transformation to ensure the normality of residuals. In the final models, skewness ranged from −2.1 to 0.5 and kurtosis from 1.6 to 6.6. After fitting the mixed-effects models, the significance of main effects was evaluated using Wald χ²-tests. All post-hoc comparisons were implemented by the package multcomp version 1.4-10⁶⁸ for R using the BH method to correct for multiple testing.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Source: Ecology - nature.com