Timing strains of the marine insect Clunio marinus diverged and persist with gene flow

Genetic divergence of populations in the presence of gene flow is a central theme in speciation research. Theory predicts that divergence can happen with full range overlap – in sympatry – driven by ecological factors, but there are few empirical examples of how ecologically divergent selection can overcome gene flow and lead to reproductive isolation. In the marine midge Clunio marinus (Diptera: Chironomidae) reproduction is ecologically restricted to the time of the lowest tides, which is ensured through accurate control of development and adult emergence by circalunar and circadian clocks. As tidal regimes differ along the coastline, locally adapted timing strains of C. marinus are found in different sites across Europe. At the same time, ecologically suitable low tides occur at both full and new moon and twice a day, providing C. marinus with four nonoverlapping temporal niches at every geographic location. Along the coast of Brittany, which is characterized by a steep gradient in timing of the tides, we found an unusually large number of differentially adapted timing strains, and the first known instances of sympatric C. marinus strains occupying divergent temporal niches. Analysis of mitochondrial genotypes suggests that these timing strains originated from a single recent colonization event. Nuclear genotypes show strong gene flow, sympatric timing strains being the least differentiated. Even when sympatric strains exist in nonoverlapping temporal niches, timing adaptations do not result in genome‐wide genetic divergence, suggesting timing adaptations are maintained by permanent ecological selection. This constitutes a model case for incipient ecological divergence with gene flow.

