http://orcid.org/0000-0002-9100-217X
Figures
Abstract
The phenomenon of chaotic genetic patchiness is a pattern commonly seen in marine organisms, particularly those with demersal adults and pelagic larvae. This pattern is usually associated with sweepstakes recruitment and variable reproductive success. Here we investigate the biological underpinnings of this pattern in a species of marine goby Coryphopterus personatus. We find that populations of this species show tell-tale signs of chaotic genetic patchiness including: small, but significant, differences in genetic structure over short distances; a non-equilibrium or “chaotic” pattern of differentiation among locations in space; and within locus, within population deviations from the expectations of Hardy-Weinberg equilibrium (HWE). We show that despite having a pelagic larval stage, and a wide distribution across Caribbean coral reefs, this species forms groups of highly related individuals at small spatial scales (<10 metres). These spatially clustered family groups cause the observed deviations from HWE and local population differentiation, a finding that is rarely demonstrated, but could be more common than previously thought.
Citation: Selwyn JD, Hogan JD, Downey-Wall AM, Gurski LM, Portnoy DS, Heath DD (2016) Kin-Aggregations Explain Chaotic Genetic Patchiness, a Commonly Observed Genetic Pattern, in a Marine Fish. PLoS ONE 11(4): e0153381. https://doi.org/10.1371/journal.pone.0153381
Editor: Heather M. Patterson, Department of Agriculture and Water Resources, AUSTRALIA
Received: January 25, 2016; Accepted: March 29, 2016; Published: April 27, 2016
Copyright: © 2016 Selwyn et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: All relevant data are within the paper and its Supporting Information files.
Funding: This work was supported by the Natural Sciences and Engineering Research Council of Canada CRO grant# 227965 awarded to P.F. Sale, a Global Environment Facility Coral Reef Targeted Research project and by a Natural Sciences and Engineering Research Council Tier I Canada Research Chair grant to DDH; a Great Lakes Institute for Environmental Research post-doctoral fellowship from the University of Windsor and Texas A&M University – Corpus Christi start-up funds to JDH.
Competing interests: The authors have declared that no competing interests exist.
Introduction
Marine organisms with dispersive pelagic larvae are expected to be characterized by little to no genetic differentiation among populations over large geographic areas due to high gene flow. However, many marine species exhibit slight, yet significant, levels of genetic heterogeneity over various local spatial scales [1,2] which may be temporally unstable. Typically in these situations, molecular markers deviate from the expectations of Hardy-Weinberg equilibrium (HWE) within some samples [3–5]. This complex pattern, termed chaotic genetic patchiness (CGP), has been attributed to non-equilibrium conditions caused by variation in reproductive success of breeding adults during larval recruitment [2,6]. By seeking to understand the biological mechanisms underlying these patterns we can deepen our understanding of the ecology and evolution of complex marine populations.
There are several mechanisms that could create CGP across populations of marine organisms [6] including: local habitat differences [1], temporal variability in ocean currents affecting dispersal [7], selection at settlement [8], differential reproductive success [9], genetic drift within isolated populations [2], and the formation of kin-aggregations [2]. In marine populations, kin-aggregations are often overlooked as a mechanism to explain CGP because the probability of sampling relatives is expected to be low for species with highly dispersive pelagic larvae [10]. However, recent studies show that larvae can remain in kin-aggregations in the plankton [7,11,12]. Kin-aggregations likely form during the larval stage through aggregations of locally-spawned eggs and larvae by ocean currents and patchiness of food resources [13]. Various life-history traits, including demersal spawning and shoaling, might also cause kin-aggregation formation in fishes [13].
Here we investigate the population genetic structure of a shoaling marine goby, Coryphopterus personatus. Despite being wide-spread throughout the Caribbean little information exists about its life history, particularly information about larval duration, dispersal, and population structure. We seek to use observed patterns of genetic structure to test the kin-aggregation hypothesis and discuss the consequences of this mechanism for the ecology of this and similar species. We then propose competing alternative hypotheses for the mechanism leading to the formation of kin-aggregations.
Methods
Fish were collected during summer 2002 from nine sites in the Mesoamerican barrier reef system (34 to 85 per site, Fig 1). Within each site, multiple shoals of C. personatus were collected by divers from a 1-ha area on the fore-reef and pooled for storage in 95% ethanol. All animals were humanely euthanized in the field using a standard fish sedative, Tricaine mesylate (MS222). Live fish were placed in a seawater bath with MS222 added for 5 minutes or until breathing ceased. All collections were approved by and performed in accordance with the ethical guidelines of the University of Windsor and the Belize fisheries department, the Mexican Secretaría de Medio Ambiente y Recursos Naturales, and the Department of Natural Resources and Environment of Honduras. Genomic DNA was extracted from pectoral fin tissue of each individual following the silica-based 96-well plate protocol of Elphinstone et al.[14]. Ten microsatellite loci were chosen from Hepburn et al. [15] and Hogan et al. [16]. PCR amplification was then performed in 12.5 µL reactions comprised of: approximately 100 ng template DNA, 200 μM of forward dye-label primer and reverse primer, 200 µM of each dNTP, 0.1 U Taq polymerase (0.025 U for CPER 184, Invitrogen, Burlington, Canada), 1x PCR buffer (provided by the manufacturer) and locus specific concentrations of MgCl2 and bovine serum albumin (Table 1). PCR conditions were 94˚C for 2 minutes, followed by locus specific numbers of cycles of 94˚C for 15 s, locus specific annealing temperatures and times (Table 1), with a final extension of 72˚C for 90 s. The size of amplicons was determined using a LiCor 4300 DNA Analyzer with GeneImagIR 4.05 software (Scanalytics, Inc).
Insets show scatter plots (and density plots in the case of two clusters) of clusters from DAPC analysis within sampling locations. The axes of the plots are the first two discriminant functions used to delineate clusters with inertia ellipses representing 67% of the variance. The end of the lines connected to the centre of each inertia ellipse represent individuals plotted on each discriminant function and denotes cluster membership. In locations with only two clusters present there is only one discriminant function, as such density plots of proportion of individuals present at each value of the discriminant function were included to show cluster separation. Numbers in parentheses indicate how many individuals were collected from each location.
One of the common hallmarks of chaotic genetic patchiness is many deviations from Hardy-Weinberg equilibrium combined with weak but significant genetic differentiation between samples. To test for deviations from HWE exact tests for goodness of fit to the expectations of HWE were performed in adegenet [17]. To test the possibility of deviations from HWE being due to elevated levels of inbreeding, Monte Carlo tests for homozygote excess, which is indicative of either high levels of inbreeding, or the presence of null alleles, were performed using the U-score implemented in HWxtest [18]. Inbreeding coefficients were calculated for each site using gstudio [19]. Both inbreeding and null alleles are supported by homozygote excess in a population. Since individuals homozygous for a null allele or heterozygous for two null alleles will present as missing data, there may be an association between the amount of missing data at a locus, in a population and deviation from HWE when null alleles are present. Therefore, we performed a linear regression between the absolute value of the difference between expected and observed heterozygosity and the proportion of individuals which failed to amplify at each locus by sample combination to determine if the observed patterns are consistent with null alleles. To test for genetic differentiation among samples the proportion of different alleles (PD) and pairwise FST between sites were calculated across all loci using adegenet [17]. A Mantel test was performed in ade4 to test for isolation-by-distance among sites using PD and FST [20]. Significance was assessed by permuting genotypes among samples 10,000 times, using sequential Bonferroni to correct for multiple testing.
To test the hypothesis that increased levels of relatedness within sites explains the pattern of chaotic genetic patchiness observed here pairwise relatedness (r) was calculated between all pairs of individuals using the triadic likelihood estimator, which has been shown to be the least biased relatedness estimator in many circumstances, using the R package related [21,22]. The arithmetic mean of the pairwise relatedness () was calculated for each geographic sample in R v3.1.1 [23]. Significance of mean within sample relatedness was calculated using a one-tail permutation test with genotypes permuted among geographic samples or clusters 1,000 times.
Our initial sampling pooled multiple shoals of individuals; if sites were composed of multiple groups of related individuals (i.e., shoals) and if those groups were pooled during collection, then the mean pairwise relatedness at the site level would be artificially deflated. We performed a Discriminant Analysis of Principal Components (DAPC) to identify genetic clusters of individuals within sites as implemented in adegenet [17,24]. All possible clustering solutions were compared using the k-means clustering algorithm. The number of clusters within each site was found by evaluating all possible values for k and choosing the most likely based on minimum BIC. The mean pairwise relatedness () was then recalculated for each cluster within sites generated with the DAPC and significance determined as above.
As further evidence of elevated relatedness and to minimize error associated with continuous metrics of relatedness, COLONY was used to determine the number of full- and half-sib dyads within clusters using the pair-likelihood score algorithm allowing for inbreeding [25]. We assumed polygamy and polyandry and a 1% genotyping error rate. Sibship was considered when dyads were supported with ≥95% confidence. Evidence of high levels of sibling pairs (both full and half-sibs) within clusters was assessed using a two-tailed permutation test, permuting individuals among clusters 10,000 times. The probability of excluding a group of three unrelated individuals from a full sibship was calculated using code written by the authors [26]. This probability of exclusion is useful in determining the statistical power the set of microsatellite markers used in this study have in determining sibship relationships.
In order to determine if elevated relatedness observed in clusters is the result of multiple dyads of related individuals or a few large groups of relatives within clusters undirected networks were created for each cluster identified using DAPC with individual C. personatus represented as nodes and sibling relationships (full- and half-sibs) represented as edges [27]. The mean local transitivity, or clustering coefficient, was used to understand the degree of interconnectedness within each cluster. Higher transitivity values indicate more interconnected networks of sibship relationships. Transitivity is a metric used to understand networks and ranges from 0 to 1 with 0 representing no clustering of relationships and 1 representing maximally clustered graphs. Biologically, transitivity values near 0 indicate the presence of multiple related pairs of individuals, which are unrelated to others in the cluster. Transitivity values close to 1 indicate the presence of a few large groups of related individuals. Significance of transitivity within each cluster was assessed by simulating 10,000 random graphs using the Erdos-Renyi model with the same number of vertices (individuals) and edges (sibling connections) as the observed cluster [27,28]. Networks were analysed using the package igraph [29].
Results
There were deviations from the expectations of HWE in 48% of site by locus comparisons (S1 Fig). These deviations from HWE were attributed to homozygote excess which was observed in 53% of comparisons (S2 Fig). The ubiquity of the homozygote excesses in all sites and in all but one locus does not support the hypothesis of null alleles but rather, indicates potentially high levels of inbreeding [30]. Additionally, we found no relationship between locus specific sample deviations from the expectations of HWE and the percent of individuals that failed to amplify at each locus (F(1,88) = 0.0018, p = 0.97, Fig 2). Additionally, all samples had significant estimates of Fis for six or more loci (Table 2). All these results taken together refute the hypothesis of null alleles driving the patterns seen here.
The plotted line is the result of a linear regression (F(1,88) = 0.0018, p = 0.97) with the shaded area indicating the 95% confidence intervals. Blue points are locus by sample combinations which were not significantly different from the expectations of HWE based on exact tests. Red points are locus by sample combinations which are significantly deviated from the expectations of HWE based on exact tests.
Significant differences in PD and FST were seen in 69% and 97% of the pairwise site comparisons respectively, even between closely spaced sites (ex. TA1 –TA2: ~5 km) (Table 3). Mantel tests for an isolation-by-distance pattern of differentiation were not significant, indicating no correlation between genetic and geographic distance (PD, r = 0.22, P = 0.19; Fst, r = 0.31, P = 0.14).
Within sample mean pairwise relatedness was 34–60% higher than the random expectation in four sites (BC1, BC2, TA1 and TA3; p < 0.05, Fig 3A). Within-site clustering analyses found 2 to 5 highly related groups of individuals (8–36 individuals per cluster) within each of the nine geographic samples (Figs 1 and S3). Group relatedness measured within these clusters was found to be 73–314% higher than expected in 68% of all clusters (p <0.05, Fig 3B).
Mean pairwise relatedness () values by geographic sample (A) and cluster (B) with standard errors. The shaded region indicates the area within the 95% confidence intervals calculated using a permutation test with 1,000 iterations. Triangles in 3B indicate proportion of full and half-sibs within the cluster with shaded symbols indicating significantly elevated transitivity when compared to randomly generated networks equivalent to those observed in the clusters.
The probability of excluding an unrelated individual from a sibship group was 0.998 when all loci are used and 0.927 (95% CI 0.738 to 0.992) when using a random combination of five loci (the minimum number of loci amplified in any individual). The proportion of full- and half-sibs within clusters was calculated to be 22.4±13.8% (S.D.; 3.0±4.4% and 19.3±11.5% respectively) and was significantly greater than clusters of random individuals generated by permutation (p <0.05; Table 4; Fig 3B). The cluster networks were found to be composed of a few large groups of siblings and have a mean transitivity of 0.504±0.213 which was significantly greater than graphs of equal order and size generated through the Erdos-Renyi model in all but four clusters (p <0.05; Table 4).
Discussion
Clusters of highly related individuals were found within geographic samples of Coryphopterus personatus (Fig 3B). This explained wide-spread deviations from HWE, genetic differentiation over short distances–as little as 5 km–and the non-equilibrium spatial pattern of genetic differentiation. While we do not have appropriate samples to determine temporal stability of these patterns, all of these results are consistent with CGP. Some geographic samples as a whole showed substantially higher than expected levels of relatedness. However, we found that within each sample there were between 2 and 5 highly related clusters of individuals. These clusters had significantly higher proportions of full- and half-sibs than expected and showed high levels of interconnectedness with a few large groups of related individuals, indicating familial relationships among individuals in these groups. This is unexpected for a species with pelagic larvae where the prevailing paradigm suggests wide dispersal and mixing via physical oceanographic processes and homogenizing the gene pool over large distances.
We propose two mechanisms driving the formation of related groups in this species. First, larvae could remain together in the plankton and settle onto a reef together. There are many potential selective benefits to larvae remaining in a group throughout the pelagic larval phase including predator avoidance [31] and maintenance of position within a food patch [32]. The presence of kin-aggregations may simply reflect the fact that group formation occurs in the egg phase. Consistent with this mechanism, related individuals have been found within a single larval cohort in several other species [7,11]. Given this mechanism, we expect to find individuals of the same cohort within an aggregation to be related to each other, but not related to individuals from other cohorts within the same sample site.
Another mechanism that could explain kin-groups on a coral reef is a lack of larval dispersal. It is possible that larvae do not enter the pelagic environment, remaining on their natal reef. Other marine species that do not have a pelagic larval stage are characterized by low levels of gene flow, frequent population bottlenecks and strong phylogeographic breaks [33]. If C. personatus has lost their pelagic life stage, we expect a similar genetic pattern. Additionally, we would expect to see highly related individuals across cohorts within a single sample site. Natal recruitment increases the chance of settling in a suitable habitat [34] and, given that C. personatus needs to feed 2 days post-hatching [35], it is plausible that rather than entering the plankton where food is sparse and patchy [36] larvae could remain on the reef where food is more abundant [37]. However, given the lack of observed isolation by distance between geographic locations this mechanism seems less likely than the possibility of individuals staying together throughout the planktonic phase.
The formation of kin-aggregations is a potentially important evolutionary process due to the likelihood of very little gene flow on extremely small spatial scales and increased rates of inbreeding, potentially leading to greater vulnerability to extirpation [38]. The benefits of kin-aggregations (protection from predation, kin-selection, etc. [39]) may outweigh the costs of inbreeding. The pattern revealed here demonstrates that benthic kin-aggregations are possible in marine species with pelagic larvae and may be more common than previously expected.
Supporting Information
S1 Fig. Deviations from Hardy-Weinberg Equilibrium shown at each geographic sample (rows) for each locus (columns).
Significant deviations at a particular sample by locus comparison indicated with a black box.
https://doi.org/10.1371/journal.pone.0153381.s001
(TIFF)
S2 Fig. Homozygote Excess shown at each geographic sample (rows) for each locus (columns).
Significant deviations at a particular sample by locus comparison indicated with a black box.
https://doi.org/10.1371/journal.pone.0153381.s002
(TIFF)
S3 Fig. Plots showing the number of clusters compared to BIC for each geographic sample.
The red points indicate the minimum BIC which was then used as the most likely number of clusters present within the site.
https://doi.org/10.1371/journal.pone.0153381.s003
(TIFF)
S1 File. Latitude/longitude coordinates for sample locations.
https://doi.org/10.1371/journal.pone.0153381.s004
(CSV)
S2 File. Microsatellite repeat data for individuals with collection location.
https://doi.org/10.1371/journal.pone.0153381.s005
(CSV)
S1 Table. Table of proportion of individuals which failed to amplify at each site/locus combination.
https://doi.org/10.1371/journal.pone.0153381.s006
(XLSX)
Acknowledgments
We thank K. Pennoyer, S. Furiness, and C. van Slyck for manuscript edits. P. Sale, R. Hepburn, J. Kritzer, S. Ludsin, P. Chittaro, C. Mora and P. Usseglio for field collections. G. Nanninga and two anonymous reviewers for helpful comments to improve the manuscript. The Institute of Marine Studies at the University of Belize for field support. NSF Award 1429518 for high power computing services.
Author Contributions
Conceived and designed the experiments: JDH DDH. Performed the experiments: JDH. Analyzed the data: JDS AMDW LMG JDH DSP. Contributed reagents/materials/analysis tools: JDH DDH. Wrote the paper: JDS JDH AMDW LMG DSP DDH.
References
- 1. Selkoe KA, Watson JR, White C, Horin TB, Iacchei M, Mitarai S, et al. Taking the chaos out of genetic patchiness: seascape genetics reveals ecological and oceanographic drivers of genetic patterns in three temperate reef species. Mol Ecol. 2010;19: 3708–3726. pmid:20723063
- 2. Broquet T, Viard F, Yearsley JM. Genetic Drift and Collective Dispersal Can Result in Chaotic Genetic Patchiness. Evolution. 2013;67: 1660–1675. pmid:23730760
- 3. Hogan JD, Thiessen RJ, Heath DD. Variability in connectivity indicated by chaotic genetic patchiness within and among populations of a marine fish. Mar Ecol Prog Ser. 2010;417: 263–275.
- 4. Marino IAM, Barbisan F, Gennari M, Giomi F, Beltramini M, Bisol PM, et al. Genetic heterogeneity in populations of the Mediterranean shore crab, Carcinus aestuarii (Decapoda, Portunidae), from the Venice Lagoon. Estuar Coast Shelf Sci. 2010;87: 135–144.
- 5. Jolly M, Viard F, Weinmayr G, Gentil F, Thiébaut E, Jollivet D. Does the genetic structure of Pectinaria koreni (Polychaeta: Pectinariidae) conform to a source–sink metapopulation model at the scale of the Baie de Seine? Helgol Mar Res. 2003;56: 238–246.
- 6. Larson RJ, Julian RM. Spatial and temporal genetic patchiness in marine populations and their implications for fisheries management. Calif Coop Ocean Fish Investig Rep. 1999; 94–99.
- 7. Selkoe KA, Gaines SD, Caselle JE, Warner RR. Current shifts and kin aggregation explain genetic patchiness in fish recruits. Ecology. 2006;87: 3082–3094. pmid:17249233
- 8. Johnson MS, Black R. Pattern beneath the chaos: the effect of recruitment on genetic patchiness in an intertidal limpet. Evolution. 1984; 1371–1383.
- 9. Hedgecock D. Temporal and spatial genetic structure of marine animal populations in the California Current. Calif Coop Ocean Fish Investig Rep. 1994;35: 73–81.
- 10. Kinlan BP, Gaines SD. Propagule dispersal in marine and terrestrial environments: a community perspective. Ecology. 2003;84: 2007–2020.
- 11. Buston PM, Fauvelot C, Wong MY, Planes S. Genetic relatedness in groups of the humbug damselfish Dascyllus aruanus: small, similar-sized individuals may be close kin. Mol Ecol. 2009;18: 4707–4715. pmid:19845858
- 12. Iacchei M, Ben-Horin T, Selkoe KA, Bird CE, García-Rodríguez FJ, Toonen RJ. Combined analyses of kinship and FST suggest potential drivers of chaotic genetic patchiness in high gene-flow populations. Mol Ecol. 2013;22: 3476–3494. pmid:23802550
- 13. Shapiro DY. On the possibility of kin groups in coral reef fishes. Ecol Deep Shallow Coral Reefs. 1983;1: 39–45.
- 14. Elphinstone MS, Hinten GN, Anderson MJ, Nock CJ. An inexpensive and high-throughput procedure to extract and purify total genomic DNA for population studies. Mol Ecol Notes. 2003;3: 317–320.
- 15. Hepburn RI, Mottillo EP, Bentzen P, Heath DD. Polymorphic microsatellite loci for the masked goby, Coryphopterus personatus (Gobiidae). Conserv Genet. 2005;6: 1059–1062.
- 16. Hogan JD, Galarza JA, Walter RP, Heath DD. Bayesian analysis of population structure and the characterization of nine novel microsatellite markers for the study of a Caribbean coral reef Gobiid (Coryphopterus personatus) and related taxa. Mol Ecol Resour. 2010; Available: http://eprints.eriub.org/1452/
- 17. Jombart T. adegenet: a R package for the multivariate analysis of genetic markers. Bioinformatics. 2008;24: 1403–1405. pmid:18397895
- 18. Engels B. HWxtest: Exact Tests for Hardy-Weinberg Proportions. 2014. Available: http://CRAN.R-project.org/package=HWxtest
- 19.
Dyer RJ. gstudio: Analyses and functions related to the spatial analysis of genetic marker data. 2014.
- 20. Dray S, Dufour A-B, others. The ade4 package: implementing the duality diagram for ecologists. J Stat Softw. 2007;22: 1–20.
- 21. Wang J. Triadic IBD coefficients and applications to estimating pairwise relatedness. Genet Res. 2007;89: 135–153. pmid:17894908
- 22. Pew J, Muir PH, Wang J, Frasier TR. related: an R package for analysing pairwise relatedness from codominant molecular markers. Mol Ecol Resour. 2015;15: 557–561. pmid:25186958
- 23. Queller DC, Goodnight KF. Estimating Relatedness Using Genetic Markers. Evolution. 1989;43: 258–275.
- 24. Jombart T, Devillard S, Balloux F. Discriminant analysis of principal components: a new method for the analysis of genetically structured populations. BMC Genet. 2010;11.
- 25. Jones OR, Wang J. COLONY: a program for parentage and sibship inference from multilocus genotype data. Mol Ecol Resour. 2010;10: 551–555. pmid:21565056
- 26. Wang J. Parentage and sibship exclusions: higher statistical power with more family members. Heredity. 2007;99: 205–217. pmid:17487215
- 27.
Kolaczyk ED, Csárdi G. Statistical analysis of network data with R. Springer; 2014. Available: http://link.springer.com/content/pdf/10.1007/978-1-4939-0983-4.pdf
- 28. Renyi A, Erdos P. On random graphs. Publ Math. 1959;6: 5.
- 29. Csardi G, Nepusz T. The igraph software package for complex network research. InterJournal Complex Syst. 2006;1695: 1–9.
- 30. Chapuis M- P, Estoup A. Microsatellite Null Alleles and Estimation of Population Differentiation. Mol Biol Evol. 2007;24: 621–631. pmid:17150975
- 31. Hewitt R. The value of pattern in the distribution of young fish. Rapp Proces-Verbaux Reun Cons Int Pour L’Exploration Mer. 1981;178: 229–236.
- 32. Hunter JR. The feeding behavior and ecology of marine fish larvae. Fish Behav Its Use Capture Cult Fishes. 1980; 287–330.
- 33. Bernardi G, Vagelli A. Population structure in Banggai cardinalfish, Pterapogon kauderni, a coral reef species lacking a pelagic larval phase. Mar Biol. 2004;145: 803–810.
- 34. Ronce O. How does it feel to be like a rolling stone? Ten questions about dispersal evolution. Annu Rev Ecol Evol Syst. 2007; 231–253.
- 35.
Thresher RE. Reproduction in reef fishes. Neptune City, New Jersey: T.F.H. Publications; 1984.
- 36. Medvinsky AB, Tikhonova IA, Aliev RR, Li B-L, Lin Z-S, Malchow H. Patchy environment as a factor of complex plankton dynamics. Phys Rev E. 2001;64: 021915.
- 37. Hamner WM, Jones MS, Carleton JH, Hauri IR, Williams DM. Zooplankton, Planktivorous Fish, and Water Currents on a Windward Reef Face: Great Barrier Reef, Australia. Bull Mar Sci. 1988;42: 459–479.
- 38. Frankham R. Conservation genetics. Annu Rev Genet. 1995;29: 305–327. pmid:8825477
- 39. Hamilton WD. The Genetical Evolution of Social Behaviour. II. J Theor Biol. 1964;7: 17–52. pmid:5875340