Population expansions dominate demographic histories of endemic and widespread Pacific reef fishes

Despite the unique nature of endemic species, their origin and population history remain poorly studied. We investigated the population history of 28 coral reef fish species, close related, from the Gambier and Marquesas Islands, from five families, with range size varying from widespread to small-range endemic. We analyzed both mitochondrial and nuclear sequence data using neutrality test and Bayesian analysis (EBSP and ABC). We found evidence for demographic expansions for most species (24 of 28), irrespective of range size, reproduction strategy or archipelago. The timing of the expansions varied greatly among species, from 8,000 to 2,000,000 years ago. The typical hypothesis for reef fish that links population expansions to the Last Glacial Maximum fit for 14 of the 24 demographic expansions. We propose two evolutionary processes that could lead to expansions older than the LGM: (a) we are retrieving the signature of an old colonization process for widespread, large-range endemic and paleoendemic species or (b) speciation; the expansion reflects the birth of the species for neoendemic species. We show for the first time that the demographic histories of endemic and widespread reef fish are not distinctly different and suggest that a number of processes drive endemism.

The existence of endemic coral reef fish species is a challenge for evolutionary biologists to explain given tropical waters are widely connected. The highest levels of endemism on coral reefs are observed near islands at the peripheries of the Indo-Malay-Philippines Archipelago (IMPA), global hotspot of reef fish species diversity [1][2][3][4][5] . As examples from the Pacific Ocean, the percentage of reef fish that are endemic in the Hawaiian archipelago is 25% 6 , 22% in Easter Island 7 , and 14% in the Marquesas Islands 8 . Even if actual patterns of distribution of coral reef fish species are now well depicted, the evolution and processes underlying the establishment and maintenance of endemic species remains unclear. Despite their unique nature and their potentially higher risk of extinction, the origin and population history of endemic species is poorly studied.
The Pleistocene era (c. 1.8-0.01 Ma), which affected the distribution and demographic history of both terrestrial and costal marine species 9 , was characterized by glacial cycles and sea level fluctuations up to 150 m below present sea level 10,11 . Large parts of continental shelves were exposed during low sea level, altering shallow water habitat and likely reducing coral reef area [12][13][14] . These dramatic changes to the environment influence the demographic history of populations, leaving a footprint in the pattern of genetic diversity 15 that will vary depending on the range extent of the species (i.e., how much of the global species was affected by the event). Population bottlenecks and expansions have often been retrieved in marine populations and shown to coincide with the last major sea level changes impacting population dynamics 9,[16][17][18] .
Species biology can affect how reef fish species respond in term of population size variation to major sea level changes such as through varying the larval phase or habitat [19][20][21][22] . However, very few studies have explored how species with different range size respond to major climatic events 23,24 . Endemic species have by definition a limited distribution 25 so are widely expected to be highly vulnerable to environmental changes that are locally disrupting 26,27 . This is particularly the case for coral reef fishes forming metapopulations where connectivity occurs only during the larval stage, i.e. colonization of new suitable habitats and migrant exchanges among already established populations. Populations of widespread species even at the edges of their geographic range can still be considered as part of a large metapopulation. The exchange of migrants through larval connectivity, even if very infrequent, can enable recovery from major environmental changes 28 . In contrast, endemic species with a limited geographical distribution cannot rely on outcrossing with source populations, being either a newly established species (neoendemism) or the remnant of an ancestral widespread species (paleoendemism). In summary, two main differences characterize endemic vs widespread species: the size of the habitat and the degree of connectivity. We therefore expect to see major differences in the demographic histories of endemic and widespread species.
We investigate the genetic diversity and demographic history of 28 reef fish species from the Gambier and Marquesas archipelagos, representing five major reef fish families and selected to include both endemic and widespread species of the same genus. In particular, we examine whether range size is a determinant of the demographic history we retrieve; i.e., do demographic histories vary with range size? Examining multiple pairs of close related species allows us to infer if co-distributed species shared the same demographic history, and test the potential influence of their range distribution. We used both mitochondrial and nuclear DNA sequences to examine demographic history. The use of several independent markers provides a replicate of the coalescent process and improves demographic estimates by reducing the coalescent variance [29][30][31] . Historically, mitochondrial markers have been widely used in population genetics because of the availability of universal primers (at least for groups of related species) and because of the lack of intralocus recombination, which may bias demographic inferences. However, nuclear markers can now be more easily typed and can also be used to examine demographic history. Using both marker types can enhance our ability to depict complex evolutionary history of species 32 . This study represents the first large-scale comparison of the demographic history of endemic versus widespread marine species. We resolve the demographic histories of the subject species and discuss the potential underlying evolutionary processes that led to present day patterns.
Past population size changes. Overall, we retrieved population expansion for almost all species, regardless of their range size or locality. Indeed, we found population expansions for 24 species and a constant effective size for 4 species (Tables 2 and 3).
Neutrality tests and Extended Bayesian Skyline Plot (EBSP) analysis allowed us to retrieve population expansion for 24 species. For 19 species, we found a population expansion with significant neutrality test F S for at least one gene for both type of marker, EBSP rejecting a constant size model in favor of one demographic change and a clear graphical expansion (Table 2, Figs 1 and S1). For 5 species, we retrieved a population expansion only from the mitochondrial markers while nuclear markers indicated a signal of constant population. Approximate Bayesian Computation (ABC) analyses were concordant with negative growth rates retrieved, indicating population expansions ( Table 2). The species for which expansions were found (19 of 28 with two types of marker, 24 of 28 with one type of marker) are from both archipelagos and include both endemic (5 large-range, 8 small -range) and widespread species (8 species), 4 of the 5 families (except Chaetodontidae), and all types of reproduction (Table 2).
We find a constant effective population size only for 4 species (Tables 2 and 3). For these, neutrality tests and EBSP analysis indicated a constant population size for either both type of markers (Plectroglyphidodon lacrymatus) or most of the analyses (Chaetodon declivis, C. citrinellus and Plectroglyphidodon leucozonus) (Tables 2 and 3). All species are from Marquesas; three are widespread species and one is a large-range endemic species. They belong to two different families (Chaetododontidae and Pomacentridae) that differ in their reproduction strategies (pelagic eggs released vs. benthic spawners species).
Timing of the demographic events. We found overall mostly signal of population expansion (24/28 species), with expansion times that varied widely among species, ranging from 8,000 YBP to 2,000,000 YBP, no matter the method used or the genes screened.
For mitochondrial markers, EBSP analysis shows that expansion time varied from 15,000 YBP to 440,000 YBP (Table 3, Fig. 1 and Figs S1 and S2). The timing of expansions retrieved does not differ significantly among species. None of the tested factors were significant: family (H 3 = 2.8578, p = 0.41), archipelago (W = 55.5, p = 0.62), range-size classification (W = 57.5; p = 0.57 and H 2 = 0.5568, p = 0.76), and reproductive strategy (H 2 = 2.6737, p = 0.26) ( Fig. 2 and Fig. S3). For nuclear markers, EBSP analysis shows that expansion times varied more widely among species than was found for mitochondrial markers, ranging from 8,000 YBP to 2,000,000 YBP (Table 3, Figs S1 and S2). As it was the case with the mitochondrial markers, the timing of expansions retrieved for nuclear markers does not differ significantly among species. Again, none of the tested factors were significant: family We compared the population expansion timing retrieved with the EBSP (T EBSP ) for the two types of marker. We found a concordant time (within 25,000 years) of expansion for some species (e.g. Epinephelus irroratus, Ostorhinchus relativus, or Stegastes aureus). For others, T EBSP varied greatly between marker types within the same species (e.g. Acanthurus nigricans, Pristiapogon kallopterus or Chromis fatuhivae) ( Table 3, Fig. 3 and Fig. S2). Considering the 24 species displaying population expansion, we found T EBSP that varied by a factor of 26 for the mitochondrial marker and 250 for the nuclear marker. When narrowed to family, T EBSP varied for Pomacentridae by a factor from 1.7 (Acanthuridae) to 26.7 (Pomacentridae) for the mitochondrial markers and from 2.27 (Apogonidae) to 250 (Pomacentridae) for the nuclear markers. Finally, comparisons between closest relatives (i.e. same genus) revealed both very similar T EBSP and very divergent ones, with a variation by a factor of 1.04 (Pseudogramma) to 13.33 (Chromis) for mitochondrial markers and 1.6 (Epinephelus) to 12.14 (Acanthurus) for nuclear markers.
Expansion times and credible intervals provided by the ABC analyses (T ABC + CI) relaxed the differences observed for each species considering solely the EBSP expansion time. We retrieved distribution of expansion time (i.e. T EBSP with EBSP and T ABC + CI with the ABC) concordant in general but not overlapping for 7 species, mostly due to values retrieved with the EBSP for nuclear marker.