The concept of population divergence with gene flow intersects with studies of local adaptation (Kawecki & Ebert, 2004;Savolainen et al., 2013), as divergence with gene flow is generally assumed to start with divergent natural selection Seehausen et al., 2014;Smadja & Butlin, 2011). For this process to result in speciation, reproductive isolation must evolve. However, it is still largely unclear how adaptive ecological divergence is maintained against gene flow at early stages of the process, when no other isolating factors exist. In principle, occupying different ecological niches can by itself reduce gene flow between diverging populations. As all abiotic or biotic factors can be part of a species' ecological niche, countless factors can potentially drive divergent ecological adaptation and many have been considered to play a role in population divergence with gene flow, e.g., soil parameters (Savolainen et al., 2006) or biotic interactions with host plants (Dres & Mallet, 2002).
An important dimension of the ecological niche is an organism's timing of activity, life history and reproduction, its temporal niche (see e.g. Häfker & Tessmar-Raible, 2020;Hut et al., 2012).
Allochrony, i.e., differences in timing between individuals, populations or species, can lead to isolation by time (IBT, Hendry & Day, 2005) and has been considered a major factor to facilitate sympatric speciation (Taylor & Friesen, 2017). Allochrony in reproductive timing automatically entails part of the reproductive isolation needed for speciation. Reproductive timing can thus be a magic trait, i.e., a trait for which ecological divergence and reproductive isolation are directly linked (Gavrilets, 2004). In the context of ecological divergence of populations with gene flow, most studies consider allochrony due to seasonal timing (Helm & Womack, 2018;Ragland et al., 2017;Tauber & Tauber, 1977;Taylor & Friesen, 2017) or daily timing (Devries et al., 2008;Fukami et al., 2003;Hänniger et al., 2017). Here, we demonstrate ecological divergence with gene flow based on daily and lunar-phase timing in the marine midge Clunio marinus. We add lunar phase as a time scale to the concept of the temporal niche and underline the importance of temporal niches in the process of population divergence with gene flow.
Clunio marinus (Diptera: Chironomidae) is primarily known as a model organism for studying circalunar clocks (Kaiser et al., 2016;Neumann, 2014), i.e., endogenous biological timekeeping mechanisms that limit reproduction and behaviour of many marine organisms to distinct lunar phases. While the existence of circalunar clocks has been repeatedly demonstrated (Neumann, 2014), their molecular basis is unknown. In C. marinus the circalunar clock has a clear ecological relevance (Kaiser, 2014) While at a given geographic location specific tidal conditions always recur at the same time during the lunar month and daily cycle, the timing of these conditions gradually shifts along the coastline. C. marinus populations from different geographic origins are adapted to the local tides in their exact circadian and circalunar emergence times (Kaiser et al., 2011;Neumann, 1967). These timing differences are particularly interesting for studying local adaptation and allochrony. They are precise within a population, but differ widely between populations. The correlation between the local time of spring low tide and C. marinus' emergence time is striking (Kaiser et al., 2011). Additionally, timing adaptations can be measured in the laboratory under common garden conditions, they are genetically determined and some of the underlying loci have been identified (Kaiser & Heckel, 2012;Kaiser et al., 2011Kaiser et al., , 2016Neumann, 1967).
In widely separated sites, the evolution of local timing adaptations is facilitated by geographic isolation, which largely arises due to the fact that C. marinus is restricted to rocky coasts (Kaiser et al., 2010). While populations in previous studies can be considered allopatric because of their occurrence in different stretches of rocky coast, we assumed that there must be parapatric Clunio populations within continuous stretches of rocky coast. The current study was motivated by the question: Can genetic timing adaptations persist in the absence of geographic isolation if ecological gradients are strong, i.e. along a continuously rocky coast with large differences in the timing of the tides? Consulting tide tables of the European Atlantic Coast, we identified the coastline of Brittany and Normandy in France as a suitable study area. The tidal amplitude in this region reaches up to 12 m, there is a steep gradient in timing of the tides, and the coastline is predominantly rocky. In a long-term project that started in 2009, we explored eight sites along the coastline (Figure 2).
We not only found the expected parapatric C. marinus populations with gradual timing adaptations, but also several sympatric populations which occupy divergent temporal niches and which we therefore introduce here as different timing types of C. marinus. As the English Channel did not exist during the last ice age (Patton et al., F I G U R E 1 Temporal niches and timing types of Clunio marinus. (a) Tidal amplitudes are bimodal across the lunar cycle, with lowest low tides occurring around full and new moon (bottom panel in a). Tides are also bimodal during the day, with two low tides a day (left panel in a). The time of low tide shifts across the lunar cycle (grey dots in central panel of a). Grey shading indicates night time. Superposition of the lunar cycle (x-dimension) and the daily cycle (y-dimension) results in four time points with extreme low tides (thick-lined squares in central panel of a), which are distinct temporal niches for reproduction of Clunio marinus. In C. marinus a circadian clock restricts reproduction to one of the two daily temporal niches, leading to Type 1 and Type 2 strains. A circalunar clock restricts reproduction to the lunar temporal niches, resulting in full moon strains (FM), new moon strains (NM), or strains emerging at both occasions (SL for "semilunar"). (b-e) Timing strains of C. marinus have specific combinations of circadian and circalunar timing. Strains can therefore be classified by temporal niche into the timing Types 1SL (b), 2SL (c), 2NM (d) and 2FM (e). Type 2NM strains generally emerge slightly earlier than optimal (d). Type 1NM or 1FM strains have not been found  (Figure 1a). Low tides are particularly low during the spring tide days, which occur just after full moon and new moon (bottom panel in Figure 1a). Additionally, there are two low tides per day (left panel in Figure 1a). The time of low tide shifts from day to day, but at full moon and new moon they always occur at about the same time of day (grey dots in central panel of Figure 1a).

Superposition of the lunar and daily temporal dimensions results in
four nonoverlapping temporal niches that are suitable for reproduction of C. marinus (thick-lined black boxes in the central panel of Figure 1a).
Timing strains of C. marinus are characterised by distinct, genetically determined combinations of circadian and circalunar timing, which synchronize reproduction to a suitable temporal niche. Timing strains were previously only reported for separate geographic sites and thus named by geographic origin. Here, we report timing strains that co-exist in the same site but reproduce in different temporal niches, which requires amending their nomenclature and distinguishing them by the temporal niche(s) they occupy.
When classified by temporal niche, the timing strains fall into four distinct groups, which we define here as ecological timing types. Firstly, as the circadian clock limits adult emergence of C. marinus to one of the two daily low tides (y-dimension in Figure 1a), we can distinguish timing strains that emerge during the first low tide after sunrise ("Type 1") from timing strains that emerge during the second low tide after sunrise ("Type 2"; see Figure 1a). Secondly, the circalunar clock limits adult emergence to a specific lunar phase (x-dimension in Figure 1a). Some C. marinus timing strains emerge only during full moon (lunar rhythm; "Type FM"), some only during new moon (lunar rhythm; "Type NM") and some during both (a so-called semilunar rhythm; therefore "Type SL"; see Figure 1a). Combining all possible circadian and circalu-  Table S1). We name timing strains based on a shorthand for their geographic origin and their timing type (Table S1). As the timing of the tides shifts gradually along the coastline, timing strains within each timing type differ slightly in circadian and circalunar emergence times in adaptation to the local tides.
On a side note, the SL type is not a mixture of FM and NM types nor the product of a genetic polymorphism. Laboratory experiments show that SL type strains have an endogenous circalunar clock with a ~15-day period, whereas NM and FM strains have circalunar clocks with a ~30-day period (Neumann, 1966). As a result, irrespective of whether SL type parents emerge during full moon or new moon, their offspring emerge during full moon and new moon in a roughly 1:1 ratio, over the course of 1-2 months.

