Introduction

In many group living animals, individuals do not socialise indiscriminately; rather, they have preferred and avoided associates (Hinde 1976; Whitehead 2008; Strickland et al. 2017; Kappeler 2019). Preferred associates are conspecifics with whom individuals choose to associate, cooperate or otherwise synchronise their behaviour: they share a social bond. Here we refer to preferential social association between a pair of individuals as a social bond, and the sum of an individuals or a group of individuals’ social bonds as their social preferences. In animals, social bonds are usually inferred from observed patterns of association or interaction. Patterns of social bonds across a population or social group give rise to complex emergent social structures, and social network analysis has emerged as a powerful tool to study these emergent patterns (Krause et al. 2015; Croft et al. 2016). Changes in social preferences—for example in response to changes in group demography—will necessarily result in changes to the emergent structure of the social network. However, observed rates of association can be affected by extrinsic factors other than social preference, which means it is methodologically difficult to quantify how social preferences change with time or in response to external events. These methodological challenges contribute to the current gap in our understanding of how individuals change their social preferences in response to the turnover of individuals in the population (Ilany and Akçay 2016; Shizuka and Johnson 2019). Here we present a new method to quantify social preferences and apply this to understand how a highly social marine mammal changes their social preference in response to the birth and death of key individuals.

A key reason why comparing patterns of social bonds is challenging is because the observed rate of association or interaction between a pair of individuals can be affected by extrinsic factors other than social decisions. Behavioural differences between populations, or within the same population at different times, may result in pairs of individuals with the same underlying social preference having different observed rates of interaction. Individuals in highly social species, for example, are likely to have higher observed pairwise rates of interaction and association compared to those in less social species. These differences in observed association rate may occur not only between populations and species but also within subgroups of a population or within the same population at different times. Moreover, data collection methodology can lead to systematic differences in observed rates of association (Croft et al. 2008; Whitehead 2008). For example, data collected by focal follows are not always directly comparable to those collected from ad hoc observations of association (Davis et al. 2018). Indeed, a key challenge when trying to understand the social structure of an animal population is to control for these sampling effects and uncover the ‘true’ association patterns in the data (Bejder et al. 1998; Croft et al. 2011; Farine and Whitehead 2015; Weiss et al. 2020). The metrics of social preference between individuals are therefore determined not only by their social bond but also by how the data were collected and features of the animals’ social behaviour. To compare patterns of social preference between systems, it is necessary to take this variation into account (Croft et al. 2008; Whitehead 2008).

In a recent study, Weiss et al. (2019) used binomial mixture models to capture the variation in social associations within a population as a way to understand social complexity. Mixture models are—in essence—a form of cluster analysis. By identifying clusters of similar associations in observed data, the models can be used to transform a continuous value of social association strengths to a categorical variable, thus increasing the robustness of the measures. For example, in a population where individuals have a strong social bond with members of their social group but occasionally associate with members of other social groups, there are two categories of social association: within-group associations and extra-group associations. Mixture models offer a way to pick out these categories of social bond from observed patterns of association. Unlike the continuous social association rate, this categorical variable is robust to variation in sampling effort, social gregariousness and network completeness (Weiss et al. 2019). An unexploited yet potentially powerful application of these models is to use these categorised associations as a tool to compare social bonds between networks of different individuals or the same individuals at different time points. This opens the potential to examine what factors are important in driving the structure of social preferences and how this changes over time or how they differ between populations and species.

We illustrate our method with the example of demographic change in resident killer whales. Demographic change—the ageing, recruitment (birth or immigration) and death of individuals in a social group—has the potential to profoundly alter the structure of animal social networks (Shizuka and Johnson 2019). However, studies to date have rarely tested the effect of demographic processes on animal social structure, partly because of the methodological challenges of doing so (Shizuka and Johnson 2019). In many mammals, individuals of at least one sex associate with their mother throughout their life (Greenwood 1980; Mabry et al. 2013). How this mother-offspring relationship changes as the offspring ages is likely to be a key axis of variation in mammalian social systems (Prox and Farine 2020). In cetaceans, for example, the strength of mother-adult offspring association is argued to be linked to the degree to which the population is split into sub-groups (Rendell et al. 2019). Thus, understanding the ontogeny of mother-offspring relationships has the potential to give important insight into how social structures develop and are maintained in social mammals. Similarly, with lifetime mother-offspring associations the death of an individual’s mother represents a dramatic change in the offspring’s social environment. How individuals adjust their social preferences after the death of key social partners, such as their mother, is central to understanding the emergent structure of animal societies (Firth et al. 2017; Shizuka and Johnson 2019).

