Genetic relatedness in social groups of the emerald coral goby Paragobiodon xanthosoma creates potential for weak kin selection

Animals forming social groups that include breeders and nonbreeders present evolutionary paradoxes; why do breeders tolerate nonbreeders? And why do nonbreeders tolerate their situation? Both paradoxes are often explained with kin selection. Kin selection is, however, assumed to play little or no role in social group formation of marine organisms with dispersive larval phases. Yet, in some marine organisms, recent evidence suggests small‐scale patterns of relatedness, meaning that this assumption must always be tested. Here, we investigated the genetic relatedness of social groups of the emerald coral goby, Paragobiodon xanthosoma. We genotyped 73 individuals from 16 groups in Kimbe Bay, Papua New Guinea, at 20 microsatellite loci and estimated pairwise relatedness among all individuals. We found that estimated pairwise relatedness among individuals within groups was significantly higher than the pairwise relatedness among individuals from the same reef, and pairwise relatedness among individuals from the same reef was significantly higher than the pairwise relatedness among individuals from different reefs. This spatial signature suggests that there may be very limited dispersal in this species. The slightly positive relatedness within groups creates the potential for weak kin selection, which may help to resolve the paradox of why breeders tolerate subordinates in P. xanthosoma. The other paradox, why nonbreeders tolerate their situation, is better explained by alternative hypotheses such as territory inheritance, and ecological and social constraints. We show that even in marine animals with dispersive larval phases, kin selection needs to be considered to explain the evolution of complex social groups.