| Fieldwork and establishment of laboratory strains
Field samples and laboratory strains for this study originate from five field trips undertaken between 2009 and 2017 (Table S2). From a total of eight sites we identified 11 different timing strains (Table   S2, Figure 2).
We applied two methods for sampling. First, we collected larvae by collecting small patches of algae and sand and transferred them Coloured dots represent the timing strains that were found in these sites. The colours correspond to the different timing types as described in Figure 1. In Plouguerneau and Roscoff several timing types co-exist in sympatry to the laboratory. Adults emerging from these samples reproduce readily and served to establish laboratory strains. Second, we sampled swarming adults in the field for genetic analysis. If these samples included copulae, we also established laboratory strains from their egg clutches. The reasoning behind this two-fold strategy and necessary deviations is described below.
We sampled larvae for all sites except for Concarneau and Camaret-sur-Mer (Table S2). As developmental time of C. marinus varies between 6-12 weeks, even for siblings (Neumann, 1966), generations of different timing types and timing strains overlap.
Consequently, at sites with sympatric timing types, larvae are present throughout the lunar cycle in ratios largely proportional to the population sizes of the different timing types. Collecting larval substrates allows us to identify all common timing types at a site, irrespective of the time point of collection.
In the case of sympatric timing types, establishment of laboratory strains from larvae regularly leads to mixed laboratory strains.
We applied two different strategies to obtain pure laboratory strains for each timing type. In cases where at least one of the temporal For type 2NM strains, which generally emerge a few days before new moon, the mean time of low tide at three days before new moon is given as an additional reference point (small grey arrows and lines). Field observations of swarming adults are given in local time as dotted circles. For the Plou-1SL strain only field observations are available F I G U R E 4 Circalunar timing of the timing strains found in Brittany and Normandy. Lunar emergence times were recorded from laboratory strains under common garden conditions (LD 16:8, 20°C). Artificial moonlight was presented during four nights (days 1-5) in a 30-day artificial moonlight cycle (x-axis). As biological rhythms are cyclic, day 30 is equivalent to the day before day 1. The bars show the fraction of emerging individuals on each day of the artificial moonlight cycle. The plots are coloured according to timing type. Within timing types the strains are sorted by geography, highlighting the gradual shift of circalunar emergence times along the coast, corresponding to the shift in the local tides. The small peaks around day 22 in the 2NM strains and day 18 in the 2FM strain are artefacts of the artificial moonlight treatment and have not been observed in the field (see Figure S1 for details) niches was exclusively occupied by one of the strains (Plou-1SL,  (Kaiser et al., 2011;Neumann, 1967). During extensive fieldwork, over the years we have only observed a handful of swarming adults at times during which such hybrids would swarm. In contrast, each day we have observed hundreds to many thousands of swarming adults during the pure timing types' temporal niches.
We conclude that timing types are distinct units in nature and that establishing pure timing-type laboratory strains is not an artefact of laboratory selection.
Field samples for genetic analysis consisted of swarming males caught on site. Females could not be found in sufficient numbers, as they are immobile and thus virtually invisible in the field. For Bria-1SL, Lou-2NM, Plou-2NM and Plou-2SL it was logistically impossible to collect adults in the field. Instead, field-caught larvae were transferred to the laboratory and the emerging adults were collected. For Bria-1SL and Lou-2NM this is unproblematic, as there is only one timing strain present in each site. For Plou-2NM and Plou-2SL this procedure was necessary because of the above-mentioned overlap in their emergence time. After rearing these field-caught larvae in the laboratory, for genetic analyses we used only those which had pure Plou-2SL type offspring or pure Plou-2NM type offspring respectively. To avoid sampling genetically related adults, in all cases several hundred larvae were collected from >10 different places in an area >5,000 m 2 and a random subset of 24 individuals was subject to genetic analysis.

| Characterization of laboratory strains
Timing adaptations, timing types and timing strains were assessed based on laboratory strains under common garden conditions.
Except for the Plou-1SL strain, all laboratory strains were maintained for a minimum of six, but usually >20 generations. They were kept in standard culture conditions (Neumann, 1966)  Lunar emergence times were recorded from the second laboratory generation onward by counting the number of adults that emerged daily. Counts were summed up over several artificial moonlight cycles. Daily emergence times were recorded from the third laboratory generation onward, after the strains had been identified based on their circalunar rhythms. Daily emergence times were recorded in 1-hr intervals with the help of a fraction collector (Honegger, 1977).
Emergence was recorded over several generations until for each strain >300 individuals had been assayed for circadian emergence times and >1,000 individuals had been assayed for lunar emergence times. The Ros-2FM and Por-1SL laboratory strains have been described (Kaiser et al., 2011;Neumann, 1966Neumann, , 1989. The re-established laboratory strains correspond to the original reports. Timing phenotypes were stable throughout the lifetime of a laboratory strain.

| Acquisition of genetic data
In order to test the evolutionary history and genetic structure of the timing strains and sympatric timing types, we obtained mitochondrial haplotypes and nuclear SNP data for 24 males from all timing strains (23 for Plou-1SL). DNA was extracted with a salting out method (Reineke et al., 1998), and subjected to whole genome amplification (REPLI-g Mini kit, Qiagen) according to the manufacturer's protocol.
For analysis of mitochondrial haplotypes we picked a fragment of cytochrome oxidase subunit I (COI), which discriminated C. marinus populations in previous studies (Fuhrmann & Kaiser, 2020;Kaiser et al., 2010). COI sequences were amplified with the primers C1-J-2183 and TL2-N-3014 (Simon et al., 1994)  (100-bp single-end reads). Only forward reads were available for all individuals, so that analysis was restricted to forward reads.
Raw data were submitted to ENA-EBI under the project accession PRJEB9361.
The reads were quality trimmed with cutadapt ( (Table S6).
With these genetic data we investigated the evolutionary history of the different timing types and strains in Brittany and Normandy, particularly the number of postglacial colonization sources, the extent of gene flow between them and the question of whether timing differences entail reproductive isolation.

| Population genetic analyses: Mitochondrial COI sequences
We inferred the number of maternal lineages that have colonized the study area by calculating a haplotype network. COI sequences were aligned in mega version 3.0 (Kumar et al., 2004) based on the clustal w algorithm (Thompson et al., 1994) and a haplotype network was calculated with network 4.6.1.1 software (Fluxus Technology Ltd), using the Median Joining method (Bandelt et al., 1999). Two multistate positions (i.e., positions with more than two alleles; positions 289 and 436) were excluded from the analysis.
We also assessed isolation by distance (IBD). Pairwise  (Table S3). We assessed correlation between the matrices of genetic differentiation and geographic distance by performing a Mantel test in arlequin version 3.5 with 1,000 permutations.

| Population genetic analyses: Nuclear SNP data
To assess the number of colonization sources and the genetic structure of populations we used principal component analysis (PCA) and the clustering algorithms implemented in STRUCTURE and ADMIXTURE. PCA was performed on the SNP data set using the R packages SNPRelate and gdsfmt (Zheng et al., 2012). Analysis was limited to 20 principal components, otherwise default options were used. The results were visualized in R (Crawley, 2007).
A script with the detailed commands for analysis and plotting is given as a supplement. Genetic admixture of populations was tested with structure 2.3.4 (Hubisz et al., 2009;Pritchard et al., 2000). Since we sampled eleven populations from eight localities, the test was carried out for K from 1 to 12. Admixed origin of populations was allowed (NOADMIX = 0). Due to the recent common origin of populations (see Section 3) allele frequencies were assumed to be correlated across populations (FREQSCORR = 1).
Burnin was set to 20,000 replications followed by 500,000 replications of data collection. All runs were iterated 10 times and convergence was assessed and visualized through the Cluster Markov Packager Across K (clumpak) web server (Kopelman et al., 2015), with default options. Additionally, best K was obtained from structure harvester (Earl, 2012) by the method of Evanno et al. (2005). Given the weak population structure (see Section 3), we also performed structure runs with location priors (LOCPRIOR = 1), but neither the model fit (LnP(K)) nor the cluster assignments changed notably (data not shown). For comparison with structure, we ran admixture (Alexander & Lange, 2011) with increased convergence stringency (-C 10 −7 ), 10-fold cross-validation (--cv=10) and otherwise default parameters. Results were again visualised through the clumpak web server.
In order to specifically assess gene flow between timing strains we finally performed structure runs for K = 11 with fixed population information (USEPOPINFO = 1), so that only the mixing between pre-defined population-specific genetic clusters was assessed. All other parameters were unchanged. We assessed convergence of the 10 replicates by calculating the standard deviation in the cluster assignments to each population (Table S4).
We further tested for gene flow between the populations by assessing if there was isolation by distance (IBD To further investigate gene flow, we attempted to detect migration events with TreeMix 1.13 (Pickrell & Pritchard, 2012). We converted the SNP data to TreeMix format with the vcf2treemix.
sh script (Ravinet). For 0 to 10 migration events, we ran five iterations each of TreeMix, with no blocks (as SNPs were already thinned for LD) and the southernmost population from Concarneau as root to the tree. As a second approach for testing migration, we tried to detect migrants in our samples from sympatric populations.
Distinguishing genetic markers were obtained from the largely undifferentiated sympatric timing strains by filtering our SNP set to those with F ST > 0.2 between sympatric timing strains (n = 12 for Roscoff; n = 88 for Plouguerneau). We then performed PCA on these subsets of SNPs as described above. Individuals clustering in a different cluster than expected could represent migrants. Genetic differentiation between timing strains was too low to obtain diagnostic markers for identifying hybrids via hybrid index and interstrain heterozygosity (data not shown).
In order to assess the contribution of geographic distance and/ or timing differences to reproductive isolation, we tested whether genetic differentiation between the populations was influenced by these factors. The 11 timing strains were grouped according to geography or emergence timing (Table 1) and analysis of molecular variance (AMOVA) was performed for the SNP data set in arlequin version 3.5 (Excoffier & Lischer, 2010) with default settings and 1,000 permutations. From these arlequin runs we also obtained nucleotide diversity π per population for the SNP data set (Table S5).
Assuming that some of the populations in our sample may be very small or young, we assessed relatedness of individuals in our sample by calculating the inbreeding coefficient F with the --het command in vcftools 0.1.14 (Danecek et al., 2011). The results were visualized as population-specific box plots in r (Crawley, 2007).  Figure 1a). These temporal niches occur at a specific lunar phase (bottom panel Figure 1a) and a specific time of day (left panel in Figure 1a). Timing strains of C. marinus are characterised TA B L E 1 Groups of populations formed for testing genetic differentiation due to timing differences or geographic site in AMOVA timing strains that emerge during the first low tide after sunrise ("Type 1") from timing strains that emerge during the second low tide after sunrise ("Type 2"; see Figure 1a). Along the circalunar time axis, timing strains emerge either only during full moon (lunar rhythm; "Type FM"), or only during new moon (lunar rhythm; "Type NM") or during both (a so-called semilunar rhythm; therefore "Type SL"; see Figure 1a). Combining the two time axes results in six possible timing types, but only 1SL, 2SL, 2 NM and 2FM strains have been found (Figure 1 b-e; Table S1). In order to distinguish sympatric timing strains, we amend the previous convention of naming timing strains based on a shorthand for their geographic origin by adding their timing type (Table S1). As an example, "Ros-2FM" is a type 2FM strain from Roscoff. As the timing of the tides shifts gradually along the coastline, timing strains within each timing type differ slightly in circadian and circalunar emergence times in adaptation to the local tides.

| In Brittany and Normandy C. marinus timing types occur in sympatry
In this study, we explored the capabilities of C. marinus to genetically adapt to divergent timing of the tides on small geographic scales by studying eight locations along the coasts of Brittany and Normandy ( Figure 2). In these eight sites we identified 11 distinct C. marinus timing strains (Figure 2), which represented all four known timing types of C. marinus, occupying all four available temporal niches (compare Figure 1 b-e with Figures 3 and 4). We characterized their daily and lunar emergence times under common garden conditions in the laboratory (Figures 3 and 4). The Plou-1SL laboratory strain persisted only briefly, so that for this timing strain circalunar data is limited and only field observations are available for the circadian rhythm ( Figure 3).
In two locations, Roscoff and Plouguerneau, different timing types co-occurred (Figures 2 to 4), which is the first description of sympatric timing types for C. marinus. In Roscoff, the Ros-2FM and Ros-2NM strains are fully separated by lunar timing, one emerging only during full moon, the other only during new moon (Figure 4).
The smaller second peaks around day 18 in Ros-2FM and around day 21 in Ros-2NM (Figure 4) are an artefact produced by the artificial laboratory moonlight (for details see Figure S1). During field work we have not observed individuals swarming at these times of the lunar month. In Plouguerneau three timing strains co-occur: Plou-1SL is fully separated from the other two by circadian emergence time (red dotted circle, Figure 3). Of the two type 2 populations, Plou-2NM emerges only during new moon, while Plou-2SL emerges during both full moon and new moon, which leads to a considerable overlap in emergence times (Figure 4).

| Assessing local adaptation to the timing of the tides
As the timing of the tides changes gradually along the coastline, we expect that within each timing type the different parapatric timing strains should be locally adapted. Indeed, within timing types, the lunar emergence peaks shift gradually when moving along the coastline (Figure 4), corresponding to the shift in the local timing of the tides. Daily emergence times also shift with the local time of the low tides ( Figure 3). However, for some timing strains, daily emergence times appear suboptimal relative to the local tides (Figure 3). 2NM type strains generally emerge about three days before new moon, when the low tide is 2-3 h earlier (see Figure 1d). When correcting for this offset in circalunar timing by plotting the daytime of low tide at three days before new moon (grey arrows in Figure 3), the match of 2NM type strains' emergence with the local tides is much better. In some timing strains laboratory emergence times and field emergence times diverge, most notably in the Eta-1SL strain (compare histograms with dotted circles; Figure 3). In the field, circadian emergence time may depend not only on the circadian clock, but also additional external stimuli (see Section 4).

| All timing types emerged from a single recent colonization event
The existence of sympatric timing types and locally adapted timing strains raises the questions of how they have evolved and how they are maintained, particularly whether sympatric ecological divergence in the presence of gene flow is plausible. In order to address these questions, we genetically analysed 24 individuals for each of the 11 populations based on mitochondrial COI haplotypes and a nuclear SNP data set obtained by RAD sequencing (summary statistics in Table S6).
In a network of the mitochondrial COI haplotypes all 11 populations and 73% of all individuals share the same major haplotype Next, we assessed population structure based on the SNP data with principal component analysis (PCA, Figure 6a) and the structure clustering algorithm (Figure 6b). Populations across Brittany and Normandy show limited genetic structure. In particular, the sympatric timing strains from Roscoff and Plouguerneau, which comprise all four timing types, appear highly similar in genetic composition in both analyses (Figure 6a and b). There is no apparent genetic distinction between the four timing types, again dismissing the hypothesis that genetically distinct and reproductively isolated timing types may have colonized Brittany independently. In conclusion, the diversification into the four timing types and the establishment of local timing adaptations must have occurred within Brittany and Normandy.

| Geographic structure and gene flow
The with groups formed according to geographic location attributes 5.96% of the variation to geography (Table 2). Similarly, genetic clusters obtained with structure weakly follow geography, especially at low numbers of clusters (Figure 6b). The most likely number of groups in the sample, as determined with the ΔK method, are K = 2 and K = 8 ( Figure S2). Results from an additional admixture analysis are very similar to those obtained with structure and have the lowest cross-validation error at K = 8 ( Figure S3). Taken together, there is weak geographic structure imposed by the eight sampling sites.
Next, we fixed the population information in the structure analysis, so that we could assess the mixing of pre-defined, population-specific clusters. In this analysis, each timing strain corresponds to a specific genetic cluster (Table S4). Additional genetic components come from sympatric timing strains or from adjacent geographic sites, which suggests that the weak geographic structure is a product of pervasive local gene flow along the coastline. We substantiated this finding by calculating pairwise genetic differentiation (F ST ; Table S7)  In contrast to the SNP data, there is no IBD in the mitochondrial COI sequences (Table S7, Figure S4). Mitochondrial gene pools show highly variable differentiation, irrespective of geographic scale, suggesting they are geographically isolated and diverge at random by genetic drift. This finding is in line with the fact that C. marinus females are wingless and basically immobile. Asymmetry in nuclear vs.
mitochondrial differentiation indicates that gene flow along the continuous rocky coast of Brittany and Normandy is primarily mediated by swarming adult males.
In order to test for migration events between the populations, we used the TreeMix algorithm, which assesses whether the observed correlation of allele frequencies between populations is better F I G U R E 5 Mitochondrial Cytochrome oxidase I (COI) haplotype network. (a) Haplotype network of an 853 bp fragment of COI for 262 individuals from eleven populations in Brittany and Normandy. Lines indicate single nucleotide polymorphisms (SNPs) separating the haplotypes. Size of the circles corresponds to the fraction of individuals carrying the haplotype. Origin of individuals is coded by geography and timing type (see legend). (b) Haplotype network of the same COI fragment for C. marinus strains from neighbouring European coasts, adapted from (Kaiser et al., 2010). Small red lines indicate SNPs. Colour of the haplotypes represents the timing type (not known for Ireland). The haplotypes from Normandy (highlighted by black lines around the circles) include the Por-1SL population, which was sampled again for this study. Across the two studies, the major haplotype and many minor haplotypes are identical. Populations from other European coasts have entirely different haplotypes

| Reproductive timing is not an effective reproductive isolation factor in C. marinus
For speciation to occur, ecological divergence between populations must be accompanied by the build-up of reproductive isolation. Given the large differences in reproductive timing of C. marinus timing types, we would assume that allochrony must play a role in this species. In order to assess to what extent temporal isolation results in genetic divergence between the timing types, we applied an analysis of molecular variance (AMOVA). Populations were F I G U R E 6 6Genetic structure and isolation by distance between timing strains of C. marinus in Brittany. Analyses are based on SNP data from restriction-site associated DNA sequencing (RAD-seq; 2,159 unlinked SNPs). (a) Principal component analysis (PCA) for 263 individuals from 11 populations. Different symbols represent timing strains, colours correspond to timing types as defined in Figure 1. A lack of clustering according to timing types suggests that timing types are not genetically distinct, reproductively isolated units. There is limited clustering according to geography. The small amount of variation explained by the principal components indicates that population structure is overall weak. (b) structure analysis for 2-12 genetic clusters (K). Results from 10 independent replicates are summarized in clumpak. The numbers of replicates that converged are given to the left of each panel ("x/10"). Replicates that did not fully converge only differ in small details (see Figure S7). There is no separation by timing type and limited geographic structure. (c) Linearized pairwise genetic differentiation (F ST ) between timing strains plotted against geographic distance indicates isolation by distance (IBD). Sympatric timing types are the least differentiated grouped according to circadian timing phenotype, circalunar timing phenotype or geography (Table 1), and we assessed which fraction of genetic variation is associated with these groups. In all three comparisons, the major component of genetic variation is within individuals (95.83%-96.62%) and there is some genetic variation between populations within each group (2.08%-7.71%). In addition to these general observations, 1.59% of genetic variation is associated with groups based on circadian timing ( Still, the results obtained in PCA on high F ST SNPs for the sympatric Plouguerneau timing strains, i.e., without any possibility for geographic confounding, also suggests that circadian timing entails more genetic structure than circalunar timing (Figure 7b). The sympatric strains from Plouguerneau are well-separated by circadian timing (PC1; 37.6% of variation), but much less separated by circalunar timing (PC2; 9.7% of variation). In summary, a small but detectable fraction of genetic variation is associated with and possibly caused by differences in circadian timing. But overall, there is no strong genome wide genetic divergence that would be attributable to circadian allochrony or circalunar allochrony.

| Local adaptation
In Brittany we found different sympatric and parapatric C.  emergence times of the different timing strains, which exist within each timing type, must be considered local adaptations. They were shown to match the local timing of the tides extremely well, with correlation coefficients of up to 0.98 (Kaiser et al., 2011). They are genetically determined and stable under laboratory conditions without any tidal timing cues (this study; Kaiser et al., 2011;Neumann, 1967).
The circalunar timing phenotypes presented here (Figure 4) fulfil the expectations for local adaptation. Within each timing type they show a gradual shift along the coastline corresponding to the strong gradient in the timing of the tides. For circadian timing the match between local tides and laboratory emergence time is still good, but shows exceptions (Figure 3). There are two major reasons for this. First, the 2NM type strains generally emerge a few days before new moon when the low tide occurs 2-3 h earlier ( Figure 1d). When correcting for this effect, the match between local low tide and laboratory emergence is much better (Figure 3).
Variation in the daytime of the low tide is larger three days before new moon compared to new moon (compare standard deviations of black and grey arrows in Figure 3). This might explain the comparatively broad spread of circadian emergence times in 2NM type strains ( Figure 3). Second, this study is the first for which emergence times have not only been recorded in the laboratory, but have also been observed in the field for many of the timing strains.
In some cases, emergence times in the field and laboratory differ notably ( Figure 3). The genetic timing adaptations, which are recorded in the laboratory without any tidal timing cues, may define a site-specific circadian emergence window within which environmental cues can trigger emergence and thereby further finetune its timing. Such phenotypic flexibility is advantageous, as the emergence of C. marinus takes place over several days during which the daytime of low tide shifts by several hours. Allowing for emergence to be triggered by tidal cues within the circadian emergence window would enable C. marinus to track this shift. Unlike other intertidal organisms, C. marinus does not possess a 12.4 h circatidal clock to achieve such tracking (Neumann, 2014).
Our field and laboratory observations indicate that environmental cues for directly triggering emergence might include light levels and exposure of the larval habitat by the tide. In the Plou-2NM, Ros-2NM and Ros-2FM laboratory strains most females emerge immediately when the light is turned off (Figure 3). In Eta-1SL field emergence was highly concentrated at the time when the larval habitat became exposed by the tide, whereas daily emergence in the laboratory is extremely broad and nonoverlapping with observed emergence in the field. Interestingly, the Eta-1SL strain has the lowest genetic diversity (Table S5)   Secondly, it is possible that timing types have evolved in separate sites in Brittany with some reduction in gene flow due to IBD, and then secondarily spread along the coast to become fully sympatric.
Given the current strong gene flow along the entire coast, this would still represent a case of divergence with gene flow, but not starting from initial panmixia. Finally, it is possible that at least some of the alleles responsible for sympatric divergence in timing predate the postglacial colonization event and upon arrival in Brittany facilitated local divergence with gene flow. The source of pre-existing timing alleles may either be standing genetic variation in the colonizing population or yet undetected introgression events. The scenario of standing genetic variation is consistent with a previous study on five allopatric timing strains, which indicated that standing genetic variation at adaptive timing loci is common in C. marinus (Kaiser et al., 2016). The scenario of undetected introgression seems conceivable in the light of C. marinus' modes of dispersal (Kaiser et al., 2010). Adult dispersal is common, but geographically limited due to the short adult life span. As females are wingless, adult dispersal is largely mediated by males, resulting in strong geographic isolation of mitochondrial gene pools as also detected here ( Figure S4). Dispersal of the benthic larvae via drifting algae appears to be rare, but can be long distance and is not expected to be sex-biased. If we assume Distinguishing between these scenarios will be very difficult given the strong gene flow and low degree of genetic differentiation. It will certainly require additional sampling and comprehensive genomic data. But taken together, all of these scenarios imply a mixing of timing types and divergence with gene flow, which poses questions to which role temporal isolation and ecological selection by the tides play in forming and maintaining sympatric timing types of C. marinus.

| Reproductive isolation and natural selection
Many prominent examples of divergence with gene flow are characterized by differential habitat or host use resulting in heterogeneous distributions within the range Filchak et al., 2000;Nosil et al., 2006;Schliewen et al., 2001), sometimes referred to as mosaic sympatry (Mallet et al., 2009) or heteropatric speciation (Maynard Smith, 1966). This does not seem to be the case for C. marinus, where we have isolated different timing types from the very same patch of algae.
In the light of C. marinus' biology and life history, reproductive timing must be expected to result in prezygotic reproductive isolation through allochrony and should eventually constitute a magic trait for speciation (Gavrilets, 2004). Given the short reproductive period of C. marinus of only a few hours, the large timing differences between sympatric timing types should minimize or even abolish any overlap in reproductive time. Our genetic analyses of sympatric timing types suggest that mating is certainly not random between timing types (Figures 6 and 7). In Roscoff we detected two or more putative migrants in a sample of 48 individuals, suggesting there is strong migration, but not panmixia. The correlation between timing differences and genetic divergence of populations is marginal for circadian timing and nonexistent for circalunar timing (Table 2). Clearly, prezygotic isolation by timing does not result in genome-wide genetic divergence.
The lack of isolation by circalunar timing may partially result from the existence of SL type strains, which emerge at full moon and new moon and may mediate gene flow between FM and NM strains. However, the strong gene flow in Roscoff is still surprising, as in the absence of an SL type strain, emergence of the sympatric timing strains at full moon vs. new moon results in no overlap in reproductive time. Field observations from Helgoland provide a potential explanation for this puzzle. In Helgoland, spring emergence of overwintering individuals was found to be triggered by the warming sea water, which overrode the effect of the circalunar clock (Krüger & Neumann, 1983). If this is also true for our study sites in Brittany, sympatric timing types could meet and interbreed every spring, i.e., every 3-4 generations, circumventing the temporal isolation imposed by their lunar timing phenotypes.
Another source of reproductive isolation may lie in selection against timing type hybrids, which are usually intermediate in timing between their parents (Heimbach, 1978;Kaiser et al., 2011;Neumann, 1967). We expect that Ros-2NM × Ros-2FM hybrids would emerge at neap tide high tides, a situation most unsuitable for reproduction. The situation in Plouguerneau would be the same, though slightly more complex given that there are three timing types. Thus, in the absence of notable temporal or geographic isolation, the existence and maintenance of distinct sympatric timing types probably relies on permanent strong selection against hybrids, an ecologically imposed form of extrinsic postzygotic isolation. As such an isolation mechanism drastically reduces the fitness of parents producing hybrids, we may expect strong selection for prezygotic isolating factors that would reinforce ecological divergence.
Assortative mating based on other factors but timing is a possibility to consider in future research. A recent study on the olfactory system of C. marinus suggests that pheromone communication exists in this species (Missbach et al., 2020).

| Perspectives
Ultimately, studies on local adaptation and sympatric population divergence will be most powerful when the loci underlying the ecological adaptations of diverging populations are identified. For C. marinus, the populations described here may be key to identifying these loci, as the observed strong gene flow is expected to homogenize their genomes except for the few ecologically relevant loci.
The current RAD sequencing does not have the resolution to identify such loci, as RAD tags in this study are on average 24 kb apart, while even in more divergent C. marinus populations genomic divergence peaks usually extend over 2-5 kb only (Kaiser et al., 2016).
But genome scans and QTL mapping have proven a powerful tool in C. marinus. They may not only help to find genes involved in local adaptation of biological clocks and sympatric population divergence with gene flow, but may also give first insights into the yet enigmatic circalunar clockwork.