Resident killer whales (Orcinus orca) live in multi-level fission-fusion societies (Bigg et al. 1990). Killer whale society is based around matrilines: neither males nor females disperse and adults of both sexes are regularly found in close association with their mothers (Bigg et al. 1990). In turn, matrilines preferentially associate in pods forming a multi-level social structure (Bigg et al. 1990). The southern resident killer whales have been subject of detailed demographic and social study since 1976. However, inter-annual sampling differences in observed association rate complicate in-depth study of temporal changes in social preference. This long-term, multi-generational data is rare, especially in long-lived animals such as cetaceans. This population therefore offers an ideal system for illustrating the use of mixture models to investigate the interplay of social and demographic change. In addition, the absence of dispersal means that the only sources of individual turn-over are the death of whales, and the birth of new calves, simplifying the processes of demographic change. And as mating is usually outside of the immediate matriline, only a female offspring joins the close social group (Ford et al. 2011, 2018). Furthermore, in resident killer whale, mother-offspring relationships are a key determinant of survival in both adult males and females (Foster et al. 2012). In this study, we use our newly developed methodological framework to understand how mother-offspring relationships change as the offspring ages and as new offspring are born. We also use the framework to understand how the social environment of individuals whose mothers are alive differs from those whose mothers have died.

Overall, we have two aims with this study. Firstly, we aim to develop and apply a novel analytical technique to allow comparison of social preferences among times, groups or populations. Secondly, we use this methodology to investigate the links between demography and social structure in resident killer whales.

Methods

Binomial mixture model method

Binomial mixture models and social associations

We use binomial mixture models to cluster association indices into one, two or more components. In these models, the observed number of associations between a pair (i and j) of individuals (x: the numerator of a social association index) is considered to be the result of a sample from the number of times the members of the pair are observed (d: the denominator of the social association index) and their social bond (b):

$$ {x}_{i,j}\sim \mathrm{binomial}\left({d}_{i,j},{b}_{i,j}\right) $$

Binomial mixture models, in essence, fit n binomial distributions to the observed associations (see ‘How many components?’ section for an explanation of how n is chosen). Each distribution is referred to as a component (k). Maximum likelihood estimation is then used to select the most parsimonious mean and deviation for each component to explain the observed distribution of the data. More information on the use and application of binomial mixture models in the study of social associations can be found in Weiss et al. (2019) and references therein.

We apply this method to categorise social bonds as belonging to a particular component. This method explicitly assumes that there are a certain underlying number of rates of social association, and the number of these association rates does not vary between years. By assigning a social association to a component, therefore, we can investigate the following: (a) how the properties of the individuals associating predict the strength of social association they will have; and (b) what social and life-history traits predict when an association between a pair of individuals will change their association rate. Overall, this allows us to understand how the social structure of a population is linked to the properties of individuals in the population.

Binomial mixture models were developed, tested and applied in R (R Development Core Team 2019), with the dplyr, ggplot2 and VGAM packages (Wickham 2016, 2019; Pedersen et al. 2019; Wickham et al. 2019). We have also developed a package—socmixmods—to facilitate the implementation of mixture models in social association data. This package can be found at: github.com/samellisq/socmixmods.

Resident killer whales

The southern resident killer whales are a population of the salmon-eating ‘resident’ killer whale ecotype inhabiting the North East Pacific ocean (Bigg et al. 1990; Ford et al. 2000). They are a closed population, with no known genetic contact with other sympatric killer whale populations (Olesiuk et al. 2005; Ford et al. 2011). The southern resident killer whales have been the focus of a long-term individual based demographic study by the Center for Whale Research since 1976. All whales in the population are individually identified based on unique appearance of their saddle patch and dorsal fin, and, in adults, clear sexual dimorphism in fin shape and body size allows identification of the sexes of individuals (Bigg et al. 1990; Ford et al. 2000). Sub-adults can be sexed based on opportunistic sightings of ventral genital markings, but some individuals die before their sex can be ascertained (Bigg et al. 1990). As almost all whales in the population are observed annually and neither males nor females disperse from the population, all births (of calves surviving until they first observed which usually occurs within the first 6 months) and deaths in the population since 1976 can be reliably inferred from the census data. Our study is based on these data collected from live animals in the field and hence it was not possible to record data blind.

In the summer months, resident killer whales regularly enter the waters around the San Juan Islands, Washington State, USA. As part of intensive surveying efforts, the Center for Whale Research has collected video and photographs of groups of whales encountered in the waters around San Juan Island since 1976. We use these images to define social associations between whales. Whales are considered to be associating if they are photographed surfacing synchronously or successively within 1 body length of each other during an encounter. For the purposes of this study, we use the social data collected by the Center for Whale Research between 1990 and 2015 over which time the median population size was 86 (range: 78–98). During this time, there has been a consistent, intensive sampling effort to quantify patterns of group structure and social associations in those groups have been recorded. We applied a chain rule to define groups of associating whales per sampling day (further details can be found in Ellis et al. 2017). We use a daily sampling period so that two whales observed in the same group at any point during a single day are considered to have been part of the same group on that day. Application of the chain rule therefore rarely, if ever, has the effect of splitting a group of whales in close association into multiple groups within the sample period (1 day). For each pair of whales i and j in the population in a given year, we calculated the numerator and denominator of the simple ratio index (Cairns and Schwager 1987; Whitehead 2008). Specifically, we calculated the following: (1) the number of times they were observed in the same group in a day (simple ratio index: x, with a daily sampling period, see also ‘Binomial mixture models and social associations’ section above) and (2) the total number of times they were observed overall including when apart (simple ratio index d: x + yi+yj + yij). Where yij is the number of sampling periods (here days) in which whales were identified but did not associate (Whitehead 2008). We do not calculate the social association for whales in either their year of birth or year of death as the partial data available for those years may affect measures of social association. In this study, we use these statistics as the basis for classifying social associations into categories (see below; Weiss et al. 2019). Overall, between 1990 and 2015, 173 whales lived in the population of whom 83 were female and 78 were male (some calves die before they have been sexed). The median number of years lived through the dataset between 1990 and 2015 by females was 18 (range 0–25) and by males was 8.5 (range 0–25). This difference is likely to be driven by the strong sex difference in lifespan in the population: very few males reach the age of 40 whereas females often live into their 80s and beyond, with many females living to have a long post-reproductive lifespan (Olesiuk et al. 2005; Foote 2008; Croft et al. 2015; Ellis et al. 2018).

