First come, first served: Possible role for priority effects in marine populations under different degrees of dispersal potential

Studying clearly delineated populations in marine lakes, islands of sea, we investigated the interplay of habitat size, dispersal potential, and priority effects in shaping marine population genetic structure.


| INTRODUC TI ON
The Coral Triangle (located between the Philippines, Indonesia and Papua New Guinea) houses the global maximum of marine biodiversity (Hoeksema, 2007;Roberts et al., 2002), yet the origins of this high diversity remain unresolved. Phylogeographical studies have shown that there are high levels of genetic structuring in populations of many marine taxa in the Coral Triangle, even at scales of less than 100 km (e.g. Kochzius & Nuryanto, 2008;Tornabene, Valdez, Erdmann, & Pezold, 2015). These observations indicate that isolating mechanisms over small spatial scales may contribute to the patterns of high global marine diversity in the Coral Triangle. However, barriers to dispersal in the sea and the timing of population divergences are elusive (Bowen et al., 2016;DeBoer et al., 2014;Peijnenburg & Goetze, 2013).
Genetic breaks in the Coral Triangle are often considered to be linked to intermittent periods of isolation caused by Pleistocene glacial cycles, which resulted in multiple smaller basins and heterogeneous environments (e.g. variation in temperature, salinity and pH) (e.g. DeBoer et al., 2014;Hoeksema, 2007).
Defining populations a posteriori based on genetic data can lead to circular reasoning because the population's identification and subsequent genetic characterization are coupled (Lowe & Allendorf, 2010;Waples & Gaggiotti, 2006). To accurately define populations and the area they reside in, marine biotopes with clear delineations are required.
Marine lakes, islands of sea, provide a setting of clearly defined marine populations residing in basins of different sizes and levels of porosity of landscape barriers, which influence dispersal potential among populations. Marine lakes are anchialine systems, small bodies of landlocked seawater, that are connected to the surrounding ocean only through underground caverns and fissures (Hamner & Hamner, 1998;Holthuis, 1973). Similar to island systems, marine lakes provide discrete and replicated microcosms of biological communities. They are connected to the sea by varying degrees, as reflected by reduced tidal fluctuation in the lakes. At one end of the spectrum are marine lakes that have high water exchange with the adjacent sea. These lakes have physical characteristics and biological communities that resemble lagoons, such as diverse reef fish communities and scleractinian corals. On the other end of the spectrum are highly isolated lakes, with limited water inflow from the sea, containing depauperate species communities. A relatively large number of marine lakes (10s-100 s) are located in Indonesia and Palau (Becking et al., 2011;Dawson, Martin, Bell, & Patris, 2009), and more are being discovered in remote areas of Indonesia (Becking, de Leeuw, & Vogler, 2015). The great majority of marine lakes are less than 50 m deep. Since seawater levels were approximately 110-140 m lower during the Last Glacial Maximum than at present (Geyh, Streif, & Kudrass, 1979;Voris, 2000), these lakes must have been dry or contained fresh water at that time (Dawson, 2006). When sea levels rose, the lakes filled up with seawater, which has been estimated to have occurred around 8,000 years before present for the great majority of lakes (Dawson, 2006;Sathiamurthy & Voris, 2006). Despite being young environments, marine lakes harbour unique, possibly endemic, genetic and species diversity (Becking et al., 2011;Dawson & Hamner, 2005;Maas et al., 2018). Selection of marine lakes with different sizes, degrees of connection to the adjacent sea and at different spatial scales offers an opportunity to assess the influence of area, geographical distance and porosity of landscape barriers in shaping marine population genetic structure. A key question is whether reductions in area of available habitat at small scales may influence genetic diversity. A larger lake could allow multiple colonization events of haplotypes or divergence of populations following colonization, resulting in a positive interaction between size of the lake and genetic diversity. A decrease in the rate of immigration due to distance or landscape barriers can lead to lower local genetic the Leiden University Fund (LUF)/ Slingelands, Singapore Airlines, the A.M. Buitendijk Fund and the J.J. ter Pelkwijk Fund (Naturalis).
Handling Editor: Richard Ladle habitat size appears to reduce genetic diversity, even at very small spatial scales. Our findings are relevant in the context of ongoing alterations to coastal hydrodynamics, which lead to habitat reduction and influence migration among populations at fine spatial scales.

