Age- and sex-dependent variation in relatedness corresponds to reproductive skew, territory inheritance, and workload in cooperatively breeding cichlids.

Kin selection plays a major role in the evolution of cooperative systems. However, many social species exhibit complex within-group relatedness structures, where kin selection alone cannot explain the occurrence of cooperative behavior. Understanding such social structures is crucial to elucidate the evolution and maintenance of multi-layered cooperative societies. In lamprologine cichlids, intragroup relatedness seems to correlate positively with reproductive skew, suggesting that in this clade dominants tend to provide reproductive concessions to unrelated subordinates to secure their participation in brood care. We investigate how patterns of within-group relatedness covary with direct and indirect fitness benefits of cooperation in a highly social vertebrate, the cooperatively breeding, polygynous lamprologine cichlid Neolamprologus savoryi. Behavioral and genetic data from 43 groups containing 578 individuals show that groups are socially and genetically structured into subgroups. About 17% of group members were unrelated immigrants, and average relatedness between breeders and brood care helpers declined with helper age due to group membership dynamics. Hence the relative importance of direct and indirect fitness benefits of cooperation depends on helper age. Our findings highlight how both direct and indirect fitness benefits of cooperation and group membership can select for cooperative behavior in societies comprising complex social and relatedness structures.

Cooperative breeding, where individuals other than breeders provide alloparental care, evolved across many different animal taxa, including insects (Boomsma 2009; Biedermann and Taborsky 2011), crustaceans (Hultgren and Duffy 2012), fishes (Heg and Bachar 2006; Taborsky 2016; Tanaka et al. 2018b), birds (Koenig and Dickinson 2016), and mammals (Solomon and French 1997; Clutton‐Brock 2016). According to Hamilton's rule, selection may favor such alloparental care through both direct or indirect fitness benefits to brood care helpers, depending on their costs and benefits, and the relatedness between actors and beneficiaries (Hamilton 1964; Taborsky 1984; Clutton‐Brock 2002; Stiver et al. 2005; Field and Leadbeater 2016; Komdeur et al. 2017). Kin selection may explain why individuals delay dispersal and help related group members in raising their offspring (Hamilton 1964; Foster et al. 2006; Bourke 2014). Indeed, in several taxa, individuals seem to preferentially cooperate with close relatives (Choe and Crespi 1997; Russell and Hatchwell 2001; Griffin and West 2003; Koenig and Dickinson 2016). However, in many cooperative breeders at least some helpers are not related to the breeding pair producing the offspring to be cared for (Dierkes et al. 2005; Clutton‐Brock 2009; Riehl 2013). These helpers will not gain indirect fitness benefits but are assumed to acquire direct benefits instead, for example, through increased tolerance by dominant individuals allowing them to remain in the group (“pay‐to‐stay”, Gaston 1978; Kokko et al. 2002; Bergmüller and Taborsky 2005; Fischer et al. 2014; Kingma 2017; Naef and Taborsky 2020a, 2020b), which reflects an exchange of different commodities (Quiñones et al. 2016; Taborsky 2016). Continued group membership may benefit helpers, for example, due to reduced mortality risk (Heg et al. 2004, Heg et al. 2005b; Bergmüller et al. 2005; Kingma et al. 2014), shared reproduction (Richardson et al. 2002; Bruintjes et al. 2011; Riehl 2013; Hellmann et al. 2015), or the opportunity to inherit the breeding position in the future (Balshine‐Earn et al. 1998; Stiver et al. 2006; Riehl 2013; Field and Leadbeater 2016; Kingma 2017).

Clarifying relatedness patterns within groups and the interplay of potential direct and indirect fitness benefits of helping in complex animal societies is crucial for a proper understanding of the evolution of apparently altruistic helping behavior. This is a nontrivial challenge, because the received indirect fitness benefits may differ between helpers, depending on the variation in relatedness with the respective receivers of help (e.g., Dunn et al. 1995; Riehl 2013). In contrast, potential direct fitness benefits can arise for related and unrelated individuals alike, which complicates drawing straightforward conclusions about the relative importance of direct and indirect fitness effects of cooperation (Zöttl et al. 2013; Carter et al. 2019).

In fishes, cooperative breeding has been described for approximately 25 lamprologine cichlids endemic to Lake Tanganyika (Taborsky 1994; Heg and Bachar 2006), where it evolved several times independently (Dey et al. 2017; Tanaka et al. 2018b; Ronco et al. 2021). Notably, the cooperatively breeding species in this clade vary greatly in within‐group relatedness levels, ranging from species where most helpers are unrelated to the breeders they support, to others where subordinates usually help their own parents (Awata et al. 2005; Dierkes et al. 2005; Tanaka et al. 2015). For example, large helpers of N. pulcher are often unrelated to the breeders they aid (Dierkes et al. 2005; Stiver et al. 2005), whereas helpers of Neolamprologus obscurus are typically closely related to the breeders they assist (Tanaka et al. 2015). Hence, in N. pulcher indirect fitness benefits of helping are apparently much less important than direct fitness benefits (Jungwirth and Taborsky 2015), and they decline with helper age (reviewed in Taborsky 2016). In contrast, cooperative behavior of subordinates in N. obscurus is probably to a larger extent driven by kin selection (Tanaka et al. 2015). This striking divergence of selection mechanisms responsible for apparently altruistic alloparental care among closely related species sharing a common ecology provides unique opportunities to elucidate the significance of relatedness and group structure for the evolution of cooperative behavior in animals. In comparison to mammals and birds, lamprologine cichlids entail the additional benefit of methodological accessibility, as they have rather small home ranges (Heg et al. 2008; Tanaka et al. 2015; Jordan et al. 2016; Josi et al. 2019) and can be observed from a short distance with little disturbance (e.g., Taborsky and Limberger 1981; Maan and Taborsky 2008; Groenewoud et al. 2016; Josi et al. 2020b). Furthermore, most of the cooperatively breeding Lamprologini defend territories throughout the year and therefore can be observed during the reproductive as well as non‐reproductive periods (Brouwer et al. 2005; Josi et al. 2020b). This enables the collection of large sample sizes on social structure, relatedness patterns, and behavior within manageable time.