Maternal pedigree has been constructed for the whales in this population based on observations of infant swimming behaviour and population genetic structure (Bigg et al. 1990; Barrett-Lennard 2000). Detailed annual population censuses mean that mother identity and age are known for all whales born in this population since 1976 (Bigg et al. 1990). In addition, for calves who had not reached sexual maturity by 1976, mother identity and approximate age were inferred from patterns of association and the onset of sexual maturity (see Bigg et al. (1990) for more details). In these pre-1976 calves, mothers were only assigned if they were known with a good degree of certainty, and ambiguous cases were not included (Bigg et al. 1990). These known and inferred mothers were used to build a maternal pedigree and to assess the relation between pairs of individuals in the population (see pseudocode workflow supplementary 1). We establish relatedness to the second degree (i.e. where coefficient of relationship r ≥ 0.25: grandparent-grandoffspring, aunt/uncle-nephew/niece), as our resident pedigree is not deep enough to regularly establish more distant relationships. We classify dyads by kinship class (sister-sister, aunt-niece, grandmother-grandson etc.) rather than relatedness (e.g. r = 0.5) because we are interested in the social association between kinship classes rather than the impact of genetic relatedness on social behaviour per se. We calculated kinship class conservatively: if the relationship between a pair of individuals cannot be established, it is marked as unknown. For example, if one of a pair of whales do not have a known mother, it cannot be established if they are siblings, so we class their kinship class as ‘unknown’ (if there is no other relationship such as mother-offspring; supplementary 1). Eighty-six percent of dyads have unknown kinship in our data and are therefore not included our kinship analysis. Matrilines in this population form long-term stable social groups called ‘pods’ (Bigg et al. 1990; Ford et al. 2000). Pods have been stable over the length of the 40+ years of this study with no observed pod fusion or movement of matrilines between pods (Bigg et al. 1990; Parsons et al. 2009), and analysis of unique within-pod vocalisations suggests that pods represent groups of common matrilineal descendants (Ford 1989, 1991). We therefore assume all dyads between whales from different pods are not maternally related and classify them as ‘out-pod non-relatives’. Dyads are classified as in-pod non-relatives if they can be demonstrated to not be related within the second-degree (i.e. r is not ≥0.25; see supplementary 1).

How many components?

To assess how many components are present in the data, we first needed to establish a criterion to determine model fit. Weiss et al. (2019) found the integrated completed likelihood (ICL) to perform best as model-fitting criteria for association data; however, they found that if the number of separate sighting (index denominator) was less than 40, the ICL could sometimes underestimate the number of components. The mean number (± std. dev. [SD]) of sightings per dyad (the simple ratio index denominator) for the killer whale data is 24 ±8. We therefore repeated the simulations used in Weiss et al. (2019) but with killer whale-like input to determine the best performing model fitting criteria for this data. We simulated a population with a known number of components (k = 1–9), with a population size chosen from a normal distribution derived from the observed killer whale data (mean = 81, SD = 4), and a sampling effort (simple ratio index denominator) chosen from normally distributed representation of the observed sampling efforts over all years of the study (mean = 24, SD = 8). We ran this simulation 500 times and scored the number of times the three tested model fit criteria—Akaike Information Criterion (AIC), Bayesian Information Criterion (BIC) and ICL—correctly assigned the best fitting model to be the inputted k value. BIC had the highest success rate, assigning the correct number of components in 293 simulation iterations (AIC success = 202, ICL successes = 271). We therefore used BIC as our model fitting criteria. We suggest others using this technique run similar simulations to assess both the utility of fitting criteria and their accuracy given the data under consideration

The aim of our analysis is to understand how social associations change in response to demographic events. To do this, we assume that the killer whale social system has a fixed underlying social structure and a corresponding fixed number of components, which does not vary between years (see ‘Discussion’ section). We model this underlying structure by fitting a model of the same number of components to each year of the data. Alternative analytical frameworks, such a fitting different numbers of components to each year of the data, may be appropriate to answer other research questions. To assess the most parsimonious number of components over all years, we fitted models with 1 to 9 components to each year of association data. The fit of each component model (BICk) was assessed by BIC. Within each year, we calculated the deviation (ΔzBICk) of each component models z-scored BIC from best fitting component models z-scored BIC:

$$ \Delta z{\mathrm{BIC}}_k=\left(\frac{{\mathrm{BIC}}_k-{\nu}_{\mathrm{BIC}}\ }{\sigma_{\mathrm{BIC}}^2}\right)-\underset{i\in K,\kern0.75em K=1,\dots, 9\ }{\min}\left\{\frac{{\mathrm{BIC}}_i-{\nu}_{\mathrm{BIC}}}{\sigma_{\mathrm{BIC}}^2}\right\} $$

