Reproductive skew in Japanese sardine inferred from DNA sequences

ICES Journal of Marine Science ICES Journal of Marine Science (2016), 73(9), 2181– 2189. doi:10.1093/icesjms/fsw070 Editor ’s Choice Reproductive skew in Japanese sardine inferred from DNA sequences Hiro-Sato Niwa*, Kazuya Nashida, and Takashi Yanagimoto *Corresponding author: tel: +81 45 788 7691; fax: +81 45 788 5001; e-mail: Niwa, H-S., Nashida, K., and Yanagimoto, T. Reproductive skew in Japanese sardine inferred from DNA sequences. – ICES Journal of Marine Science, 73: 2181–2189. Received 6 November 2015; accepted 5 April 2016; advance access publication 18 May 2016. An excess of low-frequency mutations is a ubiquitous characteristic of many marine species, and may be explained by three hypotheses. First, the demographic expansion hypothesis postulates that many species experienced a post-glacial expansion following a Pleistocene population bottleneck. The second invokes some form of natural selection, such as directional selection and selective sweeps. The third explanation, the reproductive skew hypothesis, postulates that high variation in individual reproductive success in many marine species influences genetic diversity. In this study, we focused on demography and reproductive success and the use of coalescent theory to analyse mitochondrial DNA sequences from the Japanese sardine. Our results show that population parameters estimated from both the site-frequency spectrum and the mismatch distribution of pairwise nucleotide differences refute the demographic expansion hypothesis. Further, the observed mismatch distribution, compared with the expectations of the reproductive skew hypothesis, supports the presence of multiple mergers in the genealogy. Many short external branches but few long terminal branches are found in the sardine genealogy. Model misspecification can lead to misleading contemporary and historical estimates of the genetically effective population sizes in marine species. The prevalence of reproductive skew in marine species influences not only the analysis of genetic data but also has ecological implications for understanding variation in reproductive and recruitment patterns in exploited species. Keywords: coalescent, effective population size, excess singletons, multiple mergers, recruitment variation, skewed offspring-number distribution. Introduction A growing body of molecular data show that many marine species with broadcast spawning have star-like gene genealogies. Many studies have concluded from the analysis of mitochondrial DNA (mtDNA) sequences with coalescence-based methods, such as Bayesian skyline plots (BSPs), or the analysis of nucleotide differences between sequences, that populations of many marine species experienced a recent rapid demographic expansion (e.g. Crandall et al., 2012). Under this interpretation, populations were depressed by the environmental effects of the Last Glacial Maximum, and post-glacial warming promoted population expansions. Looking back in time, a bottleneck generates a brief period of rapid coalescences in a gene genealogy just before a population expansion, and mutations appearing after the expansion occur on external branches of the genealogy and appear as singleton haplotypes in a sample. As a consequence, sequence mismatch distributions for these populations are unimodal (Slatkin and Hudson, 1991; Rogers and Harpending, 1992). It is noteworthy that the putative onset of population expansion in many marine species predated the last glacial maximum 18 –20K years ago (Janko et al., 2007; Grant, 2015; Sromek et al., 2015), suggesting that populations were unaffected by subsequent glacial cycles. The multiple-merger model based on reproductive skew among individuals (e.g. Eldon, 2011) is an alternative to the population expansion model and can account for the excesses of singleton mutations (Durrett and Schweinsberg, 2005; Berestycki et al., 2007). The reproductive skew model is based on the life history trait common to most marine fish and invertebrates with type-III survivorship curves, which spawn large numbers of eggs producing larvae with high mortality rates early in life. The heavy-tailed, power-law distribution has been used to model the distribution of offspring number among individuals (Schweinsberg, 2003). When highly successful reproduction events occur, a large fraction of the population arises from only a few individuals who win ‘reproduction sweepstakes’ (Beckenbach, 1994; Hedgecock, 1994; Hedgecock and Pudovkin, 2011). Marine fish and invertebrate populations are # International Council for the Exploration of the Sea 2016. All rights reserved. For Permissions, please email: National Research Institute of Fisheries Science, Yokohama 236-8648, Japan 2182 We used the complete sequence of Sardinops melanostictus mtDNA (NC002616) to design the SMThrL primer, 5′ -tgccccagtagctta gttcaa-3′ (forward), and the SMPheH primer, 3′ -gcttcttacggccca tctta-5′ (reverse) to amplify the mtDNA CR. Each PCR was performed in 10 mm3 containing 5 –10 ng template DNA, 1 mm3 of 10× reaction buffer, 1 mol m23 each dNTPs, 2 × 1023 mol m23 each primer, and 0.5 units of EX Taq HS polymerase (Takara, Shiga, Japan) in an ABI 9700 Thermal cycler (Applied Biosystems, Foster City, CA, USA). Initial denaturation was for 2 min at 948C; followed by 40 cycles of 30 s at 948C for denaturation, 30 s at 558C for annealing, and 30 s at 728C for extension; and a final extension at 728C for 10 min. PCR product was purified with GFX (Qiagen, Hilden, Germany) and was used as the template DNA for cycle sequencing reactions performed using Big Dye Terminator Cycle Sequencing Kit (Ver.3.1, Applied Biosystems). Sequencing was conducted on an ABI Prism 3730XL (Applied Biosystems) automatic sequencer with four primers containing PCR primers and inner primers (SMCRL: 5′ -gtagtaagaaccgaccaaccg-3′ and PheHR primer: 5′ -tcaaagcataacactgaagatgttaa-3′ ). The mtDNA CR was 1212 bp long for 106 sardine. These sequences were deposited in GenBank (accession nos. LC031518–LC031623). Methods and results Preliminary exploration of mtDNA sequence variation We identified 101 haplotypes among the 106 sardine sequences. Most haplotypes were in low frequency: 98 were observed only once, 2 were observed twice, and 1 was observed four times. Figure 1 shows the folded site-frequency spectrum (SFS) for the sardine mtDNA CR sequences. Compared with the expected SFS in a neutral equilibrium population assuming the Kingman coalescent under an infinitely-many-sites model (ISM) of mutation (Watterson, 1975; Tajima, 1989), we observed an excess of singleton sites. An abundance of low-frequency mutations likely resulted from a star-like topology with a short genealogy, whereas middlefrequency polymorphisms would be a signature of long population history. An unrooted maximum likelihood tree was constructed from the mtDNA CR sequences using program phyml (Guindon et al., 2010). The tree branched unevenly (Figure 2). There was a burst of mergers at one (...truncated)