Here, we studied within and between‐group relatedness, growth, reproductive skew, and workload in the highly social cichlid Neolamprologus savoryi. This species breeds in complex groups with breeder males monopolizing one to several breeder females. Each breeder female defends a separate subterritory and may be assisted by subordinate helpers (together referred to as “subgroup”) within the male territory (Heg et al. 2005a; Garvy et al. 2015; Josi et al. 2019, Josi et al. 2020a,b). We sampled breeding groups in two populations and used microsatellite markers to estimate relatedness. Specifically, we asked whether (1) helpers gain indirect fitness benefits by living and helping in kin structured groups and subgroups; (2) helpers gain direct benefits through reproductive share or the chance to inherit the breeding position. For this, we reconstructed the relatedness within the groups and subgroups and (3) asked if breeder and helpers’ investment in territory maintenance, brood care, and territory defense relate to kin structure and the gained direct and indirect benefits. The interplay of direct and indirect benefits on helping behavior and fitness prospects is still not well understood in cooperative breeders. Answering these questions will elucidate the relative importance of direct and indirect fitness benefits in the evolution of alloparental care.

Materials and Methods

STUDY POPULATIONS AND SAMPLING

We sampled and genotyped N. savoryi groups in two populations at the Zambian coast of Lake Tanganyika: Kasakalawe (“KK”, 8°46.849' S, 31°04.882' E) and Kasenga (“KS”, 8°42.9' S, 31°08.1' E). The KK study site features a rather homogenous sandy plain with rocks (typically ∅10‐40 cm in size) partly submerged in the sand, at 9.0‐11.5 m depth (Heg et al. 2005a; Josi et al. 2019). In contrast, the KS study site is characterized by layers of stones and large boulders (>1 m diameter), interspersed with patches of gravel and shell debris.

All data were collected by SCUBA diving from February to April and October to November 2003. First, the social status and group membership of each individual was determined based on behavioural observations, home ranges, social interactions, and breeding chamber visits (see Heg et al. 2005a). Individuals are highly territorial, and the home ranges of the breeder females’ subgroups do not overlap with each other (Josi et al. 2021). Within subgroups helpers may have private shelters and their home ranges only overlap to a certain extent, but usually not completely. This allowed us to individually identify these fish using a combination of their social status, their home ranges within the marked territories, their body size, and their individual‐specific, unique color patterns, and markings on the head and flanks (cf. Josi et al. 2020b).

Data at KK were collected from 33 breeding groups containing 495 members (33 breeder males, 60 breeder females, 386 helpers (201 males, 152 females, 33 of unknown sex) and 16 offspring (from 10 different groups)), and 21 group‐independents, summing up to a total of 516 individuals. Data at KS were collected from 10 groups containing 61 members (10 breeder males, 15 breeder females, 36 helpers (18 males, 12 females, 6 of unknown sex)) and 1 group‐independent, summing up to 62 individuals in total. Throughout this paper we refer to “groups” as units encompassing all members within the male breeders’ territory (or “harem”), whereas the term “subgroup” refers to the members of a female breeder's territory. Group members were categorized as breeder males, breeder females, helpers (>15 mm SL), and offspring (<15 mm SL; cf. Groenewoud et al. 2016; Tanaka et al. 2016). “Independents” were fish living singly or in small groups without a breeder female, but they were usually associated with a specific subgroup and occasionally visited by the respective breeder male.

GROWTH AND AGE ESTIMATION

All fish were caught using tent and fence nets, and anaesthetised using clove oil (eugenol, 1 part eugenol dissolved in 4 parts 70% ethanol; Kreiberg 2000). Body size (standard length, SL) of all individuals was determined under water either to the nearest 0.5 mm using a measuring board (KK) or to the nearest 0.1 mm using a calliper (KS). Small tissue samples of all individuals were taken from the dorsal fin and stored in 99% ethanol for subsequent DNA analyses. The fish were sexed by close inspection of the genital papilla. For verification that the sex was correctly identified, a subsample of fish (N = 61) from KS was dissected after completion of data collection. Apart from these fish used for calibration, all other fish were released back to their shelter, where they recovered within a few minutes (cf. Josi et al. 2020b).

At KK, all focal fish were caught and measured after the behavioural observations had been conducted (see below). At KS, all focal fish were caught at least 8 days before the behavioural observations started. Additionally, growth rates were estimated from a subsample of fish in this population (N = 38), which were caught a second time at the end of the second field season, within a period of 31 ± 9 days (mean ± SD; range: 18–51 days). Of these fish (9 breeder males (initial SL: 50.5‐66 mm), 11 breeder females (38‐50 mm), 7 helper males (20‐37 mm), 10 helper females (21‐34 mm), and one offspring (13 mm)) the growth function was estimated. We used the Blumberg growth curve of Skubic et al. (2004) to estimate the age (in days) of each individual in the respective population (cf. Dierkes et al. 2005). Growth rates at KS are comparable to the KK population (Josi et al. 2021).

BEHAVIORAL OBSERVATIONS

Focal observations were performed in the KK population (one 10 min observation per individual: n = 18 breeder males, 21 breeder females, 13 helper males, and 2 helper females) and the KS population (usually three 15 min observations per individual that were averaged and multiplied by 2/3 to allow joint analysis with the KK sample; n = 10 breeder males, 15 breeder females, 18 helper males and 14 helper females). Focal fish were selected haphazardly and covered the range of body sizes of group members (average ± SD: 48.1 ± 9.4, range: 23.7–70.2 mm SL, n = 114). At the beginning of each observation, the observer remained motionless in front of the territory for some minutes to acclimatize the fish to their presence. Due to the previous data collection on colony structure, the fish were habituated to the presence of an observer and behaved normally shortly after their appearance. We recorded the number of aggressive behaviors (overt aggression: biting and ramming; restrained aggression: lateral displays and opercula spreads; Josi et al. 2020a) against conspecific and heterospecific intruders, the number of visits to the breeding chamber (Josi et al. 2020b), and the frequency of territory maintenance behaviours (digging sand or small debris out of the shelter, Josi et al. 2020b). Comparable helping behaviours are shown by all cooperatively breeding Lamprologini (Heg and Bachar 2006; Taborsky 2016; Tanaka et al. 2018a), and are repeatable within individuals over time and context (Chervet et al. 2011; Le Vin et al. 2011; Hamilton and Ligocki 2012; Kasper et al. 2017). We summarized the behavioral frequencies of territory defence against conspecific and heterospecific intruders, breeding shelter visits, and territory maintenance behaviours of each group member to a composite measure of workload (cf. Balshine et al. 2001; Bruintjes et al. 2010). Behavioral frequencies were rounded to the nearest integer value to allow fitting Poisson GLMMs (see “Statistical analysis”).