where νBIC and σ2BIC are the mean and standard deviation of the BIC over all components (1–9). The best fitting model in each year will therefore have a ΔzBIC of 0, and a model fit one standard deviation from the best fitting model would have a ΔzBIC of 1. We then summed these deviations over all years to assess the parsimony of each component model over all years (Table 1). This process demonstrated that a three-component model is the most parsimonious model over the 26 years of the study (see ‘Results’ section for details; also see four-component model supplementary 2). We recommend other users of the method similarly explore the most parsimonious number of components their model contains and explore the biological implications of alternative models (supplementary 2).

Table 1 Mean numbers of intermediate k2 and strong k3 social connections in male and female killer whales in different age categories. Age categories are derived from the literature (Olesiuk et al. 2005). Female age at adulthood reflects the average age of first birth of the first offspring, whereas males represent the average age the physical maturity is attained (as opposed to sexual maturity which achieved several years earlier) (Olesiuk et al. 2005). Age of becoming post-reproductive is the age at which 95% of population female fecundity has been achieved (Ellis et al. 2018). Separating physically mature males from sexually mature males (Olesiuk et al. 2005) does not make a qualitative difference to the results (supplementary Table 6)
Full size table

Assigning bonds to components

We applied a three-component mixture model to each year of association data in our population. Here we use the same technique to assign every social bond in the population to a component. During the fitting of the mixture model, each association is assigned a probability μ of belonging to each fitted component. We consider a given social association to belong to the component which has the highest μ. Component assigning was rarely ambiguous: for 80% of associations, the difference between the highest μ and second highest μ was over 0.8, and in only 1.2% of associations was the difference in μ less than 0.1. Components, k, are numbered from 1 to 3, with k1 as indicating an association in the weakest component, k2 in the intermediate strength component and k3 as an association in the strongest component.

Applying the method

Details of the statistical analyses used in each section of the results are detailed below. For all analyses, we report posterior mean and the 95% credible intervals around these coefficients for each reported factor. Credible intervals not overlapping 0 give 95% confidence that the coefficient in question is different from 0. All statistical analyses described below were performed in R (R Development Core Team 2019), using the dplyr, stringr, ggplot2 and brms packages (Wickham 2016, 2019; Bürkner 2017; Wickham et al. 2019). The statistical models described below were evaluated with Hamilton Markov Chain Monte Carlo algorithm implemented in stan via brms (Bürkner 2017; Stan Development Team 2020). All models were checked for, and adjusted to ensure: convergence of coefficients; absence of divergent transitions; sufficient tree depths; good mixing of chains; and effective sample sizes in both the ‘bulk’ and the tails of the distribution.

We explored the social environment of individual killer whales using Bayesian generalised linear mixed effects models to test for differences in the number of associations in each component maintained by whales in different age-sex classes For each model, the number of social associations of a given component (k =1–3) a whale has was used as the response variable, with either sex or age category as the predictors. Both males and females had infant, juvenile and adult age categories, and females had an additional post-reproductive age category (Table 1). Where sex was the predictor, we used a random intercept for the age category. All models had a Poisson error-structure with a log link function, and whale id and year as additional random intercepts.

Social complexity is not a focus of this manuscript but we report Shannon’s entropy as a measure of social complexity (Weiss et al. 2019) for each annual model to highlight other uses of the model framework and to allow comparison with other systems.

Social dynamics of mother-offspring associations

We use the fitted components to explore how mother-offspring associations change as the offspring ages. Specifically, we test how the probability of a mother-offspring association being in the strongest k3 component changes with offspring age for males and females. We model the probability (P(k3)) that the relationship between a mother (m) and offspring (o) in a given year is in (P(k3)o,m) using a Bayesian generalised mixed effects models with the basic structure:

$$ {\displaystyle \begin{array}{c}\mathrm{k}3o,m\sim \mathrm{Bernoulli}\left(P\left(\mathrm{k}{3}_{o,m}\right)\right)\\ {}\mathrm{logit}\left(P\left(\mathrm{k}{3}_{o,m}\right)\right)={\mathrm{sex}}_o+{f}_1\left({\mathrm{age}}_o\bullet I\Big({\mathrm{sex}}_o=0\right)\Big)+{f}_2\left({\mathrm{age}}_o\bullet I\left({\mathrm{sex}}_o=1\right)\right)+{\varepsilon}_{o,m}\\ {}\ \varepsilon \sim N\left(0,\sigma \right)\end{array}} $$

where I is an indicator function that takes the value 1 if its argument is true, and 0 otherwise, f1 and f2 are estimated smooth functions, ageo is the offspring’s age (in years) and sexo is the offspring’s sex (0 = female, 1 = male). The effect of offspring age was modelled as a smooth function of age because we are interested in understanding potentially non-linear patterns between offspring age and P(k3). We estimated these smooth functions using penalised thin-plate splines as the basis, with the number of knots chosen by generalised cross-validation (Bürkner 2017; Pedersen et al. 2019). The term εo, m is a random intercept (group effect) for each mother-offspring pair, included to account for repeated observations of the same dyad. Models used uninformative priors for all parameters. We used the same structure to model the probability that a mother-offspring association was in component 2 (P(k2)), and the probability that an association is in either k2 or k3 (P(k3 + k2)). The coefficients of fitted smooths cannot be easily interpreted; we therefore describe the output of the model graphically, and using the predicted means and credible intervals of the relationship between component probability and age.