K E Y W O R D S
anchialine ecosystems, Brachidontes mussels, coral triangle, isolation-by-distance, marine biodiversity, phylogeography diversity, but diminished gene flow can also allow populations to diverge and ultimately form local diversity (Hendry, Bolnick, Berner, & Peichel, 2009;Vellend, 2005). When a lake has a higher connection to the sea, that is, is closer to the source pool, the potential influx of genetically different propagules could be more frequent, resulting in a positive correlation between genetic diversity and the degree of connection. If, however, a reduction in connection causes a time lag in the arrival of immigrants, then priority effects could in theory prevent a correlation between connection and genetic diversity. Priority effects arise when first colonizers have an advantage and can shape subsequent population genetic structure (De Meester et al., 2002;Fukami, 2015;Orsini, Vanoverbeke, Swillen, Mergeay, & De Meester, 2013).
Rapid population growth and local adaptation upon colonization of a new habitat could result in the effective monopolization of resources. Following population dominance, fewer genotypes will establish than expected under ongoing immigration, promoting strong genetic differentiation even at fine spatial scales (De Meester et al., 2016). Priority effects are, however, stochastic and therefore hard to prove directly, requiring assessment by a process of elimination.
One of the most dominant invertebrates in marine lakes are mussels of the genus Brachidontes, which are found almost ubiquitously in marine lakes in Palau and Indonesia (Becking et al., 2016;Goto, Tamate, & Hanzawa, 2011;Hanzawa et al., 2012;Maas et al., 2018, this study). Brachidontes spp. are sessile organisms that disperse through broadcast spawning and have a pelagic larval stage of approximately 3 weeks in the tropics (Monteiro-Ribas, Rocha-Miranda, Romano, & Quintanilha, 2006). Phylogeographical studies of Brachidontes populations in marine lakes in Palau (Goto et al., 2011) and Indonesia (Becking et al., 2016;Maas et al., 2018) indicated that there are multiple deeply diverged lineages with significantly different shell-morphology and that these lineages constitute a species complex of at least four undescribed species. Becking et al. (2016) showed that three neighbouring lakes each harboured a different lineage and found possible evidence of local divergence. With respect to the population genetic structure of Brachidontes sp., Maas et al. (2018) examined one lineage that inhabited seven marine lakes in Indonesia, which were similar in age and size but varied in environmental regimes and degrees of connection to the sea. The authors hypothesized that even incomplete dispersal barriers may cause sufficient isolation to allow priority effects to influence the population genetic structure.
Given the possible importance of priority effects for low levels of isolation (de Meester et al., 2016;Maas et al., 2018), the current study set out to examine the interplay between dispersal potential, habitat size and priority effects in marine lake populations by significantly increasing the number of sample sites. Here, we compare mussel populations of 22 marine lakes in Indonesia and Palau and 4 coastal locations from Indonesia, Papua New Guinea and Australia.
The distance between lakes ranges from 100 to 1,400 km. The lakes are of similar age (~8,000 years), but differ in size (0.04-4.7 km 2 ) and degree of connection to the adjacent sea (from high water exchange to very little exchange). By analysing population genetics and geometric morphometrics, we aimed to (a) identify the distribution of lineages of Brachidontes mussels in the Indo-Pacific, (b) test whether marine genetic diversity conforms to the prediction of accumulation of diversity with area and (c) assess which levels of dispersal potential can lead to local divergence and structure among populations at small spatial and temporal scales. AB465572 and AB509361), respectively. In most marine lakes, mussels were found in large abundances, forming dense beds in and beneath the intertidal zone along the rocky shoreline and mangrove roots. Mussels were sampled along the perimeter of the lakes, from multiple beds. Degree of water exchange between the marine lakes and the adjacent sea was assessed by placing a water level logger inside each lake, as well as one just outside in the surrounding sea, during a 48-hr period. The measurements of these Onset® HOBO U20L loggers were converted from pressure (Pa) to depth (m) using Hoboware® Pro 3.7.16 software (Fig. S5). The 'relative tidal amplitude' is calculated by (Lake max -Lake min )/(Sea max -Sea min ), where Lake max and Lake min refer to the maximum and minimum water level in the marine lake during a 48-hr period, and Sea max and Sea min represent the maximum and minimum water level in the adjacent open sea during that period. The maximum value of 1 would indicate no or very limited obstruction to water flow in and out of the lake.