MICROSATELLITE ANALYSES

Details on DNA extraction, amplification and analyses, as well as references to the 13 microsatellite loci are provided in the Appendix. Ten loci were used for the KK population and 13 loci for the KS population.

STATISTICAL ANALYSIS

All statistical analysis was conducted with SPSS26 and R (R Core Team 2017). Growth rate (mm/day) was correlated to ln[SL] in a GLM to account for exponentially diminishing growth over size (both effects of status ‐ breeder or helper, P = 0.78; and sex ‐ male or female, P = 0.78 were non‐significant). The intercept and slope were used to calculate the average relative age of an individual of a specific body size. This method has been shown to provide reliable age estimates in closely related, similar‐sized lamprologine cichlids (Jungwirth et al. in revision; Skubic et al. 2004). Group size and subgroup size were defined as the number of group members larger than 15 mm SL, including breeders and helpers (Heg et al. 2005a; Josi et al. 2019).

Dyadic estimates of pairwise genetic relatedness (Goodnight and Queller 1999) were calculated using the software Kingroup v2.1 (https://github.com/dmitryako/kingroup; Konovalov et al. 2004). Locus TmoM11 was especially prone to genotyping errors and therefore removed for the relatedness estimations (see Appendix Table 1). In the KS population, we additionally removed the loci NP101, Pzeb4, and UNH002 as they deviated from Hardy‐Weinberg equilibrium. Several microsatellite loci contained ‘rare’ alleles (that is, 19.9% of the alleles occurred in one group only), so a bias‐correction procedure to calculate background allele frequencies for each group separately is not recommended (Konovalov and Heg 2008). Therefore, we calculated for each population the population allele frequencies corrected for overall relatedness in Kingroup version 2.1, and implemented these as an allele frequency block for the analyses (Konovalov and Heg 2008).

Table 1

Number of non‐immigrants (philopatric), possible immigrants, and assured immigrants (‘Immigrant’) per status and sex (% assured immigrants in brackets)

Status, sex	Non‐immigrant	Possible immigrant*	Immigrant
Breeder male	11	21	1	(3.0)
Breeder female	31	8	21	(35.0)
Helper male	161	7	33	(16.4)
Helper female	126	5	21	(13.8)
Helper unknown sex	30	0	3	(9.1)
Independent male	7	3	7	(41.1)
Independent female	2	0	2	(50.0)
Total	368	44	88	(17.6)

*Focal individual has unique genotype inside the group, but philopatry may have been obscured due to rapid turn‐over of breeders and siblings.

The majority of fish caught at KK were genotyped (516 individuals, SL between 6.5 and 65.0 mm; for detailed size distribution see Heg et al. 2005a). The sizes of the 62 fish successfully genotyped at KS (SL between 16.8 and 70.2 mm) were comparable to those genotyped at KK. To correct for potential biotic and abiotic differences between populations that might affect genetic structuring, we tested for population effects on relatedness throughout.

To analyze how relatedness between breeders and helpers varies with helper body size, a General Linear Mixed Model (glmm) was performed. The pairwise genetic relatedness between breeder males and their helpers within the group, as well as the relatedness of breeder females to their helpers in their subgroup, were set as response variables. Helper body size (ln[SL]) was set as a continuous effect and breeder sex as a fixed effect (see also Dierkes et al. 2005). The same analysis was performed on a subset of sexed helpers, that is, by adding helper sex as a fixed effect (excluding n = 120 out of 954 comparisons of helper to breeder relatedness due to uncertain sexing of helpers). Furthermore, we analyzed how the helper relatedness changes with their body size/age. A GLMM was fitted with the pairwise genetic relatedness between helpers of the same subgroup as a response variable and the body size of the larger helper (ln[SL]) as well as the difference in helper size (i.e., ln[SL larger helper minus SL smaller helper], which includes same‐sized helpers), as continuous effects. To test the effect of helper sex, the same analysis was performed on the subset of sexed helpers. The sexes of both helper pairs were included as fixed effects. In both GLMMs ‘group’ was entered as a random effect and “population” as a fixed effect.

Next, we compared the pairwise genetic relatedness of individuals of each social status (breeder males, breeder females, helpers, and independents) among each other with respect to whether individuals were group members of the same subgroup, group members of a different subgroup within the same group, or non‐group members (see Fig. 3). Breeder males were assigned to the subgroup they associated with most of the time. If independents associated with a particular subgroup, they were assigned accordingly. We used the bootstrap data selection and loop macro provided by 26. Bootstrapping was chosen because this approach is less sensitive to the violation of the assumptions regarding the underlying sampling distribution. We calculated four separate anovas on pairwise genetic relatedness of the bootstrap dataset with Dunnett's C post hoc tests (which does not assume equal variances) for each social status (breeder males, breeder females, helpers, and independents). Every possible combination within each social status was tested against each other category (social status × subgroup member/group member from different subgroup/non‐group member; see also Fig. 3). This was repeated 100 times for each social status. The resulting 100 estimates for P‐value and the 95% confidence intervals of the Dunnett's C‐test (lower‐bound and upper‐bound) were averaged. The averaged Dunnett's C‐test results were interpreted as significant if the confidence interval did not include zero.

Figure 3

Average relatedness (± SE) of (A) breeder males (Bm), (B) breeder females (Bf), (C) helpers (He), and (D) independents (In) to subgroup members (black dots), group members of a different subgroup (grey dots), or non‐group members (white dots). Note that for breeder males the subgroup was defined as the one where he spent most of his time, and independents were associated with one subgroup only. Note that within each of the four panels, different letters indicate significantly different categories. The sample size of each bootstrap ANOVA per panel with post hoc Dunnett's C‐test was averaged across 100 repetitions and indicated as mean at the top of the graph.

For the KK population, we furthermore reconstructed parent‐offspring relationships and full‐sib and half‐sib relations contained within each group. The Simpson‐assisted descending ratio algorithm in Kingroup version 2.1 (Konovalov 2006) was used to split the data into groups. Once the dataset was split, we checked for parent‐offspring, full‐sib and half‐sib relationships within each group. That was done in Kingroup using the likelihood‐method of Goodnight and Queller (1999), where the respective relations were compared against the null hypothesis of no relatedness. The utilized likelihood‐method ensured that if an individual was misplaced by the group‐splitting step, it would be automatically filtered out as a non‐relative by the corresponding P‐value of no relatedness. This reconstruction was used threefold: (i) to detect matrilineal or patrilineal inheritance of the group, that is, whether the current breeder had inherited the group; this was possible as usually the other‐sex parent or similar‐aged full or half‐siblings were still present in the group. Nevertheless, this must be regarded as a minimum estimate, because potential mortality of siblings and parents may make inheritance undetectable; (ii) to estimate the minimum and maximum numbers of days each breeder had occupied the breeding position in their respective group, where the minimum is the estimated age of the oldest group member produced (set to zero if no offspring was produced), and the maximum is the estimated age of the next older non‐immigrant (see below) group member not produced by the respective focal breeder. These minimum and maximum estimates were averaged per individual breeder (‘tenure’) and entered into a Kaplan‐Meier life‐table analysis to compare the tenure between breeder males (n = 33) and breeder females (n = 60). The terminal event was set as ‘censored’ when no next older non‐immigrant group member was available (i.e., the tenure must be regarded as a minimum, “open‐ended” estimate, n = 8 breeder males, and 14 breeder females); (iii) to categorize all group members as either non‐immigrants (matrilineal/patrilineal inheritance, or full/half‐siblings present inside the group), possible immigrants (unique genotype, but no older and younger kin‐groups present inside the group to compare against, so rapid breeder turnover may have obscured philopatry), or immigrants (unique genotype plus older and younger kin‐groups present).

Parentage was assigned pairwise using Kingroup to all adult group members (potential producers) versus any group member smaller than these adults (potential offspring), taking the growth rates of both fish into account (i.e., whether the potential producer was actually adult at the time the potential offspring was produced). Adulthood was estimated based on benchmark data derived from the dissection of KS fish (n = 61): (1) males had developed testes from 32.5 mm SL onward, with clearly developed testes from 48.0 mm SL onwards; (2) females had active ovaries from 31.5 mm SL onward, with the ovaries containing large, developed oocytes from 43.9 mm SL onward. Parentage was assigned based on zero or one mismatch. The latter occurred only among nine breeder males (out of 161) and four breeder females (out of 89) with one of their offspring each. The number of offspring produced by male and female helpers within subgroups was set in relation to their relatedness to the breeders (same or opposite sex) with a Poisson GLMM. Only helpers that were sexually mature (i.e., females larger than 38.5 mm and males larger than 43.5 mm) were included in these analyses (see table 2) and population was set as a random effect. The effect of within‐subgroup kinship on the workload (n = 114, n = 50, n = 64) was tested for breeders and helpers separately using Poisson GLMMs with log‐link and “group” as a random effect. The workload was set as a response variable in all models. For the breeders, we included parentage (number of produced subgroup members), subgroup size, and mean relatedness to all helpers not produced per subgroup, as well as sex (male or female) as a predictor. For helpers, subgroup size, sex, and mean between‐helper relatedness of the respective subgroup (which tests for kinship effects) were included as a predictor. Variables were removed from the model if the difference between the Akaike Information Criterion (AIC) values was equal to or larger than 2.

Table 2

Parentage in N. savoryi: total number of group members produced inside the focal individual's main subgroup or other subgroups within the same territory, in relation to the total number of group members

			Group members produced in:		Number of group members b
Role and sex	n Producers	Body size producers mean (mm) ± SD (range)	Subgroup	Group a	Subgroup	Group a
Breeder male	31	59.1 ± 5.2 (50.0‐67.9)	109 (32.1%)	18 (11.6%)	340	155
Breeder female	37	45.5 ± 3.9 (40.0‐56.8)	66 (25.1%)	17 (8.5%)	263	199
Helper male	12	52.4 ± 4.8 (43.5‐62.2)	24 (6.1%)	9 (3.1%)	392	292
Helper female	5	41.3 ± 3.5 (38.5‐47.2)	6 (7.9%)	0 (0.0%)	76	37
Independent male	1	52.0	1 (5.0%)	0 (0.0%)	20	27
Independent female	0	−	‐	0 (0.0%)	0	3
Total	86	51.1 ± 7.9 (38.5‐67.9)	206	44

a In other subgroup(s) than the main subgroup in which the focal individual spent most of its time.

b Sample sizes of subgroup or group members (from other subgroup(s)) they could have potentially produced based on age difference of at least 210 days; that is, the focal fish and thus potential parent was approximately 31.0 mm SL or larger at the estimated time point of spawning.

Results

Growth rate per day significantly declined with body size of the fish (Fig. 1A, β ± SE: intercept 0.38 ± 0.129, P = 0.003; slope ‐0.087 ± 0.033, P = 0.011). These growth data were used to estimate the relative age of each fish based on its body size. Note that we do not have growth estimates for very large helpers (>40 mm SL), and the two outliers with very high growth rates only marginally influenced the coefficients, which remained significant after the outliers were excluded from the analysis (intercept 0.34, P < 0.001; slope −0.081, P < 0.001).

Figure 1

(A) Growth rate per day (mm) of offspring (triangle), helper females (white circles), helper males (white squares), breeder females (black circles), and breeder males (black squares), depending on their initial body size (abscissa). Regression line shows the average body size vs growth rate relationship. Negative values originate from measurement errors within and between observers. (B) Relatedness of the helpers to either the breeder male (black dots, bold line) or breeder female from the same subgroup (white dots, thin line). Depicted are means ± SE per 1 cm size class with the two regression lines from a glmm and estimated age in days (derived from A). Sample sizes are indicated inside the graph.

HELPER TO BREEDER RELATEDNESS

Helper to breeder relatedness significantly declined with body size and thus the age of the helpers (Fig. 1B, glmm, data from 41 groups, effect of ln[SL]: F 1,945.2 = 26.5, P < 0.0001; estimate ± SE: −0.137 ± 0.027). Helpers were more closely related to the breeder male than to the breeder female (Fig. 1B, glmm: F 1,934.8 = 8.04, P = 0.005; estimate ± SE to the males: 0.050 ± 0.018, with relatedness to the females set to zero; intercept ± SE: 0.618 ± 0.108, F 1,890.6 = 42.43, P < 0.0001). This effect was similar in both populations (P = 0.72). These results did not change when helper sex was included in the model, which by itself showed no significant effect (P = 0.34).

HELPER TO HELPER RELATEDNESS

Helper to helper relatedness decreased with both absolute helper size (ln[SL]: glmm, (F 1,3803.2 = 80.0, P < 0.0001; estimate ± SE: −0.172 ± 0.019) and difference in helper size (ln[difference in SL]: GLMM, (F 1,3801.4 = 42.6, P < 0.0001; estimate ± SE: 0.086 ± 0.013), (see Fig. 2A and B). There was no difference between populations (GLMM, F 1,36.3 = 1.5, P = 0.22; data of 39 groups; intercept: F 1,922.2 = 123.8, P < 0.0001; estimate ± SE: 0.865 ± 0.09). When helper sex was included in the model, it did not reveal significant effects (both P = 0.46).

Figure 2

(A) Average relatedness (± SE) of the helpers per 1 cm size class to other subgroup member helpers dependent on body size; also indicated is their estimated age in days. (B) Model predicted estimates were used to depict helper‐to‐helper relatedness within subgroups as a landscape.

RELATEDNESS TO GROUP MEMBERS AND SUBGROUP MEMBERS

Individuals were more related to members of their own subgroup than to group members belonging to different subgroups, whereas relatedness to non‐group members was virtually zero (see Fig. 3 for average sample sizes per bootstrap anova, 100 bootstrap anovas per status, average P‐values: breeder males, df = 9, P < 0.001; breeder females, df = 8, P < 0.001; helpers, df = 5, P < 0.001; independents, df = 1, P < 0.001). Post hoc bootstrap Dunnett's C‐tests comparing the different categories pairwise were only significant for those combinations where we had reasonable sample sizes (see Fig. 3). These involved mainly combinations with the helpers: (i) Breeder male versus helper relatedness did not differ between group members belonging to the same or different subgroups, and both were significantly higher than relatedness to non‐group members (Fig. 3A). (ii) Breeder female versus helper relatedness declined significantly from group members of the same subgroup to group members from other subgroups, and to non‐group members (Fig. 3B). (iii) Helper vs helper relatedness declined significantly from group members of the same subgroup to group members from other subgroups, and to non‐group members (Fig. 3C). (iv) Helper versus independent relatedness was significantly higher for subgroup members than for both group members from other subgroups and non‐group members (Fig. 3C). Helper versus helper relatedness did not differ from helper vs independent relatedness from the same subgroup (Fig. 3C). (v) Finally, independents were highly related to independents from the same subgroup, whereas relatedness to independents from other groups was zero (Fig. 3D; this result was significant despite the small sample sizes).

BREEDER TENURE AND TERRITORY INHERITANCE

The tenure of breeder males (median estimate 257.5 days, n = 33; range: 36.5–737 days) was significantly longer than that of breeder females (191.5 days, n = 60; range: 11–331 days, Kaplan‐Meier Wilcoxon‐Gehan statistic = 5.3, P = 0.021). There was a non‐significant trend that females inherited the breeding territory more often than males (matrilineal inheritance (MLI): 31 out of 60; patrilineal inheritance (PLI): 11 out of 33; χ2‐test, χ2 = 2.89, df = 1, P = 0.089). PLI did neither change with group size nor with male body size (Logistic Regression (GLM): intercept: 9.06 ± 6.09; effect of group size, df = 1, β ± SE: −0.101 ± 0.094, P = 0.28; effect of body size, df = 1, β ± SE: −0.156 ± 0.121, P = 0.2). In contrast, MLI was more likely to occur in large groups (Fig. 4, Logistic Regression (GLM): intercept: −2.369 ± 5.446; effect of group size, df = 1, P = 0.04, β ± SE: 0.08 ± 0.039; effect of female body size, df = 1, P = 0.85, β ± SE: 0.025 ± 0.036).

Figure 4

Matrilineal inheritance of groups depending on group size. Multiple cases are indicated by increasing symbol sizes (1 to 10 cases). Also depicted is the model‐predicted logistic regression equation.

IMMIGRATION

Overall, 16.5% of the group members were immigrants from different groups, with a further 8.6% of individuals being potential immigrants, but their unique genotypes could also have been due to the rapid turnover of breeders and siblings from their own group (Table 1). To investigate the relationship between immigration and status (breeder, helper, or independent) as well as body size, a Multinomial Regression model (MR) was fitted including data of all individuals with known sex (n = 467). As a response variable, the immigration state (non‐immigrant, possible immigrant, or immigrant) was included and status as well as body size were set as predictors. The interaction between status and body size was significant (P < 0.001), indicating that the relationship between the likelihood of being an immigrant and body size varied depending on the status. Based on this significant interaction we subsequently calculated MRs for each status separately, with the likelihood of being an immigrant as a response variable and body size as a predictor. Larger breeders were less likely to be immigrants (P = 0.024, Fig. 5A), whereas larger helpers were more likely to be immigrants than smaller ones (P < 0.001, Fig. 5B). There was no relationship between the likelihood of being an immigrant and body size for independents (P > 0.1, Fig. 5C). The likelihood of being a “possible immigrant” was higher for both larger breeders and helpers (P < 0.001), supporting the interpretation that a failure to resolve the origin of unique genotypes in older individuals may be due to group member turnover, and that some true immigrants may be hidden among these individuals.

Figure 5

Proportion of non‐immigrants (black), possible immigrants (grey), and immigrants (white) per size class for (A) breeders, (B) helpers, and (C) independents. Sample sizes are indicated inside the bars.

PARENTAGE

We determined either the father (n = 121 individuals), the mother (n = 49), or both parents (n = 40) for individual group members based on one or zero mismatches in the pairwise parent‐offspring assignments. Parentage was skewed toward the dominant pair, but helpers of both sexes had some reproductive share (Table 2). The dominant male and the dominant female produced a significantly higher proportion of the group members inside their “own” subgroup compared to other subgroups (Table 2, χ2 = 23.3 and 21.0; both df = 1, both P < 0.001, n = 109(own) vs. 18(other), n = 66(own) vs. 17(other)). Helper males, but not helper females, also produced offspring in other subgroups (Table 2). The proportion of offspring produced within the own subgroup compared to other subgroups of the same harem did not significantly differ (χ2 = 3.4 and 3.1; both df = 1, P = 0.066 and 0.079, n = 24(own) to 9(other) and n = 6(own) to 0(other)). Helper males had a higher reproductive share with increasing relatedness to the breeder female (GLMM: n = 48; χ2 = 7.51; df = 1, P = 0.006). Only few helper females had a reproductive share, which was higher if they were unrelated to the breeder male (GLMM: n = 24; χ2 = 5.24; df = 1, P = 0.02). The reproductive share of female helpers correlated negatively with their relatedness to the same‐sex dominant (GLMM: n = 24; χ2 = 3.8; df = 1, P = 0.05), whereas there was no relation between the reproductive share of male helpers and their relatedness to breeder males (GLMM: n = 48; χ2 = 0.08; df = 1, P = 0.78).

KINSHIP AND WORKLOAD

Overall, the mean kinship within subgroups was significantly higher than zero (n = 81 subgroups, r = 0.135 ± 0.125 SD, pairwise relatedness values, range: −0.12 ‐ 0.52, one‐sample t‐test t = 9.72, df = 80, P < 0.001), but with high variation of relatedness between subgroup members (mean within‐group standard deviation 0.244 ± 0.079 SD, range: 0.01 ‐ 0.57, n = 75, excluding six subgroups with only two individuals genotyped). The workload was determined for 64 breeders and 50 helpers. It varied between 1 and 36 behavioral events (rounded to the nearest integer, average ± SD: 8.46 ± 6.05). Breeders worked harder with increasing numbers of own (genetic) offspring present in their subgroup (Fig. 6A, Table 3, β ± SE: 0.063 ± 0.023, P = 0.003). Breeder sex (P = 0.25), relatedness to helpers not directly produced (P = 0.68) and subgroup size (P = 0.99) did not affect breeder workload (for details, see Table 3). The helpers’ workload decreased with increasing subgroup size, that is, when there were more individuals sharing duties (Fig. 6B, Table 3, β ± SE: −0.092 ± 0.029, P = 0.002). Additionally, with increasing relatedness among helpers of a subgroup, focal helpers worked less (Fig. 6C, Table 3, β ± SE: −1.763 ± 0.795, P = 0.02). Overall, male helpers took on a higher workload than female helpers (Table 3).

Figure 6

(A) Parentage effects (number of produced subgroup members) on the per capita workload of breeders. Per capita workload of helpers is represented in relation to (B) the size of the subgroups and (C) the relatedness among helpers. The regression lines are based on the final model in Table 3 and predicted 95% confidence intervals are represented in grey.

Table 3

Workload and kinship; results of the final model (Poisson GLMM, dependent variable: workload; n = 114 focal individuals). Females are the reference for the factor of sex. Significant P‐values are marked in bold

Fixed effects	Estimate	±	SE	df	z‐value	P‐value
Final model breeders:
Intercept	2.056	±	0.174		11.832
Parentage	0.063	±	0.022	1	2.888	0.003
Sex	0.112	±	0.097	1	1.162	0.25
Relatedness to subgroup helpers	−0.134	±	0.329	1	−0.407	0.68
Dropped variable from full model:
Subgroup size	0.00006	±	0.018	1	0.003	0.99
Final model helpers:
Intercept	1.66	±	0.199		8.343
Subgroup size	−0.092	±	0.029	1	−3.131	0.002
Sex	0.43	±	0.152	1	2.837	0.004
Relatedness to subgroup helpers	−1.763	±	0.795	1	−2.218	0.022

Discussion

The social and genetic group structure of a population can provide important information about evolutionary mechanisms underlying complex social organization (Cockburn 1998; Riehl 2013; Tanaka et al. 2018b). Here, we show that in a cooperatively breeding cichlid, breeder tenure, territory inheritance, immigration, and reproductive skew are linked to the level of within‐group relatedness, which in turn bears upon the amount of alloparental care provided by subordinate helpers. Group composition in social cichlids is a dynamic process causing relatedness patterns within groups to change with time; the older helpers get, the lower on average is their relatedness to other group members (Dierkes et al. 2005). Due to this dynamic, helper age determines the relative importance of direct and indirect fitness benefits of group living and cooperation.

N. savoryi lives in polygynous groups or “harems” divided into subgroups containing female breeders and subordinate helpers. Relatedness differs between group members, declining with the increasing age of helpers (Fig. 1B). Overall, breeder females and helpers were more related to each other in their own subgroup compared to members of other subgroups within the same harem, let alone to members of different groups (Fig. 3). Genetic relatedness of helpers to their breeder male did not differ between the subgroups, but it was higher to the male breeder of their own group than to breeder males of other groups (Fig. 3). Remarkably, helpers were significantly more closely related to their breeder male than to the breeder female of their subgroup (Fig. 1B). This is consistent with the finding that breeder males had a significantly lower turn‐over rate than breeder females (due to death, emigration, or usurpation by unrelated individuals). Helper to helper relatedness declined with their age, and it was significantly higher within than between subgroups (Figs. 2 and 3C).

Dispersal patterns of individuals living in complex social systems are often obscure and difficult to unveil. Recent field evidence suggested that N. savoryi helpers may disperse to a neighboring female subgroup within the same male harem (Josi et al. 2021), and a breeder female may establish a subgroup within the group of a closely related male (either the father or the brother; Josi et al. 2019). In line with these findings, the present results indicate that only 16.5% of the group members were immigrants from outside the male's group (or 25.1% if we count all individuals with unique genotypes; see Fig. 5 and Table 1). Immigration was more difficult to detect for larger, that is, older group members (particularly breeders), because the likelihood of finding unique genotypes inside a group increases progressively over time when all close, non‐offspring relatives have disappeared from the group. The low immigration rate implies that both males and females have a good chance to inherit a breeder position within a group. Nevertheless, if we look at the acquisition of breeder status within groups, about 48% of breeder females’ and 67% of breeder males’ positions were acquired by non‐group members, which might indicate that neighboring individuals of both sexes compete with each other, for example, to take over territories of higher quality.

Similar to other cooperative breeders (e.g., Dierkes et al. 2005; Leadbeater et al. 2011) female inheritance was more likely to occur in large groups (see Fig. 4). Such a biased probability of territory inheritance may result from a greater number of potential candidates or a higher survival rate of group members in large groups. It might further indicate that large groups are particularly valuable for breeders as they render higher reproductive rates, enhanced survival for group members, and longer group persistence (Heg et al. 2005b; Jungwirth and Taborsky 2015).

We should like to point out that estimating immigration rates, classification of immigrants, tenure times, and territory inheritance from relatedness patterns may be compromised by reproductive parasitism from outside the group and emigration or death from previous producers of offspring. However, as the relatedness of helpers to the breeder male was higher than to the breeder female within their subgroup, it seems likely that extra‐pair parentage as a potential source of error can be regarded as rather low.

In N. savoryi, group membership is essential, as individuals of all sizes face a high risk of predation (Heg et al. 2008; Josi et al. 2020a). Consequently, helpers gain direct fitness benefits through increased protection from predators due to safety in numbers, vigilance of other group members and shared defence, as suggested also in other cooperative breeders (Taborsky 1984; Heg et al. 2004; Caro 2005; Santema and Clutton‐Brock 2013; Groenewoud et al. 2016; Teunissen et al. 2020). Additionally, especially older helpers, can gain further immediate direct fitness benefits by sharing in reproduction (see Table 2), and future benefits by queuing for the breeding position (Griffin and West 2002; Leadbeater et al. 2011; Taborsky 2016; Kingma 2017; for matrilineal inheritance see Fig. 4). Indeed, large helper males sired group members in 20.5% of the cases where the father was successfully detected (n = 161), while helper females had produced young in 6.7% of the cases where the mother of the offspring could be successfully detected (n = 89). Such reproductive share may also explain why N. savoryi helpers assist their breeders in offspring care (Josi et al. 2020a,b). Alternatively, helping can serve as a rent paid to dominants to be allowed to remain in the group (“pay‐to‐stay”; Kokko et al. 2002; Zöttl et al. 2013; Fischer et al. 2014; Naef and Taborsky 2020a, 2020b; Trapote et al. 2021). These results add to the increasing evidence that direct fitness benefits can foster the evolution of cooperative breeding and complex sociality (Clutton‐Brock 2002, 2009; Griffin and West 2002; Riehl 2013; Field and Leadbeater 2016; Komdeur et al. 2016; Taborsky et al. 2016; Kingma 2017; Downing et al. 2018).

The mixed kin‐structure of N. savoryi groups and subgroups resembles that of many cooperatively breeding birds and mammals in which complex relatedness and parentage structures have been observed (e.g., Painter et al. 2000; Riehl 2013; Shen et al. 2016; Vitikainen et al. 2017; Leedale et al. 2018; Lukas and Clutton‐Brock 2018). Within cooperatively breeding lamprologine cichlids, the kin structure of N. savoryi groups lies between species characterized by high relatedness among group members (e.g., N. obscurus; Tanaka et al. 2015) and those featuring groups of largely unrelated individuals (e.g., Julidochromis ornatus; Awata et al. 2005). This raises the question which direct and indirect benefits individuals gain from group membership and helping. At one end of the spectrum, helpers mainly gain direct fitness benefits by paying for valuable group membership, participation in reproduction and territory inheritance (J. ornatus). At the other end, they may mainly benefit from indirect fitness gains by helping to raise close kin (N. obscurus). Hitherto, existing data suffice to estimate the relative importance of direct and indirect fitness effects of delayed dispersal and alloparental care in only few cooperatively breeding vertebrates (Kingma et al. 2011; Jungwirth and Taborsky 2015; Green and Hatchwell 2018). Therefore, N. savoryi is a good model to study potential selection mechanisms underlying complex social organization and cooperation. Age seems to be the crucial parameter determining the relative importance of direct and indirect fitness benefits to helpers in N. savoryi. Relatedness between helpers and breeders declines with increasing helper age due to breeder mortality and territory takeover of immigrants. Therefore, kin selected benefits are skewed towards young helpers, if they assist their parents in raising their full or half siblings (see Fig. 1B; for a similar situation in N. pulcher, see Bruintjes and Taborsky 2011). Young helpers face a particularly high predation risk due to their small size (Heg et al. 2004). Especially these small individuals benefit from remaining with close kin, as breeders invested more in territory defense and brood care when the young in their territory were closely related (see Fig. 6A). Such increased investment in groups containing highly related subordinates further results in additional benefits for helpers, as they are able to reduce their individual workload (Fig. 6C). Finally, relatedness to the dominants might further affect the parentage of the helpers, either through negative effects of inbreeding or through mutual kin selected benefits of co‐breeding. The risk of inbreeding and consequential inbreeding avoidance when the opposite sex dominant is a relative may select for increased helper dispersal (Beck et al. 2008; Nelson‐Flower et al. 2012). In contrast, an increase in reproductive opportunities for helpers can select for delayed dispersal, cooperative brood care, and participation in reproduction (Hellmann et al. 2016). Indeed, we found that the reproductive share of helper females declined with increasing relatedness to the breeder male. As helper females always reproduced with males (breeders and subordinates) of the same group, this effect might be explained by inbreeding avoidance. In contrast, male helpers that were related to the dominant female had a higher reproductive share than unrelated males. It is important to note that male helpers did not reproduce with the breeder female in their own group. This suggests that breeder females do not prevent related male helpers from reproducing, which yields indirect fitness benefits to them. From our data, it was not possible to reconstruct the genetic mother of the male helpers’ offspring, which may indicate that the breeder females evicted those mating partners.

A functional relationship between intragroup relatedness and reproductive skew among sexually mature group members has been predicted by transactional skew theory (Reeve and Keller 1997; Johnstone 2000). In particular, concession models predict a positive association between relatedness and reproductive skew. This is based on the logic that if helpers can gain only little indirect fitness benefits by alloparental investment due to low relatedness between donors and beneficiaries, incentives to stay and participate in brood care need to be provided to raise their direct fitness, such as taking a share in reproduction (Vehrencamp 1983; Keller and Reeve 1994). This relationship between relatedness and reproductive skew is compatible with patterns observed in cooperatively breeding lamprologini. For example, in J. ornatus and J. transcriptus relatedness within groups is low, and male and female helpers gain substantial parentage (41‐56% in J. ornatus: Awata et al. 2005; Heg and Bachar 2006; 20–100% in J. transcriptus: Kohda et al. 2009). In contrast, in N. pulcher, where relatedness within groups is mixed, the share in reproduction of subordinates is much lower (2.5‐19% Dierkes et al. 1999; Heg et al. 2006; Heg and Hamilton 2008; Hellmann et al. 2015). Finally, in N. obscurus, where relatedness levels are particularly high within groups, subordinate group members have underdeveloped gonads suggesting reproductive austerity (Tanaka et al. 2015). Similar to N. pulcher, N. savoryi is apparently positioned in between the extreme cases, with 6–8% of young produced by helpers within their own subgroup. Nevertheless, this rough concordance between the predictions of concession models and the described patterns does not necessarily imply that the underlying hypothesis explains reproductive participation in all lamprologines alike. This note of caution can be illustrated by the existence of outliers from the general trend; in N. multifasciatus, for example, preliminary data suggest that high intragroup relatedness goes along with a rather high reproductive share of helpers (Kohler 1998). Hence, it seems that conventional skew models cannot fully account for the complexity of evolutionary mechanisms involved in reproductive skew among members of fish groups (Taborsky 2009).

There is surprisingly little evidence for a positive relationship between intragroup relatedness and reproductive skew in other taxa. In hymenoptera, there are some examples showing the predicted positive relationship between intragroup relatedness and reproductive skew (Reeve et al. 2000; Reeve and Keller 2001; Lucas et al. 2011; Andrade et al. 2016), and the same holds for some birds (Jamieson 1997; Koenig et al. 2009). However, a comparative analysis of all 83 species of cooperatively breeding birds from which genetic data existed suggests that inbreeding avoidance rather than relatedness‐driven reproductive concessions can explain the relationship between kinship and helper reproduction (Riehl 2017). In cooperatively breeding mammals, reproductive skew may be comparatively low because of viviparity hampering full control of reproduction by dominants (Raihani and Clutton‐Brock 2010), and the relationship between intragroup relatedness and reproductive skew may again be driven mainly by inbreeding avoidance (Cooney and Bennett 2000; Nichols 2017; Wikberg et al. 2017). In lamprologine cichlids, dominant group members are largely in control of reproduction due to inherent size differences, and relatedness patterns vary between closely related species sharing the same ecology. Hence, they seem to be a highly suited target for futures studies on the significance of reproductive concessions in the evolution of cooperative breeding and complex social organization. In addition, N. savoryi represents a suitable model for studying the ability of brood care helpers to fine‐tune their support to the demands of breeders and to their own costs, as well as to the benefits breeders gain from the presence of helpers in dependence of the relatedness between them.

AUTHOR CONTRIBUTIONS

D.H., T.T., M.K., and M.T. conceived the study. D.J., D.H., D.B., and D.A.K. analyzed the data; D.H., M.T., D.J., and J.G.F. wrote the manuscript, which was revised and approved by all authors.

DATA ARCHIVING

Data are archived in the Dryad Digital Repository: https://doi.org/10.5061/dryad.t1g1jwt30.

ETHICAL NOTE

Data collection followed the ASAB/ABS guidelines for the treatment of animals in behavioural research and teaching, and the regulations of the “Zambian prevention of cruelty to animals” act.

CONFLICT OF INTEREST

The authors declare no conflict of interest.

Associate Editor: M. E. Maan

Handling Editor: A. G. McAdam

Table A.1

Number of adults successfully typed (n), number of different alleles (Alleles), observed (H O) and expected homozygosities (H E), tests for Hardy‐Weinberg equilibrium (HWE), estimated null allele frequencies (Null), respectively, per locus for the two study populations. Based on the Arlequin 3.01 analysis of adult breeders

Locus	n	Alleles	H O	H E	HWE	Null	Reference*
Kasakalawe population (KK)
NP007	64	8	0.484	0.52	0.44	0.041	1
NP101	64	16	0.625	0.645	0.36	0.018	2
NP773	64	12	0.734	0.686	0.7	–0.049	1
Pzeb3	64	12	0.516	0.571	0.03	0.041	3
Pzeb4	64	6	0.453	0.514	0.51	0.051	3
TmoM11	64	19	0.563	0.927	<0.001	0.237	4
TmoM25	64	19	0.953	0.911	0.44	–0.027	4
ULI2	64	16	0.563	0.522	0.8	–0.060	1
UME003	64	24	0.875	0.935	0.38	0.028	5
UNH154	64	4	0.109	0.121	1	–0.019	6

1 = Schliewen et al. 2001; 2 = Brandtmann et al. 1999; 3 = Van Oppen et al. 1997; 4 = Zardoya et al. 1996; 5 = Parker et al. 1996; GenBank accession numbers for 6 = G12306 and 7 = G12259.