We also explore how mother-offspring associations change in response to social events, namely, (1) the birth of an offspring to the daughter in the mother-daughter association (daughters’ offspring model); (2) the birth of an offspring to the mother in the mother-offspring association (mother’s offspring model). To examine the effect of either of these events, we compare the social preferences of the actors before and after the event. We focus on the 5 years either side of the event, not including the year of the event. Five years (mean ± SD: 5.48 ± 2.58) is the mean inter-birth interval in this population (derived from (Nattrass et al. 2019)), and results do not qualitatively differ with small changes in this threshold (SE unpublished results). The data are treated as categorical with the years before the event coded as 0, and the years after the event coded as 1. It is important to note that the analysis is therefore limited to the associations where the event occurs: so in the daughters’ offspring model, mother-daughter pairs where the daughter does not produce her own offspring by 2015; and in the mother’s offspring model, offspring who do not have a younger sibling by 2015 are not included in the model. As in the previous models, the response variable is whether a particular mother-offspring association in a given year is in component k3. The predictors indicate whether the year is before or after the event. Both models have dyad-id and the age (in years) of the offspring (in the mother-offspring association) in the year of the event as random intercepts. The mother’s offspring model has mothers age (in years) as an additional random intercept. The models have a Bernoulli error structure and logit link function.

Impact of mother death on social preferences

Classifying associations into components provides the opportunity to investigate how social associations with other kinship classes differ when a whales’ mother is dead compared to when they are alive. This analysis uses all whales with a known mother in the population, classifying in any given year if their mother is alive or dead. A smaller sample size precluded the use of the before-after framework we used for the changing mother-offspring association analysis described above. Computational and conceptual limitations lead us to model the association with each kinship class separately. In each model, the probability of the association with a given kinship class being in component 3 (P(k3)) is modelled as the response variable in a Bayesian generalised mixed effects model, with mother alive as a predictor. Offspring age (in years) and offspring id are included as random intercepts.

We also included offspring sex as a predictor interacting with mother status if the number of observations in all categories (mother alive or dead vs offspring sex male or female) was greater than 50 (supplementary Table 2). Fifty was arbitrarily chosen as the threshold but corresponds approximately to the number of observations that are required for these models to allow a model with an interaction term to fit to the data. Offspring sex was included in models of: same-sex siblings; aunts, in-pod non-relatives and out-pod non-relatives. If there were fewer than 50 observations in one or more sex-status categories (supplementary Table 2), the models did not include sex as an interacting effect. It is also important to note that some matrilineal kinship classes—namely offspring and grandoffspring—are only possible for females. Small sample sizes precluded modelling of the relationships between a whale when their mother was alive or dead and their uncle and grandoffspring (supplementary Table 2); these kinship classes are therefore not included in the results. A lack of variation in the response values (no k3 associations are between whales in different pods) prevented the out-pod non-relative response to mother death being modelled, so this kinship class is not included in this analysis.

Results

Social mixture model summary

Three-component mixture models are the most parsimonious fit to the data over the 26 years of the study (summed ΔzBIC of 0.74; Table 1). In support of this interpretation, the three-component mixture model was the most commonly the best fitting model (11/26 years). The one-component mixture model was the least parsimonious model and was not the best fitting model in any year of the study. The four-component model was the second most parsimonious and second most commonly best fitting model (supplementary Table 1). The results presented below are qualitatively similar when the analysis presented below is repeated with a four-component model (supplementary 2).

We fitted three-component binomial mixture models to each year of 26 years of resident killer whale association data, consisting of a total of 86343 social (or potentially social) associations between 173 whales. Of these associations, 77562 were classified as weak k1 component associations (this includes both weak and completely absent associations), 7338 were classified as intermediate strength k2 component associations and 1443 were classified as strong k3 component associations (supplementary Fig. 1; supplementary 3). The median social complexity in the three-component models was 1.06 (range 0.69–1.10; supplementary Table 3).

In a given year, each whale maintains, on average, 1.90 (±1.10, mean ± SD) of the strongest k3 associations and 7.10 (± 5.34) intermediate strength k2 associations. All other social connections are of the weakest k1 level. Males do not differ from females in their number of k2 or k3 associations while either immature or adults (supplementary Table 4). However, whales of different sex-age categories differ in the number of k2 and k3 associations they maintain (Table 1). In particular, males have fewer k2 associations as infants than as adults, and females have fewer k2 associations as infants and juveniles than as adults (supplementary Table 5).

Kinship class is known for 12541 (14%) annual whale-pair dyads, representing 2917 unique pairs of individuals (because many pairs are observed over multiple years). We found a strong link between kinship and social association in resident killer whales: closer relatives tend to have associations classified in a higher component (Fig. 1). For example, 44.6% of mother-offspring pairs, and 71.4% of mother-immature offspring pairs, are found in the strongest k3 category (Fig. 1), while 90.1% of associations between non-relatives are found in the weakest, k1 association category (Fig. 1). Under the alternative four-component model, the additional category separates the mother-immature offspring relationships into an ‘extra strong’ association component (supplementary 2).