| Sample collection
Conversely, the minimum of 0 would mean that there is no water exchange with the sea.

| DNA extraction and gene amplification
DNA was extracted from the posterior adductor muscles of the mussels, using the DNeasy Blood&Tissue kit (Qiagen), following the manufacturers protocol. Using gel electrophoresis, DNA extracts were viewed under UV light on 1% agarose gel stained with ethidium bromide to check for quality. Subsequently, polymerase TA B L E 1 Sample locations with genetic diversity indices. Per locality is provided: number of individuals sequenced for cytochrome oxidase 1 (N), number of CO1 haplotypes, nucleotide diversity (π), haplotype diversity (h), mean pairwise differences and Tajima (aligned length), using the primers D23F/D6R (Park & Foighil, 2000) and 22F/1789R (Medlin, Elwood, Stickel, & Sogin, 1988

| Sequence assembly and alignment
Forward and reverse sequences were aligned and edited using CodonCode aligner 5.1.5 (CodonCode Corporation), using the Muscle algorithm (Edgar, 2004). Chromatograms of the resulting alignments were checked by eye and primer ends were trimmed.
The bivalve origin of the obtained sequences was verified through blast searches (http://blast.ncbi.nlm.nih.gov). To create a second dataset with only unique sequences, the final alignment was collapsed using DAMBE 5.5.29 (Xia, 2013). Many bivalve molluscs, in- M-type sequences resulting from DUI were filtered out using methods described in Goto et al.(2011) and Becking et al. (2016), and subsequent analyses were focused on F-type sequences only.
Genbank numbers of Mytilidae samples that were used in the analysis are shown in Figures S1-S3.

| Model selection, maximum likelihood and Bayesian inference
To classify the different lineages, phylogenetic trees were made. The Net nucleotide divergence (Nei & Li, 1979) between the main lineages in the resulting phylogenetic tree were calculated in Mega 6.06, using the model with the best fit available (Tamura & Nei, 1993)

| Morphometric analysis
To identify morphospecies, geometric morphometric analyses were processed for a total of 419 digital images of 13 mussel populations in Indonesia (lineage A-D, F). The shells of lineage E were not available for analysis. The analyses were performed following the procedure described in Becking et al. (2016). Briefly, mussels were positioned in a standardized way and photographed. Shell outlines were drawn via curves starting from the beak of the mussel using tpsDig (Rohlf, 2010). Curves were converted to semi-landmarks and standardized for size and orientation in tpsRelw (Rohlf, 2010).
We obtained centroid size and relative warp scores for each specimen. Furthermore, repeatability was confirmed for the first three relative warp axes; thus, subsequent analyses were performed only on the first three axes. To study differences among groups, we performed a nonparametric permutational multivariate analysis of variance (PERMANOVA, 1,000 permutations) implemented in past 2.11 based on Euclidean distances (Hammer, Harper, & Ryan, 2001). We corrected for multiple testing by calculating Bonferroni-corrected P values.

| Genetic diversity and demographic history
Genetic diversity was estimated using haplotype (h) (Nei, 1987) and nucleotide (π) (Tajima, 1983) diversities and mean pairwise differences in Arlequin 3.5.1.2 (Excoffier & Lischer, 2010) for each sampling location. To assess recent demographic histories, two methods were used to test the data for signatures of recent population expansions. First, Tajima's D (Tajima, 1989) was calculated, with significance determined by 10,000 random permutations using Arlequin. A significantly negative Tajima's D is consistent with recent population expansion, or a selective sweep (Tajima, 1989). Second, mismatch distributions were calculated. For each sampling location, the observed distribution of pairwise differences between sequences was compared with a theoretical distribution, as expected under a sudden expansion model (Rogers & Harpending, 1992). This theoretical distribution was computed in DnaSP (Librado & Rozas, 2009) with values for τ and θ calculated and tested in Arlequin. Harpending's raggedness index (rg; Harpending, Sherry, Rogers, & Stoneking, 1993) was determined as well. This is a measure for the smoothness of the observed mismatch distribution, which can be used to distinguish between expanded and stationary populations. The value of the raggedness index will be low and non-significant in expanding populations, while it is usually high and significant in stationary populations (Harpending, 1994;Harpending et al., 1993).

| Population structure of lineage A
To investigate the population structure within lineage A, pairwise We transformed the Φ ST values to calculate pairwise genetic distance between populations (Φ ST /(1 − Φ ST )) (Slatkin, 1995) and correlate this with geographical distance. The geographical distance matrix was calculated using lake coordinates and determining the minimum pairwise geographical distances between lakes via the R function distm implemented in the geosphere package (Hijmans, 2016). We used Mantel tests (Mantel, 1967;Slatkin, 1993) to assess the correlation between genetic and geographical distances. The test was run with 1,000 permutations, as implemented in the function mantel (vegan package). A correlation of r > 0.6 was considered strong, and significance was assigned when the P value was smaller than 0.05. to the genus Mytilus (Fig. S1). Lineages D and F (clade two) are more diverged from each other (nucleotide divergence 42%, Table S1) than the lineages within the first clade (net nucleotide divergence 9.3%-21%) and likely represent distinct species. Lineages A, B, C and E are relatively closely related, and fit in the Brachidontes s.s.

| Demographic history
The were unimodal for the majority of the lake populations, as supported by low and non-significant raggedness index and SSD values (Fig. S4).

| Population structure
Because of the high divergence between lineages, population genetic analyses were not performed on the complete dataset but only within lineage A, which was most thoroughly sampled (12 populations, n = 389). Within lineage A, a total of 51 haplotypes was found. Haplotype sharing between marine lakes is limited: two haplotypes are shared between Palau and Papua, two within Palau and one within Papua (Figure 2). Estimates of population differentiation between lakes (pairwise Φ ST ) range from 0.015 to 0.99 (average 0.64, Figure 5a, Table S2), indicating strong genetic structuring.
All the values are significant (P < 0.01), except between Palau6 and Palau7, which share their dominant haplotype ( Figure 2) and are on the same island (Mercherchar), less than 500 m apart (Figure 1).
Many estimates are close to 1, indicating that marine lake populations are not connected by gene flow. There was no significant correlation between genetic distance and geographical distance ( Figure 5b).

| D ISCUSS I ON
Marine lakes represent relatively controlled biotopes where each lake can be seen as an independent replicate of eco-evolutionary dynamics F I G U R E 4 Haplotype diversity in relation to (a) area of lake and (b) degree of connection to adjacent sea. Linear regression was used to test effects, values for R 2 and P are shown in the bottom right corner. Relative tidal amplitude (b) is used as a proxy for connectivity with higher ratios (maximum = 1) indicating a higher degree of connection to the adjacent sea over time. By comparing mussel populations in Indonesian and Palauan marine lakes at different spatial scales, and varying in area and degree of connection to the sea, we were able to assess local drivers of genetic diversity and differentiation among marine lake populations. Our study shows that marine genetic diversity conforms to the prediction of accumulation of diversity with area and our results indicate that even a small reduction in area (10s of m 2 ) can lead to decreased population genetic diversity. On a local scale, we find repeated evidence that even when there is high water exchange, with the potential influx of larvae, there is strong population structure and local divergence of populations. We conclude that a reduction in connection can cause a time lag in the arrival of immigrants, which allows first colonizers to benefit from priority effects and gain a competitive advantage. This leads to an apparently random assignment of colonizing lineages and reduced gene flow between populations. Subsequent local divergence of founder populations results in rapid genetic differentiation.
Six divergent lineages, that likely represent at least six distinct species, were detected based on the genetic and morphological characteristics of mussels found in 22 marine lakes and four coastal locations in the Indo-Pacific. The genetic distances between lineages are larger than previously reported between species in the genus Brachidontes (Lee & Foighil, 2004;Sirna Terranova, Lo Brutto, Arculeo, & Mitton, 2007), and comparable to interspecific distances in other bivalves (e.g. Shearer, Oppen, Romano, & Wörheide, 2002).
Shell shape was found to be congruent with the pattern of genetic differentiation, with each lineage having distinctly different shell outlines. Morphological variation can be related to underlying genetic variation but can also be a result of phenotypic plasticity and thus reflect differences in environmental regimes (Burridge, Goetze, Raes, Huisman, & Peijnenburg, 2015;Dawson & Hamner, 2005;Luttikhuizen, Drent, & Baker, 2003;Mariani, Peijnenburg, & Weetman, 2012). In marine lake jellyfish populations, morphological variation was found to correlate with environment rather than genetic lineage (Swift, Gómez Daglio, & Dawson, 2016).
Three similar patterns are observed in the marine lake populations: a single lineage per lake, accumulation of genetic diversity with area and demographic patterns of founder populations that have expanded. The similarity of these patterns found for different lineages and different lakes suggests general processes are at play.
First, only one genetic lineage was found in each lake, and even lakes that are less than 5 km apart could harbour highly divergent lineages. We assume that the floodwaters that filled marine lakes during the Holocene sea level rise (<8,000 years before present) allowed for independent colonization of lakes by propagules from the surrounding sea (Dawson, 2006;Dawson & Hamner, 2005).
Preliminary studies of sediment cores from the marine lakes indicate that Brachidontes was present at the onset of the lakes in Raja Ampat (Klei, 2015;Maas et al., 2018). A variety of taxa in marine lakes (jellyfish, fish, bivalves) show the same consistent pattern of a single genetic lineage per lake (Dawson & Hamner, 2005;Goto et al., 2011;Gotoh, Sekimoto, Chiba, & Hanzawa, 2009;Rose, Masonjones, & Jones, 2016). Inter-lineage competition could explain this pattern, and is expected for closely related species that occupy similar niches. Whichever lineage colonizes a lake first can have a competitive advantage from priority effects, resulting in dominance over any later arriving immigrants Fukami, 2015;. By significantly reducing successful immigration of other lineages, priority F I G U R E 5 Genetic structure of populations of Lineage A from marine lakes and relationship between genetic distance and geographical distance. (a) Principal coordinate analysis ordination of the ΦST distance matrix (based on Table S2). The amount of variation captured per axis is indicated in the axis labels. Distances between points reflect relative genetic distances between populations. Location codes provided in Table 1. (b) Pairwise genetic differentiation [(ΦST/(1−ΦST)] of mussels as a function of log 10 transformed geographical distance (m) between marine lakes. The correlation was tested with a Mantel test, values for r and P are shown in the top left corner effects are expected to substantially reduce gene flow. The initial chance arrival of propagules results in a largely stochastic outcome, since whatever lineage colonizes first can proliferate and out-compete any subsequent arrivals (Gillespie & Baldwin, 2009).
To strengthen the indirect evidence for priority effects found in the current study, future research should analyse sediment cores to assess shell morphometric changes over time, which can give further insights into colonization history (Cole, 2017;Klei, 2015).
Second, there is an accumulation of genetic diversity with increasing area. It is striking that even a reduction at scales of a few hundred m 2 in marine lakes showed a decrease in haplotype diversity, irrespective of the degree of connection to the adjacent sea.
Typically, a larger area will support a bigger population size. Given a set mutation rate, it is expected that genetic diversity will increase with population size (Hague & Routman, 2016;Kimura, 1983).
Another possible effect of habitat size is the potential for multiple niches (Warren et al., 2014), which would enhance the maintenance of locally diverged haplotypes.
Third, each lake appears to show the genetic signature of a founder event followed by rapid population expansion (Rogers & Harpending, 1992;Slatkin & Hudson, 1991). Each lake, except Papua 2, shows a single colonization event with a single dominant haplotype and multiple private haplotypes separated by one or a few mutational steps.
It is remarkable that of the 51 haplotypes of lineage A, only five are shared by two or more lakes. We therefore suggest this pattern did not result from random redistribution of haplotypes that existed in the sea during the formation of each lake, but that the majority of the within lineage diversity is due to local divergence within the marine lakes during the approximately 8,000 years of their existence. This would correspond to a genetic structure of evolution in peripatry, similar to oceanic islands and satellite lakes of ancient rift lakes (Becking, Cleary, & de Voogd, 2013, Becking et al., 2016Chen & He, 2009;Dawson & Hamner, 2005;Emerson & Gillespie, 2008;Genner et al., 2007;Gotoh et al., 2011;Goto et al., 2011;Rosindell & Phillimore, 2011 Depending on the degree of connection to the sea, there is high water flow with potential gene flow into many of the marine lakes in this study. The most isolated lakes have small fissures connecting the lakes to the surrounding sea, which likely obstruct the entrance of larvae. In contrast, the well-connected lakes have short, wide tunnels that presently allow the entry of even large organisms such as adult fish and turtles. These lakes harbour typical reef flat species communities and have similar tidal regimes compared to the adjacent sea, which would indicate the possibility of waterflow as a vector for larvae. Yet, we observed significant genetic structure even in these well-connected lakes. This structure, resulting from genetic drift and/or adaptation, would not be expected under ongoing dispersal. Remarkably, there was no relationship between genetic diversity and the degree of connection to the surrounding sea. Priority effects potentially provide an explanation for the population structure we observe (De Meester et al., 2016;Fukami, 2015). A key element of priority effects is that early-arriving colonizers are able to grow rapidly in population size, thereby substantially occupying or modifying niches before others arrive. As a result, novel migrants are not able to establish. Small passive dispersers with short generation times and high population growth, such as marine mussels, could monopolize resources and, subsequently, habitats (De Meester et al., 2016). Given the high growth rates and fecundity of Brachidontes mussels, a slight barrier to dispersal could be enough for a dominant role of priority effects in shaping populations (Cilia & Deidun, 2012;Fukami, 2015;Morton, 1988).
The scale at which we see structure and divergence is at the level On a final note, marine lakes are ecologically valuable systems that are currently not considered in the conservation zoning of the MPAs in Raja Ampat (West Papua, Indonesia) (Maas et al., 2020).
Marine lakes are, however, threatened marine ecosystems due to aquaculture and unregulated tourism (Becking et al., 2011;Mangubhai et al., 2012). Some of the lakes from the current study are filled with jellyfish that are gaining the attention of tourists as the region becomes more accessible (Becking et al., 2015;Maas et al., 2020). In Palau, increased tourism pressure has been associated with the collapse of jellyfish populations and the introduction of invasive species (Dawson, Martin, & Penland, 2001). The data of the current study indicate that marine lakes in Raja Ampat contain unique genetic diversity and isolated populations.

PERMITS
We are grateful to the Indonesian Institute of Sciences (LIPI) and the

ACK N OWLED G EM ENTS
We would like to thank R. Moolenbeek and F. Wesselingh for identification of the samples and suggestions for morphometric characters, Alice Burridge for suggestions with morphological analyses