### Study species

*V. moorii* is a small herbivorous cichlid endemic to Lake Tanganyika in East Africa. Male-female pairs defend shallow rocky territories ranging in size from 1 to 4 m^{2}, where they deposit eggs and raise broods of fry [26,27,28]. Parents raise one brood at a time for approximately 100 days [57], during which time they defend their young against territorial intrusions from con- and hetero-specifics [27]. Although the *V. moorii* mating system is characterized by extreme paternity loss, paired males engage in very little cuckoldry themselves as the vast majority of cuckoldry is perpetrated by unpaired males in the population [31].

### Field collections and microsatellite genotyping

We identified *V. moorii* territories within a study quadrat (~ 100 m × ~ 50 m; depth range, 1.7–12.1 m) in the south of Lake Tanganyika, Zambia (8° 42′ 29.4″ S, 31° 07′ 18.0″ E) over the course of two sampling excursions, from September 22 to October 28, 2015 (corresponding with the dry season), and from April 4 to 20, 2016 (corresponding with the rainy season). The genetic samples collected from these field seasons have been previously analysed by Bose et al. [31] to investigate patterns of extra-pair paternity between paired and unpaired males. In the current study, however, we use these data to answer a very different set of research questions. See ref. [31] for more detailed descriptions of sampling methods and genetic analyses.

In brief, we identified parental care-giving pairs of *V. moorii* in a study quadrat in the field while on scuba. The parents were captured using gill nets (though several individuals evaded capture), measured for total length (to the nearest 0.1 cm), fin clipped, and released back to their territories. Furthermore, for each adult pair, we measured the distances from their territory to their nearest 2–3 neighbors’ territories (center to center) and used this to generate a map showing the locations of, and relative distances between, territories. During the dry season, we sampled paired males from 74 territories and paired females from 85 territories. During the rainy season, we sampled paired males from 45 territories and paired females from 49 territories. We also captured fry from each care-giving pair and sacrificed them in MS-222 (1 g/L lake water). For this study, we included 821 fry from 32 pairs in the dry season, and 1129 fry from 38 pairs in the rainy season. Fry and fin clips were stored in 99.9% ethanol and transported back to the lab for parentage analyses using microsatellite genotyping. This work was carried out with the permission of the Fisheries Department of Zambia and under a study permit issued by the government of Zambia.

We extracted DNA from fry tissues using a standard Chelex protocol [58] and from fin clip tissue using an ammonium acetate precipitation protocol [59]. All adult fin clips and most fry were genotyped at 14 microsatellite loci though 4 broods (138 fry in total) were genotyped at 9 microsatellite loci only. The microsatellite loci used were Pmv17 [60], TmoM11 [61], Pzeb3 [62], UNH2075 [63], Ppun21 [64], Ppun9 [64], Hchi59 [65], Hchi94 [65], UME002 [66], Pmv13 [60], UME003 [66], UNH908 [67], Ppun5 [64], and Ppun20 [64]. These markers were highly polymorphic, with an average of 21.2 alleles per locus and a mean expected heterozygosity of 0.88. Exclusion probability (assuming the mother was known) calculated across loci amounted to 0.9999997 (calculated in GERUD2 [68]). All loci complied with Hardy–Weinberg expectations after correcting for multiple testing (see ref. [31] for additional details).

### Parentage analyses and construction of cuckolder genotypes

Parentage analyses were conducted with the aid of COLONY (v 2.0.6.1, ref. [69]), using population allele frequencies obtained from the adults captured in each season (*N* = 219 in the dry season, *N* = 98 in the rainy season). COLONY used the multi-locus genotypes of the care-giving adults as well as the sampled fry to identify fry as either within-pair or extra-pair, and then split the extra-pair fry into full-sib groups that were each sired by a different extra-pair male (see ref. [31] for details). With the maternal genotypes known, we could also construct the multi-locus genotypes of cuckolder males based on the extra-pair fry that they had sired. Generally, for each group of extra-pair full-sibs, the non-maternal alleles at each locus were assigned to the cuckolder father. However, when cuckolders sired too few offspring, their genotypes could not be fully reconstructed, and we focused only on cuckolders for which we could unambiguously identify both alleles for at least 10 of the 14 loci (although note that our relatedness estimates between parents and offspring take all cuckolders into account: see below). When only one paternal allele at a locus could be identified in a given a full-sib group, it was unclear whether the cuckolder father was homozygous at that locus or whether the fry inherited only one of his alleles. In this case, our assignment of alleles depended on how many extra-pair fry the cuckolder had sired; if the male had sired fewer than eight fry, we scored his genotype as being potentially heterozygous (i.e. ‘paternal allele observed in offspring’/‘unknown’), but if he had sired at least eight fry, we scored his genotype as being homozygous for the observed paternal allele. Furthermore, when the mother and an extra-pair fry shared the same heterozygous genotype at a locus, for example X/Y, it was not possible to identify which allele was inherited from the mother or the cuckolder father. In this case, the paternal allele had two possibilities, either X or Y, and examination of additional extra-pair fry genotypes was required to resolve this ambiguity. However, if the ambiguity persisted after assessing all of the cuckolder’s fry, we assigned the paternal allele as the one that had a higher frequency within the population (unless the difference in population allele frequencies was < 1%, in which case the allele was recorded as ‘unknown’). Thus, we attained multi-locus genotypes directly for paired individuals and indirectly for cuckolders, provided that the cuckolders fertilized a sufficient number of offspring. We then used the ‘Demerelate’ r package (v. 0.9 – 3, ref. [70]) to calculate two symmetrical estimates of relatedness between all individuals in our dataset: *r*_{QG} (following ref. [71]) and *r*_{LR} (following ref. [72]). We used both estimators, *r*_{QG} and *r*_{LR}, because simulations based on the allele frequencies in our dataset showed that *r*_{LR} outperformed other indices for unrelated pairs, whereas *r*_{QG} did better for closely related pairs (as shown previously in refs [73, 74].).