This central explanation was complemented by other reasons for nonbreeders to remain peacefully in the group, such as that nonbreeders could wait to inherit breeding positions within the group or on a nearby breeding territory (future selection hypothesis, Kokko & Johnstone, 1999;Williams, 1966;Woolfenden & Fitzpatrick, 1978) and that nonbreeders stay within a group rather than leave to breed elsewhere because of ecological constraints (ecological constraints hypothesis, Emlen, 1982), or wait peacefully within a group rather than contest to breed due to social constraints (social constraints hypothesis, Koenig & Pitelka, 1979, Cant et al., 2001. Most of the evidence for these ideas comes from cooperatively breeding birds (e.g., Emlen & Wrege, 1988;Koenig & Pitelka, 1979;Pruett-Jones & Lewis, 1990;Woolfenden & Fitzpatrick, 1978Komdeur, 1992), mammals (e.g., Cant et al., 2001;Clutton-Brock, 2002;Lukas & Clutton-Brock, 2013) and eusocial insects (e.g., Field et al., 1999;Korb, 2009;Ross & Keller, 1995). Other taxa with social group structures remain understudied. Of the four commonly studied drivers of social group formation, kin selection is considered primary and has received the most attention.
For kin selection to operate, individuals must interact with relatives. There are three widely recognized dispersal mechanisms that can lead to relatives finding themselves in the same social group: delayed dispersal, kin associations during dispersal and limited dispersal. Delayed dispersal, when offspring delay the departure from their native group for one or several breeding seasons, is well known in birds and mammals (see reviews in Emlen, 1994;Stacey & Koenig, 1990). Delayed dispersal leads to social groups including parents and their offspring, which can secondarily lead to nonbreeding subordinates providing alloparental care and receiving kin selection benefits (Koenig et al., 1992). Another way for dispersal patterns to lead to social groups made up of relatives is kin association during dispersal, when, for example, siblings leave the natal territory together and travel part or all of the dispersal distance together to settle close to one another. Small kin groups staying together during dispersal has been demonstrated in social insects, birds and mammals (Metheny et al., 2008;Peeters & Ito, 2001;Sharp et al., 2008;Williams & Rabenold, 2005), and it has been shown to have the potential to lead to cooperative kin groups, as it not only leads to high relatedness within groups but also reduced local competition among kin (Kümmerli et al., 2009). Limited dispersal (or population viscosity) occurs when organisms only travel short distances from their place of birth or when they return to their place of birth after the dispersal period (Queller, 1992(Queller, , 1994. This pattern has been documented in a wide range of organisms including lichen, arthropods, gastropods and vertebrates (e.g., Chapuisat et al., 1997;Grosberg, 1987;Huyvaert & Anderson, 2004;Prince et al., 1987;Walser, 2004). Limited dispersal increases the probability of relatives interacting with one another and can thus shape group structures without relying on kin discrimination (Aguillon et al., 2017;Cornwallis et al., 2009;Finn et al., 2006). Different dispersal patterns may lead to close relatives within the same groups or other small-scale relatedness patterns, yet for some taxa this relationship has yet to be disentangled.
While the kin selection hypothesis has generally been overlooked in marine species based on the assumption that the dispersal larval phase breaks up kin aggregations (Leis, 1991;Shanks, 2009;Victor, 1984), recent empirical and theoretical results suggest that delayed dispersal, kin cohesion and/or limited dispersal may lead to some close relatives finding themselves in close proximity even in marine populations (D'Aloia & Neubert, 2018;. In the marine environment, evidence of delayed dispersal comes from very few species that do not have a larval dispersal phase, such as the damselfish Acanthochromis polyacanthus (Miller-Sims et al., 2008), and/or are organized in eusocial societies, such as snapping shrimp of the genus Synalpheus (Duffy, 1996). Support for the kin cohesion hypothesis is found in a number of examples of full and half siblings near each other in marine species (e.g., Bernardi et al., 2012;Buston et al., 2009;Dubé et al., 2020;Herrera et al., 2016;Riquet et al., 2017;Selkoe et al., 2006). Finally, evidence for limited dispersal comes from observations of high levels of self-recruitment to local populations for a range of species (e.g., Almany et al., 2007;Beldade et al., 2012;D'Aloia et al., 2013;Jones et al., 1999Jones et al., , 2005Rueger et al., 2020;Saenz-Agudelo et al., 2012;Selwyn et al., 2016;Swearer et al., 1999) and from dispersal curves showing that the probability of successful dispersal declines as a function of distance (Almany et al., 2013(Almany et al., , 2017Buston et al., 2012;D'Aloia et al., 2015;Williamson et al., 2016). These patterns raise the possibility that relatedness may after all be a factor in the formation of social groups in marine animals. Despite its potential importance, for most marine taxa we do not know whether groups of individuals include close relatives.
A suitable model organism for assessing the potential for genetic relatedness and kin selection in marine fishes is the emerald coral goby Paragobiodon xanthosoma (Bleeker, 1853). It is an obligate coral-dwelling fish found throughout the Indo-Pacific, where it exclusively inhabits colonies of the needle coral Seriotopora hystrix (Dana 1846, Lassig, 1976. Obligate coral-dwelling fishes often have mutualistic relationships with their host cnidarian, providing nutrients and defence against corallivores and encroaching algae (Dixson & Hay, 2012;Dirnwoeber & Herler, 2013;Holbrook et al., 2008), and in turn receiving protection and suitable breeding sites (Lassig, 1976). P. xanthosoma lives in social groups of two to 16 individuals composed of one monogamous breeding pair and zero to 14 nonbreeding subordinates (Wong et al., , 2008a. They are organized in a strict size hierarchy with the male being the largest individual, the female the second largest and the subordinates progressively getting smaller (Lassig, 1976;Wong et al., 2007). Rank ascension only occurs when the larger individual disappears, as lower ranks that get too close in size are evicted . Therefore, subordinates adjust their growth to remain at a smaller size and wait to ascend rank to eventually inherit the breeding territory . Most of the parental care in P. xanthosoma in the form of fanning and cleaning the demersal brood is performed by the male (Lassig, 1976;Wong et al., 2008b). After hatching, Paragobiodon spp. have a pelagic larval duration of 36-47 days (Brothers et al., 1983), before they settle into coral colonies. Little is known about what happens to coral gobies during the pelagic larval phase or how far they disperse during that time. However, observations of limited dispersal and kin aggregations in other Gobiiformes (D'Aloia et al., 2015;Rueger et al., 2020;Selwyn et al., 2016) make it one of the most likely candidates for the presence of relatives in social groups in a marine fish, and thus for kin selection to potentially play a role in the formation of social groups.
To study whether genetic relatedness and kin selection might play a role in the social evolution of coral gobies, we test two specific hypotheses: (i) pairwise relatedness is higher within groups than among groups, and higher within reefs than among reefs; and (ii) dyad relatedness is dependent on the nature of the relationship (breederbreeder; breeder-nonbreeder; nonbreeder-nonbreeder). Elevated relatedness is predicted within groups and within reefs based on any of the three dispersal mechanisms. Evidence for these mechanisms is based on the status and size of the relevant dyad members, with the caveat that size is not a perfect predictor of age in P. xanthosoma due to the size hierarchies in which growth is regulated. Elevated relatedness is predicted for breeder-nonbreeder pairs under delayed dispersal, because parent-offspring dyads may occur. Elevated relatedness is predicted for breeder-breeder, nonbreeder-nonbreeder and similarly sized pairs under sibling cohesion, because here high relatedness might be found among individuals that are close in age.
Elevated relatedness is predicted irrespective of the nature of the relationship under limited dispersal with overlapping generations, as offspring leave groups, but higher order relatives, such as grandoffspring, nieces and nephews, may return.

| Sampling
Field sampling was conducted on inshore reefs near Mahonia Na Dari Research and Conservation Centre, Kimbe Bay, Papua New Guinea (5°30′S, 150°05′E), in May 2019. We collected tissue samples from 78 individuals in 16 groups on three inshore reefs. Groups were collected at "Bob's Knob" (reef 1, N groups = 10, N individuals = 39), "Kilu" (reef 2, N groups = 2, N individuals = 19) and "Lui" (reef 3, N groups = 4, N individuals = 20). Of these, reef 2 lies farthest to the north in Kimbe Bay and is ~560 m offshore, reef 3 is located ~3.5 km farther south and ~630 m offshore, and reef 1 is ~1.5 km farther south and ~730 m offshore. No direct information about ocean currents connecting these reefs is available, but larval connectivity among inshore reefs has been found in other reef fish within the same area (Rueger et al., 2020). Groups comprised 2-12 individuals (mean group size ± SE = 4.69 ± 0.62); two groups had N = 2 group members, two groups N = 3, five groups N = 4, three groups N = 5, three groups N = 7, and one group N = 12. All individuals in a group were caught on SCUBA using hand nets and diluted clove oil solution as a mild anaesthetic (Munday & Wilson, 1997). Groups were clearly distinguishable because they live in one distinct coral head. Coral heads were thoroughly checked upon collection and again before returning the sampled fish to the coral, and to the best of our knowledge all group members were sampled. The standard length (SL) of each fish was measured to the nearest 1.0 mm using calipers underwater, and a fin clip was taken from the caudal fin. Maturity and breeder status were assessed according to rank, which was based on SL: The two largest individuals represent the monogamous breeding pair (mean SL ± SE, SL R1 = 18.44 ± 0.88 mm; SL R2 = 15.19 ± 0.99 mm), all other individuals are nonbreeders (SL R3 = 10.68 ± 1.33 mm; SL R4 = 7.79 ± 1.24 mm; SL R5 = 7.71 ± 1.61 mm; SL R6 = 8.63 ± 2.21 mm; SL R7 = 7.75 ± 2.36 mm; SL R8 = 10 mm; SL R9 = 8 mm; SL R10 = 7 mm; SL R11 = 5.5 mm; SL R12 = 5 mm) (Wong et al., 2008b). Tissue samples were preserved in 99% ethanol for genetic analysis. After sampling, the fish were released back into their host coral. Contigs were scanned for microsatellite repeats (minimum perfect repeat length =5) and primers were designed with msatcommander (version 1.08-beta) software. Primer lengths ranged between 18 and 24 bp, primer Tms between 57°C and 62°C, and optimal PCR (polymerase chain reaction) product size was 325-350 bp. Over 4000 unique tetrameric microsatellite loci with designed primers were recovered from the assembled contigs.

| Marker development and multiplex PCR
All individuals were sequenced at 54 microsatellite loci using a multiplex PCR protocol for targeted amplicon sequencing. For detailed description of multiplex PCR and Nextera barcoding methods, see D' Aloia et al., (2017). We amplified the loci with multiplex PCRs using QIAGEN Multiplex kits and the primers are listed in Table S1.
Samples were pooled across multiplexes and Illumina's S5 and N7 Nextera primers were used to run a barcoding PCR. The sequencing library was prepared by pooling all barcoded individuals, which were then size-selected with Ampure XP (Beckman Coulter). The library was sequenced on an Illumina MiSeq with paired 250-bp reads at the BioResource Center, Cornell University.

| Data processing
We used a python script (amplicon.py, https://bitbu cket.org/corne ll_bioin forma tics/ampli con/src/maste r/) to call genotypes at each microsatellite locus. Default commands were used except the following: -c 1, -a 0.005, -l 150. We also explored two reads ratios (-r command) for calling heterozygotes: the default of -r 20 as well as -r 40. A minimum of two reads were required for each allele; otherwise the diploid genotype was recoded as missing data. To retain only the highest quality markers and individuals, we first excluded loci missing >20% of individuals, then individuals with >20% missing loci. After the filters were applied, 73 individuals remained, analysed at 37 loci. Seventeen of these loci were in Hardy-Weinberg equilibrium (HWE) at -r 20 and 20 of these loci were in HWE at -r 40.
TA B L E 1 Polymorphic microsatellite markers developed for Paragobiodon xanthosoma: primer sequence, repeat motif and repeat count, number of alleles (Na), observed (Ho) and expected (H E ) heterozygosities, deviation from the exact test of Hardy-Weinberg-Equilibrium (HWE), and fixation index (F IS ) are presented for each locus Statistical analysis (see below) was conducted using relatedness values based on our final data set. The analyses were repeated to test for robustness with the complete marker set (37 markers including those out of HWE), as well as for each reef separately using only those markers in HWE for the samples on each reef (reef 1, N = 20 markers; reef 2, N = 32 markers; reef 3, N = 30 markers). The results of these supplemental analyses can be found in the Supporting Information part III (Tables S1, S2, Figures S3-S11).

| Population structure
Before investigating relatedness, we investigated whether there was any population structure, because relatedness of any two individuals is the probability that they share an allele relative to the probability that two individuals from the same population share the allele (Queller, 1994). We estimated pairwise F ST between different reefs using genepop (Raymond & Rousset, 1995;Rousset, 2008).

| Relatedness
Pairwise relatedness was calculated for all pairs of individuals (hereafter called dyads) in the population using the R package Demerelate (Kraemer & Gerlach, 2017). A simulation to compare different relatedness estimators was conducted in the R package related (Csillery et al., 2006;Pew et al., 2015). Individuals of known relatedness were simulated (1000 each of full siblings, half siblings, parent-offspring and unrelated dyads) and estimated relatedness values were com-  Figure S1). Given this, we decided to perform dyad relatedness calculations using the Queller-Goodnight estimator because it has recently been used in similar species and it allows comparison with other reef fish species where relatedness has been investigated (Buston et al., , 2009Rueger et al., 2020).
Overall mean degree of relatedness in the sample was calculated as the average of pairwise relatedness among all sampled individuals.
Mean degree of relatedness among group members was calculated as the average of the pairwise relatedness within each group.

| Statistical analysis
We used R (R Core Team, 2019) and lme4 (Bates et al., 2012) to perform linear mixed model analysis (LMM). Pairwise relatedness (Queller-Goodnight estimator) was used as the response variable.
Individual tags of both dyad members were used as crossed random terms (random intercept with fixed mean, formula: ..+1|individual 1 tag +1|individual 2 tag) to account for nonindependence of relatedness values involving the same individuals. We used Akaike's information criterion (AIC) (Sakamoto et al., 1986) for model selection. If ΔAIC was below 2, a likelihood ratio test was performed to decide which model was the best fit, with preference for the simpler model.

Significance tests for LMM were performed by likelihood ratio tests
of the full model with the effect in question against the model without the effect. No obvious deviations from normality and homoscedasticity were detected by visually inspecting the residual plots. Models were also tested for outliers and collinearity between variables with performance (Lüdecke et al., 2020). No outliers or high variance inflation factors (VIFs) were detected in any of the best-fit models.
Conditional and marginal R 2 were calculated using Nakagawa's R 2 in performance (Nakagawa & Schielzeth, 2013;Nakagawa, Johnson & Schielzeth, 2017). We performed post-hoc tests using the emmeans package (Lenth, 2020). p-Values were adjusted for multiple testing using a multivariate t distribution correction.
To test the hypothesis that pairwise relatedness estimates are elevated within groups, location of individuals in a dyad was added to the LMM as a fixed factor with three categories: individuals in the same group; individuals on the same reef but not in the same group; or individuals on separate reefs. To test the hypothesis that pairwise relatedness estimates differed with dyad composition, dyad maturity (categorical variable, breeder-breeder; breeder-nonbreeder; nonbreeder-nonbreeder) and dyad size ratio (continuous variable, larger individual: smaller individual) were added as fixed factors.

| Relatedness
The prediction that pairwise relatedness should be dependent on the spatial relationships of the dyads (higher within groups and reefs than among groups and reefs) was supported by the study.
Observed mean relatedness in the overall sample population was −0.010 (±0.001 SE). There was a significant difference in relatedness based on dyad location (same group, same reef but different group, and different reef) (χ 2 (2) =50.911, p < .001), with estimated mean relatedness within groups being 0. for dyads on the same reef, compared to dyads on different reefs (CI =0.008, 0.022, p < .001) (Figure 1).
The prediction that dyad relatedness will be dependent on the nature of the relationship (pairing or size ratio) was not supported by the study. There was no significant difference in pairwise relatedness between dyad type (breeder/breeder, breeder/nonbreeder, nonbreeder/nonbreeder) (χ 2 (2) = 3.356, p = .187), but the best fit model did include dyad type. In contrast, the best fit model did not include size ratio ( Figure S2). The best fit model probably included dyad type because the interaction between location and pairing was marginally significant (χ 2 (4) = 9.111, p = .058) (Figure 2). Post-hoc contrasts reveal that the breeder/breeder dyads are more closely related when they are in the same group, compared to when they are on the same reef but in different groups, and to when they are on separate reefs from one another (Figure 3). Similarly, nonbreeder/ nonbreeder and breeder/nonbreeder dyads show the highest relatedness estimates when they are in the same group, although the differences are smaller (Figures 2 and 3). The two fixed factors (location and pairing) together explained ~2% of the variance in relatedness estimates (R 2 c = .231, R 2 m = .022) and the variance of the random intercept was 0.0005 (±0.022 SD) and 0.0002 (±0.013 SD) for the two individual tags per dyad.

| DISCUSS ION
Assessing small-scale relatedness patterns and their connection to dispersal patterns is crucial for our understanding of group formation in marine systems. Here, we found that pairwise relatedness was fractionally higher within groups than within reefs, and fractionally higher within reefs than within the population at large. We found that relatedness did not differ with dyad breeding status or size ratio in groups. These results are consistent with the predictions of the limited dispersal hypothesis. The elevated relatedness within groups may indicate an opportunity for weak kin selection to operate in Paragobiodon xanthosoma.
Limited dispersal explains the relatedness patterns found in P.
xanthosoma. Similar to other reef fish species, we found relatedness to decline with larger spatial scales, indicating that dispersal success in this species attenuates with distance. Dispersal was not limited enough, however, to lead to isolation-by-distance patterns between the different reefs in our population. For the pajama cardinalfish Sphaeramia nematoptera very similar small-scale relatedness patterns were found and, using parentage analysis, more recruitment back to the natal group and reef was confirmed to be the mechanism causing these patterns (Rueger et al., 2020). Another strong example for limited dispersal operating in a small coral reef fish is the sponge-dwelling goby Elacatinus lori; for this species dispersal declined exponentially with distance and median dispersal distance was just 1.7 km (D'Aloia et al., 2015). In contrast to the strong evidence for limited dispersal, we found no evidence for delayed dispersal or kin cohesion operating in P. xanthosoma: relatedness was not higher in breeder-nonbreeder dyads, as would be predicted by the delayed dispersal hypotheses (e.g., Duffy, 1996); relatedness was not higher in similarly sized dyads, as would be predicted by the sibling cohesion hypothesis (e.g., Buston et al., 2009). Based on growing evidence from coral reef fish populations, limited dispersal is often sufficient to explain relatedness patterns occurring at small spatial scales. How might small differences in relatedness among groups or populations contribute to social group formation in coral reef fishes?
From the perspective of subordinate P. xanthosoma, kin selection is unlikely to play a large role in answering why they forego reproduction. To date, there is no evidence that subordinate P. xanthosoma contribute to the reproductive output of the dominants, consistent with another reef fish with similar social organization (Buston, 2004;Buston & Elith, 2011), which means there is nothing they gain by being related to the dominant breeders in the group. Rather, natural selection seems to have favoured subordinates that forego reproduction due to a combination of alternative factors: they stand to inherit the breeding territory ; ecological constraints keep them from dispersing (Wong, 2010); and social constraints keep them from contesting for a breeding position within the group .
However, from the perspective of dominant P. xanthosoma, small levels of relatedness may be sufficient to confer enough benefits to help explain why they tolerate subordinate nonbreeders. The within-group relatedness recorded here indicates that dominants may gain inclusive fitness by accepting nonbreeding group members as long as they remain small and inflict no costs. Both P. xanthosoma and A. percula subordinates are known to regulate their growth to remain under a size threshold to be tolerated in the group (Buston, 2003;Wong et al., 2007) and, at least in A. percula, subordinates that did this were found to inflict no fitness costs on dominants (Buston, 2004;Wong et al., 2016). In the absence of costs inflicted by the subordinates, the estimated within-group relatedness of 0.025 may be high enough to confer kin-selected benefits to the dominant breeders as their distant relatives go on to inherit their territories and breed when they die. While relatedness values tend to be low, small-scale relatedness patterns at the group level may indicate that weak kin selection plays a role in social group formation of some coral reef fishes.
Is an average pairwise relatedness of 0.025 sufficient for behaviours to evolve via kin selection? According to Hamilton's rule, relatedness does not have to be high for social behaviours to evolve (Foster et al., 2006;Hamilton, 1964). Small increases in fitness can be profitable as long as the ecological cost does not become too great and therefore selection could favour evolution of behaviours based on small elevations in relatedness, as long as they are above the population average (West-Eberhard, 1975). In the cichlid Variabilichromis moori, inclusive fitness benefits conferred by an average relatedness of 0.034 between male breeders and cuckolding males helps to mitigate the cost of cuckoldry (Bose et al., 2019). Even in some mammals performing alloparental care, mean within-group relatedness was found to be low; the mean relatedness ± SE of females in groups of wild boar, Sus srofa, was 0.08 ± 0.1 (Briga et al., 2012;Poteaux et al., 2009), and in spear-nosed bats, Phyllostomus hastatus, was only 0.01 ± 0.01 (Briga et al., 2012;McCracken & Bradbury, 1981).
West-Eberhard (1975) laid out three scenarios in which the benefit to cost ratio may be large enough to make altruism beneficial even at low values of r: when the beneficiary has a lot to gain; when the aid provided is cheap; and when small amounts of aid have large effects.
In the case of P. xanthosoma, the "beneficiary," the subordinate nonbreeder has a lot to gain, that is safety and the potential inheritance of the breeding territory (Wong, 2010;Wong et al., 2007) and the "donor," the dominant breeder, is providing aid that is cheap, that is tolerating the subordinate on their territory as long as it remains small and inconsequential (Wong et al., 2008a). Thus, the fractional elevation in relatedness within groups of some species may help fill a large gap in our understanding of social group formation.
Other mechanisms are likely to contribute to the evolution of complex groups in P. xanthosoma, but our study adds to a growing body of evidence that illuminates the importance of evaluating relatedness patterns when investigating group structures and the evolution of mating and social systems of marine animals. While delayed dispersal is rare in marine fishes and has not been found to lead to stable groups of first-degree relatives, limited dispersal may lead to relatedness values high enough to confer subtle kin selection benefits and help explain social group formation in marine environment. F I G U R E 2 Pairwise relatedness estimates for nine types of dyads of Paragobiodon xanthosoma in Kimbe Bay. Dyads are categorized based on location of individuals within the dyad (different reef = light green; same reef but different group = turquoise; same group = blue) and relative maturity (large silhouette: breeder; small silhouette: nonbreeder). Central bar represents the median; boxes represent lower and upper quartiles; whiskers represent ±1.5 interquartile range [Colour figure can be viewed at wileyonlinelibrary.com] Even in marine animals with larval dispersal phases, kin selection needs to be considered as a possible contributing factor to the evolution of complex groups.

ACK N OWLED G EM ENTS
We thank the Tamare-Kilu communities for granting access to their reefs for this study. We thank Siobhan Heatwhole, Catherine Froehlich and Nelson Sikatura for field assistance, Chancey MacDonald for help with the artwork and three anonymous reviewers for suggestions that helped to improve the manuscript. Fieldwork