Discussion
Population expansions dominate the demographic history of endemic and widespread Pacific reef fish. The timing of the population expansions we retrieved varies greatly, between 15,000-440,000 YBP for mitochondrial markers and between 8,000-2,000,000 YBP for nuclear markers (Table 3 and Fig. 3). There are no consistent patterns among the predictor variables with respect to when the expansions likely occurred. Population expansion timing even varied both among and within the endemic and widespread species, which was unexpected. Such expansions have usually been associated with Pleistocene interglacial periods, when sea level variations likely profoundly affected habitat distribution 10,33 . Indeed, correlations between population expansion and such  39,40 , and also for other marine organisms such as gastropods 41,42 . Nine glacial cycles have been recorded in the last 800,000 years, all of which resulted in major sea level variations 43,44 that could have produced departures from genetic equilibrium conditions. However, if Pleistocene climate changes drive species demography, we should observe both expansions and population contractions following climatic oscillation. Strikingly, we found two clear patterns only in our survey; (i) population expansion; and (ii) a constant effective size, with the former largely dominating the demographic history of endemic and widespread Pacific reef fish (Table 2). Neither neutrality tests nor EBSP provided any evidences of population contraction for any markers in any species. The statistical power to detect a bottleneck is dependent on the number of markers used 45 . Unfortunately, reconstructing complex demography can be a very difficult task even when analyzing whole genome data (see for example Boitard et al. 46 ). Using relatively few markers like in this study, it is possible that no coalescence survived the bottleneck preceding the last expansion, erasing all traces of more ancient expansions 47 .
We expected then to recover only the most recent expansion, occurring after the Last Glacial Maximum (LGM), i.e., between 26,500 and 19,000-20,000 years ago 48 . Therefore, we considered that a species expanded because of the climatic change after the LGM if the estimated expansion time was < 26,500 years B.P. When looking at T EBSP we found only 5 species for which at least one of the two marker types (i.e., mitochondrial or nuclear) was strictly consistent with an expansion after the LGM (Abudefduf conformis, Chromis flavapicis, Epinephelus irroratus, Pristiapogon kallopterus, and Stegastes emeryi; Table 3 and Fig. 3). T EBSP is a point estimates inferred visually and a posteriori after the likelihood based analysis performed with BEAST. We therefore performed an ABC  are computed for each type of marker (mitochondrial and nuclear). Significance (ie. rejection or not of a constant population size model) is reported with an asterisk (*) and the trend of the curve of Ne through time is then reported (Exp: expansion; cst: constant). The mode and the credible interval (9(% low-95% up) of the growth rate retrieved with the ABC analysis are reported for each type of marker. Superscripts next to species names refer to range size classifications (W for widespread, L for large-range endemic, and S for small-range endemic).
estimate of the expansion time to compute its 95% credible interval. We found the distribution of T ABC consistent with the LGM for 9 more species with at least one of the two marker types (Abudefduf sordidus, Chromis agilis, Chrysiptera galba, C. glauca, Dascyllus strasburgi, Epinephelus fasciatus, Ostorhinchus relativus, Stegastes aureus, S. fasciolatus; Table 3 and Fig. 3). These 14 species (5 detected with T EBSP and 9 with T ABC ) probably experienced variations of their population size during the Pleistocene interglacial periods and started expanding soon after their habitats were restored.  However, we retrieved expansion times for both endemic and widespread species older than the LGM for almost half of the species. These results are totally in agreement with previous phylogeographic studies based on mitochondrial markers. Population expansions have been found for reef fish in Marquesas Islands of 15,400-23,900 YBP for Chaetodon ornatissimus 49 ; 180,000-370,000 YBP for Lutjanus kasmira and 120,000-240,000 YBP for L. fulvus 38 . We propose two other scenarios as potential explanations of the large range we find in the timing of population expansions and the lack of differences between endemic and widespread species.
(a) Signals of population expansion observed before the LGM correspond to the colonization of the archipelagos by the species through connectivity. This scenario seems to be more dominant as most expansions were older than the LGM (14/24 expansions retrieved considering expansion and mixed categories (  54 and so does C. abrupta and C. margaritifer (many C. abupta COI haplotypes matches (100% means identity) C. margaritifer COI haplotypes (GenBank numbers: FJ583159; FJ583158; and FJ583162)) while P. sagmarius present little genetic divergence with its sibling species, less than 2% (P. imparipennis COI sequence number JQ350225 (GenBank)), equivalent to a divergence date of less than 2 Mya. The expansion retrieved for these species would then correspond to their birth as endemic species. All our results are based on the hypothesis that we correctly specified the mutation rate for all our species and all the genes here considered. Rates may vary greatly from one species to another and the estimates obtained rely heavily on the statistical methods used to calibrate the clock 56,57 . Nonetheless, molecular clock rates for COI and cytochrome b for reef fish are well known and we used the 1-2% divergence rate generally used as a consensus in  the literature 23,49,56,58,59 . Concerning nuclear markers, we used a divergence rate for S7 calibrated on reef fishes 60 and polar fishes 20 . They constitute the only calibration so far obtained for this intron while no calibration has been proposed yet for the GnRH intron. However, evolution rates for nuclear DNA are usually slower than those for mitochondrial DNA 57,61 and such rate seems to be in accordance with rates of mitochondrial markers we used, though those rates can vary greatly in vertebrates 56,62,63 . We acknowledge that variability in the molecular rate of evolution among our species may exist and inflate the differences we have observed in the estimated expansion times. Nevertheless, if we consider T EBSP , we note that it varies between the youngest and the oldest estimate by a factor of 26 for the mitochondrial markers and 250 for the nuclear markers. Such wide distribution was confirmed when looking at closely related species: expansion times varied for example for Pomacentridae by a factor of 19 for the mitochondrial markers and 200 for the nuclear markers and for Chromis genus by a factor of 13,3 for the mitochondrial markers and 5,7 for the nuclear markers. Similar values are found when considering the mode of T ABC . In summary, the width of the distribution of those estimated expansion time (both T EBSP and T ABC ) for both mitochondrial and nuclear markers is so large that it cannot be explained by the uncertainty of the molecular rate or the variance associated with the estimation process. In summary, we cannot reconcile all the dates to a single environmental change. It is much more likely that each species (or at least several groups of species) was affected differently by successive specific events.
The use of two types of molecular markers allowed us to reconstruct the demographic history of reef fishes from two remote archipelagos hosting high level of endemism. We find that the demographic history of most Pacific reef fish species, despite their variety of range size and life traits, is dominated by population expansions. Consistently with the absence of difference in the demographic signature, we find that the genetic diversity of endemic and widespread species is similar. This matches the conclusions of our previous study, which was only based on one mitochondrial marker (cytochrome b) 64 , but see also Eble, et al. 23 and Hobbs, et al. 65 . We strengthen this result here by adding another mitochondrial locus and two unlinked nuclear loci. Genetic diversity and signature of effective population size change are therefore not correlated with the range-size distribution of reef fish species. More loci will be needed to refine the demographic trajectories of each species, both to describe events that could have been missed due to a lack of statistical power and to improve the parameter estimation process. Overall, this study highlights that demographic histories of endemic species do not differ from widespread species; both have complex and similar histories.

Specimen collection. The Gambier archipelago (centered on 23°S, 134°W) is southeast of the Tuamotu
Archipelago; it includes 11 high islands spreading only over 40 km and enclosed by a wide barrier reef. The Marquesas archipelago spread over 500 km between 8°-11°S and 141°-138°W and includes 12 high volcanic islands surrounded by fringing reefs 66 ; these are the northeastern-most islands of French Polynesia. We sampled a total of 1,244 reef fishes using polespears or anesthetic while SCUBA diving and exploring all islands of the Gambier Islands and of the Marquesas Islands respectively in 2010 and 2011. This dataset is composed of 28 species (~45 individuals per species), with 20 and 8 co-distributed species respectively in Marquesas and in Gambier Islands (Fig. 4a), and includes a wide range of life history traits (Table 1). Species range distribution (Fig. 4b) varied from widespread (range size > 12 000 km, 13 species), to large-range endemic (1000-8000 km, 6 species) and small-range endemic (less than 500 km, 9 species, see justification in Delrieu-Trottin et al. 64 ). Five different families are represented: Acanthuridae, Apogonidae, Chaetodontidae, Pomacentridae and Serranidae. Several reproductive strategies are represented meaning larval dispersion potential varies greatly ( Table 1). Difference of sampling between the archipelagos relates to the uniqueness of Marquesas Islands. Third highest region of endemism for coral reef fishes in the Indo-Pacific 8 , they host many more cases of endemism.
Ethics Statement. The study protocol was approved by the National Center for Scientific Research and is in accordance to the laws of the French Republic and of the collectivity of French Polynesia.

Laboratory procedures and genetic analyses.
We extracted whole genomic DNA from fin tissues preserved in 96% EtOH at ambient temperature using QIAxtractor (QIAGEN, Crawley) according to manufacturer's protocols. Fragments of the cytochrome b for 1,237 individuals (Cyt b, 739-999 bp) and cytochrome C oxydase subunit 1 for 1,005 individuals (COI, 571-688 bp), both parts of mitochondrial genome, were amplified using PCR. These were sequenced using the universal primers GLUDGL-CB3H of Palumbi et al. 67 and different combinations of primers from Ward et al. 68 . In addition, part of the third intron in the Gonadotropin-Releasing Hormone gene for 912 individuals (GnRH, 267-417 bp) and part of the first intron of the S7 ribosomal protein gene for 973 individuals (S7, 485-783 bp) were amplified. For these, we used primers GnRH3.3F-GnRH3.3R 69 and primers S7RPEX1F-S7RPEX2R 70 . Fragments were amplified using PCR protocols and sequencing as described by Williams et al. 51 . Differences in the number of samples are due to our inability to amplify cytochrome b for 3 species (Epinephelus irroratus, Epinephelus fasciatus and Ostorhinchus relativus), GnRH intron for 7 species (Dascyllus strasburgi, Chromis flavapicis, Chromis abrupta, Apogon lativittatus, Ostorhinchus apogonoides, Ostorhinchus relativus and Pristiapogon kallopterus) and S7 intron for 3 species (Apogon lativittatus, Epinephelus fasciatus and Ostorhinchus apogonoides). Overall, we worked with a dataset comprised of at least one mitochondrial marker and one nuclear marker for 26 species. Four markers were available for 19 species, 3 markers for 5 species and 2 markers for 4 species. Only mitochondrial data were available for 2 species (Table 1).
Sequences were edited using GENEIOUS PRO v.6.1.7 (Biomatters) and aligned with Clustal W 71 as implemented in GENEIOUS. Alignments were unambiguous with no indels or frameshift mutations. For the nuclear markers, allelic state from sequences with multiple heterozygous single nucleotide polymorphisms (SNPs) was estimated using PHASE v.2.1.1 as implemented in DnaSP 72-75 . We tested for recombination using SBP and GARD Scientific RepoRts | 7:40519 | DOI: 10.1038/srep40519 methods implemented online on the Datamonkey webserver 76,77 . The GARD and SBP tests failed to find recombination for any of the nuclear markers of either species.
Haplotype diversity (h) and nucleotide diversity (π ) were calculated using the software package DnaSP v.5.1 75 . To detect departures from a neutral Wright-Fisher model, we used Fu's F S neutrality test 78 . Assuming selective neutrality, significant negative values of F S indicate population growth while significant positive values are a signature of either genetic subdivision or population contraction. The test was implemented in DnaSP v.5.1 77 and its significance was determined from 1000 coalescent simulations.
Changes in effective population size (Ne) through time were estimated using the Extended Bayesian Skyline Plot (EBSP) implemented in BEAST v.1.7 79,80 . EBSP is a non-parametric model that does not specify any prior hypothesis on the tempo and mode of changes of the Ne. EBSP allows the analysis of multiple loci, reducing the stochastic variance of the coalescent process and improving the reliability of demographic inferences 47 . Mitochondrial and nuclear markers have different mechanisms of evolution and inheritance, which are also affected by reproductive strategies and sex-biased processes. The species investigated use a range of reproductive strategies (e.g. sequential hermaphrodites vs gonochoristic species, couple vs harem) so we analyze mtDNA and nuclear markers separately in the EBSP analyses.
The two nuclear loci are unlinked so have different genealogies, substitution models and clock rates but the same underlying demography. As the mitochondrial genes are in linkage by structure, we used different substitution models and clock rates but the same genealogy. Sequence divergence estimates for the two mitochondrial markers (COI and Cyt b) in reef fish range from 1% to 2% per million years 56,58,59,81 . We set a strict clock with a uniform prior distribution for the clock rate with an upper and lower mutation rate ranging from 0.5 × 10 −8 to 1 × 10 −8 per site per year. We retrieved various sequence divergence estimates for S7 intron from 0.28% to 1.7% 20,56,60 while no divergence rate was available in the literature for the third intron in the Gonadotropin-Releasing Hormone (GnRH). For these two introns, we set a strict clock with a uniform prior distribution for the clock rate with an upper and lower mutation rate ranging from 0.14 × 10 −8 to 0.85 × 10 −8 per site per year. To take into account possible site-specific variation in the mutation rate we used the HKY + G model of mutation for all genes. We ran 10 million MCMC iterations with a thinning interval of 1,000. We checked convergence by visually inspecting the trace and computing the effective sample size (which was always higher than 100) for each parameter in two independent runs using the program TRACER v.1.5 (http://tree.bio.ed.ac.uk/software/ tracer/). EBSP were plotted using R v.3.2.3 82 (R Development Core Team 2015) setting the burn-in to 10%.
Generation time is unknown for most species selected for this study and it potentially varies greatly among them as we selected a large range of families. To remove that uncertainty, we therefore set the mutation rate in units per site per year to obtain time points expressed in calendar years; independent of the generation time 80 . That method allows comparing potential expansion times between close related species but also different families. The effective size being scaled by the generation time (unknown for most species), only the timing of expansions between species will be discussed in the present work.
Finally, we compute the posterior distribution of the number of demographic changes occurring along the gene genealogies. This method formally tests how many times Ne changed through the history of the gene genealogy, allowing to reject a constant size population model. We rejected a constant population size model if the lower bound of the 95% high posterior density (HPD) of this distribution was higher than zero. Expansion times were identified from the skyline output as the point when the Ne started to increase.
The data did not meet assumptions of normality so we computed the non-parametric tests Wilcoxon-Mann-Whitney and Kruskal-Wallis to examine differences in genetic diversity (haplotype and nucleotide) and demography (TMRCA, time of expansion) among four factors: range size (endemic vs. widespread and all three range types), archipelagos, families, and reproductive strategies. We also computed Wilcoxon signed-rank test to compare the genetic diversity indices for endemic species that had a widespread congener (i.e. same genus) caught in the same archipelago. All statistical analyses were performed on R, using the vegan package 83 , and ggplot2 for graphics 84 .
To compute the full posterior distribution of the expansion time (T) for each species we developed an Approximate Bayesian Computation (ABC) approach 85 . ABC summarizes the probability of the data through a vector of observed summary statistics and it is therefore very flexible as it can investigate any demographic model provided it can be simulated 86 . EBSP showed that all species were either expanding (following approximately an exponential growth) or constant (see results). EBSP uses a likelihood function to compute the probability of the observed data and it is therefore to be preferred over methods which approximate the likelihood. However, the EBSP does not compute explicitly the T exp , which is inferred visually a posteriori and has thus no credible interval associated. We therefore implemented the following three parameters model with ABC: an ancestral constant population characterized by an effective population size (N anc ) starts to increase (or decrease) exponentially at T exp to reach at time 0 the modern effective population size (N mod ). We assigned the following uniform priors to the three parameters: N anc and N mod : {1,000 ÷ 10,000,000}; T exp : {100 ÷ 2,000,000}. Mutation rate was set in years as in the EBSP; we used a value of 10 −8 per site per year for the mtDNA genes and 8 * 10 −9 for the nuclear genes. We performed 500,000 coalescent simulations with parameter values extracted from prior distributions using fastsimcoal2 v.2.5.1 87 . We computed as summary statistics the number of haplotype, the number of segregating sites, the mean pairwise difference and Fu's Fs using alrsumstat 88 . Parameters were estimated from the 5,000 simulations closest to the observed dataset using a local linear regression according to Beaumont et al. 85 as implemented under the R environment in the library abc 89 . Analyses were performed separately for the whole mtDNA (one genealogy) and the two nuclear genes (two independent genealogies with the same underlying demography) as in the EBSP.