Fig. 1
figure 1

The percentage of each category of relatedness pair classified into each component (k1 component = weakest category of association; k3 component = strongest category of association) of a binomial mixture model from 26 years of social association data. y-axis represents the categories of matrilineal kinship distinguishable from the pedigree up to (matrilineal) relatedness of 0.25, with sample size in brackets. Mother-offspring relationships are split into those between immature (infant and juvenile) offspring and their mother and adult offspring and their mother. All other kinship categories include all ages of partners. Numbers in each tile show the percentage of a given category of social bonds classified in that component. Tile colour also represents the percentage: from a low percentage of associations in that component (white/pale yellow) to a high percentage of associations in that component (dark blue). This figure represents only pairs of individuals with a known kinship, and ‘non-relatives’ are individuals known to have a relatedness of <0.25

Full size image

Social dynamics of mother-offspring associations

Between 1990 and 2015, we observed 1146 mother-offspring association years from 201 unique dyads. 96 of the associations (covering 588 association years) are between a mother and her male offspring, and the remainder are between mothers and their female offspring. The mean (± SD) age of sons in the data is 9 (± 6) while the mean age of daughters is 13 (± 9). Mother-offspring relationships change as offspring age, and the change in this relationship differs for male and female offspring (Fig. 2). Predicted probability of mother-daughter association being k3 drops sharply between the ages of 12 and 18 from 0.84 (post. mean, cred. int. = 0.73–0.92) to 0.40 (post. mean, cred. int. = 0.29–0.52). In contrast, the probability of mother son-association being in k3 does not change between 1 (post. mean = 0.72, cred. int. = 0.60–0.82) and 26 (post. mean = 0.59, cred. int. = 0.19–0.84). The probability of mother-offspring bonds being in k2 mirrors the probability of the association being in the k3: there is no change in the probability of mother-son’s bonds being in the k2 with age, and a sudden increase in the probability of mother-daughter associations being in the k2 between ages 12 and 17 (supplementary Fig. 2). There is no change in the probability of mother-offspring associations being in either the k2 or k3 (i.e. not the k1) with age for either sex (supplementary Fig. 3). Biologically similar dynamics are found in the four-component model with the extra step of associations moving from the strongest ‘mother-immature offspring’ association category to the next strongest category over the first 10 years of the offspring’s life (supplementary 2).

Fig. 2
figure 2

The posterior predicted (± credible intervals) probability of a mother-offspring association being in component 3 as a function of age for male and female offspring. Predictions are from a Bayesian generalised linear model with offspring age modelled as a four-knot spline. Lack of overlap in credible intervals between points suggests that the whales with those properties have different probabilities of their association being in k3

Full size image

The daughters in 53 of the 105 mother-daughter dyads gave birth to their own offspring. In the 5 years after the birth of their first offspring females (daughters’ offspring model, see ‘Methods’ section) have a lower probability of their association with their own mother being in component 3 than they have in the 5 years prior to the birth of their first offspring (post mean = −1.17, cred. int. = −2.06 to −0.38; Fig. 3). In contrast, the 5 years either side of the mother giving birth (mothers’ offspring model, see ‘Methods’ section), there is no difference in the probability of a mother-daughter association being in k3 (post mean = −0.44, cred. int. = −1.20–0.31; Fig. 3). In both models, offspring age (as year) is included as a random intercept.

Fig. 3
figure 3

Number of mother-daughter associations in k3 and not in k3 in the 5 years before and after the birth of a new offspring (n = 105). Left tile shows the 5 years before and after the birth of an offspring to the daughter, and the right tile shows the 5 years either side of the birth of an offspring to the mother

Full size image

Impact of mother death on social preferences

The predicted probability of the association between siblings of opposite sex being in component 3 is higher when their mother is dead than when their mother is alive (post. mean. 0.99, cred. int. = 0.44–1.52; Fig. 4). Similarly, grandmother-grandoffspring associations have a higher probability of being in component 3 when the grandoffspring mother is dead than when she is alive (post. mean. = 3.11, cred. int. = 1.18–5.98; Fig. 4). The probability of a whale’s social association with other relatives, including same sex siblings, and non-relatives does not change when the whale’s mother dies (Table 2; Fig. 4).

Fig. 4
figure 4

Predicted change in the probability of the association between a whale and their relations when the whale’s mother is alive and when she is dead being in the strongest social component (k3). Predictions are from a Bayesian linear model, described in the text. Each panel represents a separate model. In each panel points shows the posterior mean and error bars the 95% credible intervals. An interaction between offspring sex and mother status is included if each of the four categories (mother alive, dead vs sex: male, female) have more than 50 data points, if not sex is not included in the model. Note, as all kinship classes are matrilineal only females have offspring as relations. The probability of the association between opposite sex siblings and grandmothers and their grandoffspring being in component 3 is if a whale’s mother is dead. All other relationships are unchanged. Sample sizes were too small to investigate how relationships between a whale and their uncles or grandoffspring changed after the death of their mother, and a lack of variation in the response prevented the out-pod non-relative model fitting, so these relationships were not modelled and not included in this figure