### Do *V. moorii* exhibit more reproduction among relatives than expected under random mating?

To begin, we tested whether relatedness estimates differed between the two seasons (i.e. sampling excursions). We used two linear models to test whether relatedness within pairs (estimated by *r*_{QG} and *r*_{LR}) differed between the dry and rainy seasons. We similarly used two linear mixed effects models (LMMs, ‘lme4’ r package, v. 1.1 – 16, ref. [75]) to test whether relatedness between paired females and their extra-pair (cuckolder) males differed between the seasons. We included ‘Female ID’ as a random intercept term to account for some females mating with multiple cuckolder males. Lastly, we also used two LMMs to test whether the relatedness between paired males and their cuckolders differed between the seasons. We included ‘Paired male ID’ as a random intercept to account for some males being cuckolded by multiple extra-pair males. We applied Yeo-Johnson power transformations where necessary to improve normality of the models’ residuals. Overall, we found no significant differences between the seasons (all *p* ≥ 0.37). We therefore pooled our samples from both seasons for the remainder of our analyses unless otherwise specified.

We used a series of randomization tests to test for elevated pairwise relatedness between all three parties in the *V. moorii* mating system: paired males, paired females, and cuckolder males. In particular, we were interested in whether the (i) mean, (ii) variance, and (iii) skewness of our observed distributions of pairwise relatedness estimates were greater than would be expected under random mating. For example, a distribution of observed relatedness estimates with an elevated mean would suggest that *V. moorii* pair with (or cuckold) relatives more often than expected by chance. On the other hand, an observed distribution with elevated variance and positive skewness would specifically indicate that there were more *high-relatedness* pairings in the observed data than would be expected by chance. In combination, all three measures (mean, variance, and skewness) give a comprehensive understanding of the pairwise relatedness structure between individuals in a population.

To test for elevated relatedness within social pairs, we first calculated *r*_{QG} and *r*_{LR} between the paired males and females that we observed in the field during the two seasons. We then recorded the mean, variance, and skewness of the observed distribution of relatedness estimates (‘moments’ r package, v. 0.14, ref. [76]). Next, we simulated random pairings between the males and females. We randomly assigned our observed males and females to new partners (i.e. re-sampling without replacement) creating the same number of social pairs as in our observed dataset though we ensured that random pairs were only formed between individuals sampled during the same season. Again, we calculated *r*_{QG} and *r*_{LR} from each new pairing and recorded the mean, variance, and skewness from the resulting distribution. We repeated this randomization process 10,000 times to create null distributions against which to compare our observations. We computed *p* values as the proportion of the randomized trials that yielded mean, variance, and skewness values *more* extreme than our observed data.

We followed a similar process to test for elevated relatedness between paired females and cuckolder males. We first calculated *r*_{QG} and *r*_{LR} between our observed paired females and each one of their extra-pair males (as some females spawned with multiple cuckolders). Note again that the cuckolders used here had each sired a sufficient number of offspring for us to reconstruct their genotypes. We then took the average of the relatedness estimates between each female and her *set* of extra-pair males (i.e. cuckolder males). From this data, we recorded mean, variance, and skewness. We then randomly assigned a new set of extra-pair males to each paired female in our observed data. Again, we ensured that individuals could only be assigned together if they were sampled during the same season. Furthermore, it is exceedingly rare for individual cuckolder males to successfully cuckold in more than one territory over the time span of a single brood cycle (approximately 100 days, ref. [31]), and so for each permutation of our data here, we assumed that each cuckolder male would sire offspring with only a single female. We also ensured that each reshuffling of the data preserved the distribution of extra-pair males per female that we originally observed in each season. We took the average of the relatedness estimates between each female and her new set of randomly assigned extra-pair males, and then, we recorded mean, variance, and skewness. This randomization process was repeated 10,000 times, and *p* values were computed as the proportion of randomized trials yielding mean, variance, and skewness values *more* extreme than our observed data. It is important to note that for this randomization test, we omitted broods where the paired male had 0% paternity (6 out of 70 broods). This was done because we could not be sure whether these incidents indicated 100% cuckoldry or a territory takeover. Omitting these broods ensured that, in the event of a takeover, we would not mistake the previous paired male for a cuckolder. We followed a similar procedure to test for elevated relatedness between paired males and their set of cuckolder males. We similarly omitted broods of 0% paternity in order to guarantee that our analyses compared relatedness between males and cuckolders that were both present at the time of spawning. Due to multiple comparisons, we used the Benjamini–Hochberg procedure for controlling false discovery rates following ref. [77], setting the false discovery rate to 10% [78].