Full size image
Table 2 Output of models investigating how the probability of a whale’s social association with a given ‘Relation’ changes after the death of the whale’s mother. In all models, probability of an association being in k3 is used as the response variable with ‘Variable’ as the fixed effects in a Bayesian linear model, with a Bernoulli error structure and offspring age and offspring identity as random effects. Offspring sex and an interaction between offspring sex and mother status are included if there are more than 50 data points in all four possible categories (mother status vs. sex); otherwise, sex is not included in the model. Whales of unknown sex are not included in the analysis. All other relationships are unchanged. Note, as all kinship classes are matrilineal, only females have offspring as relations. All other relationships are unchanged. Sample sizes were too small to investigate how relationships between a whale and their uncles or grandoffspring changed after the death of their mother, and a lack of variation in the response prevented the out-pod non-relative model fitting, so these relationships were not modelled and not included in this table. Estimates with a credible interval not overlapping 0 are italicised
Full size table

Discussion

Comparing social preferences between different groups of individuals represents a significant methodological challenge. Here we have developed the use of mixture models to compare social preferences over time or between different groups of individuals. By classifying associations using mixture models, much of the noise derived from sampling variability and other extrinsic factors can be eliminated allowing us to compare across networks through time (or across populations) using relatively simple models and metrics. Previous work comparing social preferences and structures between networks has tended to rely on data being collected in a comparable way or by using statistical and randomisation methods to control for known drivers on inter-network differences (e.g. Matsuda et al. 2015; Meise et al. 2019; Wilkinson et al. 2019). The advantage of the methodology demonstrated here for comparing social preferences between networks is it transforms the continuous social association index—the values of which are susceptible to various extrinsic factors—to a categorical variable. This categorical variable allows us to develop a comparable measure of social preference independent of variation in association metrics driven by methodological and sampling effects. These categories then allow an investigation into how properties of individuals, such as their ages or life-history stage, predict properties of their association, and how that changes with time. In addition, the association categories are directly comparable between networks and can therefore be used to compare social preferences in different contexts and systems. We have illustrated the utility of this methodology be comparing over annual social networks in resident killer whales which has led to four important insights into their social system: (1) resident killer whales society can be considered as consisting of multiple classes of social bond defined by their association frequency; (2) the classes generally reflect the kinship structure of the population; (3) mother-son associations stay constant with age into adulthood, whereas mother-daughter associations decline after the birth of the daughters first calf; and (4) whales have higher association with their opposites sex siblings and their grandmothers if their mother is dead compared to if their mother is alive.

Encouragingly, the binomial mixture model method has provided results that concord with previous work in the southern resident killer whales. For example, in this study, we found that southern resident killer whales have three categories of association, and that the categories are associated with kinship. In agreement with this, using long-term observations and association indices, Bigg et al. (1990) proposed that the resident killer whale society has three tiers: matrilines, pods and communities. Further observations of patterns of association, analysis of genetic relatedness and study of vocalisations have largely confirmed this structure in the southern residents (Ford 1991; Barrett-Lennard 2000; Parsons et al. 2009), and this schema has been widely adopted by researchers studying this population (e.g. (Baird 2000; Olesiuk et al. 2005; Ford 2014)).This study has shown that the strongest associations are between close maternal relatives, the matriline, the second strongest bonds are between more distant relatives and some other members of their pod, with the weakest associations between whales in different pods. In our analysis, we assume (after analysis) that the population has an underlying structure with three types of social association. Changing this assumption to four types of association does not change our conclusions about killer whale biology (supplementary 2). In common with previous work on a non-overlapping dataset (Bigg et al. 1990) and anecdotal knowledge—using our new methodology, we found that while mother-daughter associations in resident killer whales decrease in strength as the daughter ages, mother-son associations remain strong throughout the sons life.

Our new methods allowed us to examine, for the first time, how social preferences change in response to demographic processes in resident killer whales. While males associate strongly with their mother their whole life, females associate less with their mother after the birth of their first offspring, whereas after the birth of a sibling, there is no change in the strength of relation between mothers and their daughters. This suggests that it is the daughters driving the change in mother-daughter association patterns with age rather than mothers. This accords with theoretical expectations. In species with bisexual philopatry and non-local mating—like resident killer whales—reproductive competition is expected between mothers and their daughters who are competing for the same resources (modelled in Johnstone and Cant (2010) and tested in Croft et al. (2017)). Competition for reproductive resources is not present between mothers and their sons because the sons offspring are outside the social group (Johnstone and Cant 2010; Foster et al. 2012). The changing social preferences of mothers and daughters could be part of a response to this reproductive competition. Interestingly, this contrasts with patterns found in a study of common bottlenose dolphins (Tusiops truncatus)—who have a different social structure—where previous work has suggested that mother-offspring associations decline sharply when their mothers become pregnant (Connor et al. 2000). More research is needed to understand the behavioural ecology of mother-daughter associations in cetaceans. A consequence of females associating less with their own mother is that they begin to separate their matriline from their mother’s matriline, which in turn will drive the development of a multi-level society. Our results support the argument that degree of matrilinearity is a key driver of the emergent structure observed in cetacean societies, and that increasing matrilineal affiliation can lead to increasing separation of the society into distinct social modules (Rendell et al. 2019).

The presence or absence of an individual’s mother is reflected in differing social preferences for other relatives. Whales whose mothers are dead have a higher probability of associating with their opposite sex siblings and their grandmother. Males are therefore associating more with their sisters (and vice versa). This suggests that at least part of the social support is being provided by close female kin. Reproduction in this population is dominated by the oldest males (Ford et al. 2011, 2018). By increasing their association with their brothers, females may also increase their inclusive fitness if an increased association with siblings helps their brother to survive (for example by increasing their social centrality (Ellis et al. 2017)). Grandmothers may also gain inclusive fitness benefits by increasing their association with their grandoffspring (we could not distinguish between association with grandsons and granddaughters in our analysis). Female resident killer whales have an extended post-reproductive lifespan, ceasing reproducing in their late thirties (Ellis et al. 2018). The evolution of this unusual phenomena is driven by the benefits that older females can provide their offspring and grandoffspring (Foster et al. 2012; Brent et al. 2015; Nattrass et al. 2019) and reproductive conflict between older females and their reproducing daughters (Johnstone and Cant 2010; Croft et al. 2017). By associating with their grandoffspring, grandmothers can increase their inclusive fitness both by increasing the reproductive success of their granddaughters and helping their grandsons survive long enough to reproduce (Hawkes et al. 1998; Ford et al. 2011, 2018; Nattrass et al. 2019).

These insights into the social structure of killer whales are driven by our novel methodology allowing us to classify social associations. Classifying social associations into components removes inter-annual variation due to sampling and reveals the underlying social structure of the system. The method of classifying bonds also potentially has other uses in this and other systems. For example, in this study, we were able to use the classified associations to quantify the social environment experienced by individual killer whales. Variation in the number of social associations of different strengths within a population, or between populations, may represent differences in cognitive capacity (Taborsky and Oliveira 2012; Dunbar and Shultz 2017) or a trade-off between socialising and other activities (Dunbar et al. 2009; Marshall et al. 2012) or lack of social opportunities (e.g. Brent et al. 2017; Goldenberg and Wittemyer 2017). Understanding the variation in individual social environments will give important insights into the links between cognition and social behaviour. Another potential use of this method is to identify types of association at the population level. Recently, there has been interest in understanding the role that ‘strong’ social bonds have on individual fitness (Silk et al. 2009, 2010; McFarland et al. 2017; Ellis et al. 2019). However, this approach has been criticised because it requires choosing either a threshold over which to consider an association ‘strong’, and changing this threshold can change the interpretation of the underlying biology (Ellis et al. 2019). Rather than choosing a threshold, applying mixture models would allow strong bonds to emerge from the data. This methodology may also provide a useful way to understand hierarchical or multilevel societies. Distinguishing the levels of social organisation in a multi-level society is a challenging analytical problem (Morrison et al. 2019; Grueter et al. 2020). Classifying associations may provide a new and general way in which to understand social association in a complex society. Lastly, classifying associations allows networks to be compared. Although here we compare the same population at different times, the same methods could be used to compare different populations, species and systems.

There are limitations to the methodology we describe here. Firstly, applying a mixture model to social data relies on several assumptions. Perhaps the most important assumption is that the diversity of social relationships in the population can be approximated to a fixed number of distinct association rates. This is a common assumption in studies of animal sociality, particularly in studies of social complexity, and is based around an understanding of cognition (Bergman and Beehner 2015; Fischer et al. 2017; Weiss et al. 2019). In a multi-level society—such as that of resident killer whales—the binomial components can be linked clearly to biology but in other types of social system, the link may be less clear. Secondly, we assume number of components or types of social relationship is constant through time. In biological terms, we are assuming that the system has a constant underlying social structure. This is not to say that there can be no annual variation from this structure: changes in demography (for example a year which happened to contain no mother-offspring pairs) or ecology (for example high levels of structure in good ecological conditions) may add or remove components in a given year. In this study, we use information criteria to choose the most parsimonious number of components to identify the core of the social structure which stays constant through the study period. However, in other systems, perhaps with more fluid social structures or more ecologically driven social variation, the assumption of an unchanging core of social structure may not be a valid. To answer other research questions, it might be more appropriate to make other assumptions or other analytical decisions. For example, to identify ‘strong bonds’ in a population, it might be more appropriate to combine associations for all years and apply the mixture model to this combined dataset. A third limitation of this approach is its complexity. An alternative approach to answering the questions addressed here is to control for effects of sampling either within the structure of a statistical model or using forms of data permutations (e.g. (Matsuda et al. 2015; Meise et al. 2019; Wilkinson et al. 2019)). However, these approaches also have limitations, particularly in being able to robustly quantify and account for sources of variation in the data. Overall, we have demonstrated that this method can be used to answer important questions and that, alongside other methods, can become a tool for researchers trying to understand the structure and function of animal social relationships.

In the study, we have demonstrated that using mixture models to classify social associations can lead to important insights into the underlying social structure in complex animal societies. We used this methodology to show the links between kinship structure and social structure in resident killer whale societies, and to show how theses associations change as individuals age, give birth and their relatives die. We expect that applying this methodology to other systems and other research questions will lead to further important insights into animal sociality.