### Does parent-offspring relatedness deviate from expectations?

We calculated *r*_{QG} and *r*_{LR} between paired adults and each of the within-pair and extra-pair fry under their care. If paired males were related to their female partners, then we would expect their relatedness estimates to within-pair offspring to exceed 0.5. To test this, we first subtracted 0.5 from the relatedness estimates between the paired males and each of their within-pair offspring. We then tested for a significant intercept effect using two intercept-only LMMs, one for *r*_{QG} and the other for *r*_{LR}. We included ‘Brood ID’ as a random intercept in each model. Next, if paired females were related to their extra-pair mates, then we would expect female relatedness to their extra-pair offspring to exceed 0.5. To test this, we subtracted 0.5 from the relatedness estimates between the paired females and each of their extra-pair offspring. Again, we fit two intercept-only LMMs. Here, we included the identity of each offspring’s genetic sire as a random intercept nested within Brood ID. Finally, if paired males were related to their cuckolders, then we would expect their relatedness to extra-pair fry to exceed 0 *after* accounting for any deviation from random allele sharing occurring within their pair bond. Here, we subtracted ½*R* from each relatedness estimate between the paired males and the extra-pair offspring in their brood, where *R* is the relatedness estimate between the paired male and female. Note that R could be positive or negative depending on whether the social pair shared more or fewer alleles than expected by chance. We then fit two intercept-only LMMs, similarly including the identity of each offspring’s genetic sire as a random intercept nested within Brood ID. We omitted broods of 0% paternity for these analyses for the same reasons outlined above.

### Is there spatial structure in the relatedness between males in the field?

Bose et al. [31] recently showed that the vast majority of cuckoldry in the *V. moorii* system is perpetrated by unpaired males in the population. We therefore tested whether unpaired males in the population were caught in close proximity to any of their male relatives that were pair bonded. Our sampling regime during the dry season (October sampling excursion) also included the opportunistic capture of 40 unpaired males in addition to the paired males that we sampled. Unpaired males that were observed to be in close proximity to certain territories were captured opportunistically, fin-clipped, and genotyped. We used the distances that we measured between territories within the study quadrat and a field-drawn sketch of nest locations to calculate the relative locations of all nests using trigonometric principles. From these locations, we generated a matrix describing the pairwise distances between territories and between males (in Mathematica version 11.2, Wolfram Research, Inc.). For the purposes of this analysis, the spatial positions of our unpaired males were assumed to be the same as the territories that they were captured next to. Next, we generated another matrix to describe the pairwise relatedness estimates, *r*_{QG} and *r*_{LR}, between each unpaired male and each paired male in the study quadrat. In order to test whether unpaired males remain in close proximity to paired male relatives, we conducted a series of permutation tests. First, we took the maximum relatedness estimate (*r*_{QG} and *r*_{LR}) that we could find between each unpaired male and the paired males within an X m radius. We then subtracted the maximum relatedness estimate between each unpaired male and the paired males *beyond* the X m radius. We took the average of these Δ*r* values as our observed test statistic. We then permuted the spatial positions of unpaired males in our sample by randomly assigning each unpaired male to a territory within the quadrat. In each permutation round, territories were able to receive up to one unpaired male, if we had originally sampled no unpaired males next to them, or up to the number of unpaired males that we had originally captured by them. We then calculated Δ*r* at radius *X* m around each unpaired male as described above and repeated this process 10,000 times. We computed *p* values as the proportion of randomized trials yielding Δ*r* values more extreme than our observed data. We performed this test starting at *X* = 1 m and then repeated it every 1 m increment up to a radius of 10 m. Note that for two territories, the paired male was replaced by a new paired male during our sampling period. Because we could not know the exact date when this switch occurred in relation to when we captured most of the unpaired males in the quadrat, we opted for a conservative approach and omitted these two paired males as well as the four unpaired males that were captured alongside them from these analyses. Thus, these analyses were conducted on 36 unpaired males and 72 breeding territories/paired males that we had genotyped. This omission does not qualitatively change our results. Because of the multiple comparisons here, we used the Benjamini–Hochberg procedure and set the false discovery rate to 10% [77].

Lastly, we generated relatedness matrices describing the pairwise relatedness estimates, *r*_{QG} and *r*_{LR}, between all paired males in the study quadrat during the dry season. Using these matrices as well as our distance matrix, we tested for any spatial structure in relatedness between the paired males in our study quadrat (i.e. increasing or decreasing relatedness with spatial separation between territories). For this, we used Mantel tests (vegan r package, ref. [79]), each based on 10,000 random permutations of the data, and we performed separate Mantel tests for *r*_{QG} and *r*_{LR} estimates.