A novel integrative approach elucidates fine-scale dispersal patchiness in marine populations.

Dispersal is one of the main determining factors of population structure. In the marine habitat, well-connected populations with large numbers of reproducing individuals are common but even so population structure can exist on a small-scale. Variation in dispersal patterns between populations or over time is often associated to geographic distance or changing oceanographic barriers. Consequently, detecting structure and variation in dispersal on a fine-scale within marine populations still remains a challenge. Here we propose and use a novel approach of combining a clustering model, early-life history trait information from fish otoliths, spatial coordinates and genetic markers to detect very fine-scale dispersal patterns. We collected 1573 individuals (946 adults and 627 juveniles) of the black-faced blenny across a small-scale (2 km) coastline as well as at a larger-scale area (<50 kms). A total of 178 single nucleotide polymorphism markers were used to evaluate relatedness patterns within this well-connected population. In our clustering models we categorized SHORT-range dispersers to be potential local recruits based on their high relatedness within and low relatedness towards other spatial clusters. Local retention and/or dispersal of this potential local recruitment varied across the 2 km coastline with higher frequency of SHORT-range dispersers towards the southwest of the area for adults. An inverse pattern was found for juveniles, showing an increase of SHORT-range dispersers towards the northeast. As we rule out selective movement and mortality from one year to the next, this pattern reveals a complex but not full genetic mixing, and variability in coastal circulation is most likely the main driver of this fine-scale chaotic genetic patchiness within this otherwise homogeneous population. When focusing on the patterns within one recruitment season, we found large differences in temperatures (from approx. 17 °C to 25 °C) as well as pelagic larval duration (PLD) for juveniles from the beginning of the season and the end of the season. We were able to detect fine-scale differences in LONG-range juvenile dispersers, representing distant migrants, depending on whether they were born at the beginning of the season with a longer PLD, or at the end of the reproductive season. The ability to detect such fine-scale dispersal patchiness will aid in our understanding of the underlying mechanisms of population structuring and chaotic patchiness in a wide range of species even with high potential dispersal abilities.

Introduction

In marine ecosystems, larval dispersal and recruitment are the main processes determining population structure of species and understanding how they maintain connectivity is crucial for effective management[1]. Although our understanding of larval dispersal in the marine environment has greatly improved over the last several decades[2], much is still unknown. Most marine organisms exhibit a larval phase where small larvae can, in theory, travel large distances with ocean currents. Marine populations were thus generally thought to have elevated gene flow and little population structure[3]. This paradigm has been largely disproven as many recent studies have revealed that genetic population structure can exist across even small spatial scales and self-recruitment is often larger than previously believed[4,5]. In addition, direct measurements of dispersal using pedigree reconstruction analyses have revealed that self-recruitment and recruitment patterns are often heterogeneous and complex[6,7].

Dispersal patterns depend on pelagic larval duration (PLD)[8], reproductive output[9], ocean currents[10] or localized eddies and oceanographic retention patterns[11,12]. The duration of the larval period is known to influence larval dispersal distance[13] and can have a direct relationship with connectivity[9,14]. Information on PLD has been used in numerous studies modelling dispersal capabilities, connectivity and for the establishment of networks of Marine Protected Areas[7,15-17]. Several studies have shown an influence of ocean temperature on PLD, decreasing its duration with increasing temperature[13], thus a change in temperature could have a direct influence on population connectivity[18]. Considering the potential increase of sea temperatures due to global climate change[19], studies analysing the relationships among temperature, PLD and connectivity are essential in marine ecology.

In the last few years, studies of connectivity have used observations of kin, such as parent-offspring pairs or full siblings, to infer the trajectories and distances of larval movement at small scales, but require large sampling efforts and genotype datasets[7,17,20-23]. Further factors determining patterns of dispersal are availability of habitat, environmental selection, behaviour or kinship[24]. The population structure of marine organisms can be influenced by the degree to which larvae from different populations or demes are mixed in the plankton[25,26]. Some studies have demonstrated the existence of familial structure with marked genetic distances between recruitment cohorts; thus the lack of larval mixing during the planktonic period can lead to chaotic genetic patchiness[27-29]. Other studies, however, showed that related individuals did not necessarily disperse cohesively, but limited dispersal was the cause of non-random patches of closely related individuals[30]. Hence, many questions remain to be answered about dispersal, settlement and survival patterns of related individuals, how these are affected by environmental conditions such as temperature and how in turn they affect the structure of populations in the marine environment. Particularly genetic structure on a very fine-scale has not been investigated and we here present a novel approach to understanding recruitment patterns leading to population structuring.

In the last decade, many fields of ecology have experienced a data revolution[31,32], molecular ecology not being an exception. The existence of large-scale structured datasets with fine scale metrics related to environmental data, behavioural phenotypes, or genetic information, challenges traditional ways to analyze data structure. Unsupervised clustering methods[33-36] are useful to get insights about data structuration when there is minimal or no information about it. Several first principled methods can be used to highlight structural properties, some of them particularly focusing on the interpretability of such structures[37]. Here, we apply novel and easy-interpretable clustering procedures to genetic relatedness information in order to infer fine-scale genetic structuring and dispersal in large populations.

We aim to detect fine-scale structure and dispersal patterns of a rocky shore fish species, the black faced blenny (Tripterygion delaisi), using patterns of kinship inferred from single nucleotide polymorphism (SNP) genotypes. We used a novel approach of binary data clustering to consider these fine-scale dispersal patterns. Specifically, we combined larval traits, such as pelagic larval duration, and evaluated selective dispersal patterns via genetic marker heterogeneity in space. In addition, the detected fine-scale patterns allowed studying the effect of temperature on dispersal throughout the recruitment season with lower temperatures at the beginning and higher temperatures at the end. In summary, this approach is used to detect fine-scale patterns of dispersal and potential spatial assemblages of recruits, as well as differences between early and late season individuals in a well-connected population of T. delaisi, revealing patterns of relatedness that were not previous detected using traditional methods.

Methods

Study system and sample collection

T. delaisi is a small, common, rocky-shore fish from the Mediterranean Sea and the eastern Atlantic[38]. Individuals of this species can live to a maximum of three years and adults display high territorial fidelity[39,40]. Furthermore, 1-year old individuals are the most abundant component of the reproductive population. The spawning period starts in April, when water temperatures in the study area are ca. 15 °C and finishes at the end of July, with temperatures of ca. 23 °C[41]. The planktonic larval duration (PLD) of this species was estimated to be about 2–3 weeks, depending on the water temperature[42]. Therefore, this species is a good study organism to analyse the relationships between temperature and dispersal patchiness.

Overall, 946 adult and 627 juvenile T. delaisi were collected on SCUBA from the coast of Blanes, Spain in the western Mediterranean Sea (41°40′30.7″N 2°48′14.4″E). Adults were caught with hand nets, body length was measured, and dorsal fin clips were taken non-lethally underwater, and then preserved in 95% ethanol. Adults were collected during the length of their breeding season from April 2010 to July 2010. Juveniles were also captured with hand nets on SCUBA between July and September 2010 and subsequently used for genetic and otolith analysis. Sea surface temperature readings were taken 5 m below the surface twice weekly with a Conductivity, Temperature and Depth probe (CTD) at Las Medes Marine Reserve (approximately 40 km north of the focal sampling area in Blanes). Collection and field procedures followed the Spanish Laws (Royal Executive Order, 53/2013) for Animal Experimentation, in accordance with the European Union directive (2010/63/UE) and approved by the Consejo Superior de Investigaciones Cientificas (CSIC).

Two different spatial scales were studied. A ~2 km stretch of coastline in the vicinity of Blanes (Fig. 1) is the small-scale area in which exhaustive searches were performed for all dominant adult males protecting small nests (N = 793), and randomized searches for non-dominant males and females (N = 153), as well as juveniles (N = 382). The large-scale area encompasses the coast towards the southwest and northeast of the small-scale area, spreading over a total of 42 km, where juveniles (N = 245) were collected (Fig. 1). Juveniles were collected in this large-scale area from La Pilona (2.7922, 41.6701; N = 1), Palomera (2.8078, 41.6781; N = 62), Blanes north (2.8181, 41.686; N = 13), Clotilde (2.8214, 41.6887; N = 23), Lloret de Mar south (2.8402, 41.6929; N = 50), Lloret de Mar north (2.86125, 41.699; N = 23), Tossa de Mar (2.934, 41.7146; N = 74). The exact location of collection of every sample in both areas was geo-referenced.

Figure 1

Map of sampling area including the small-scale intensive sampling area (Blanes) and the large-scale collection area. (1) La Pilona, (2) Palomera, (3) Blanes north (4) Clotilde, (5) Lloret de Mar south, (6) Lloret de Mar north, (7) Tossa de Mar. The map was created using ArcGIS 10.5 (www.arcgis.com).

Otolith readings

Juvenile otoliths (lapilli) were extracted and mounted on microscope slides using thermoplastic glue (Crystalbond 509). To expose all daily increments, the otoliths were mounted and polished using two different grained sandpapers (3 μm and 1 μm Imperial lapping film, 3 M) to obtain a thin section exposing the nucleus and the first growth ring. Reading of otolith bands was performed using a high-powered microscope with transmitted light (AxioPlan, Zeiss) connected to a ProgRes C10 camera (Jenoptic) and an image analysis system[43]. To verify the first growth increment of the otolith, fertilized T. delaisi eggs collected in the field were reared in aquaria. First growth bands were checked at day 1 and day 3 post-hatching in ten individuals. This yielded both the timing of formation of the first band and the daily pattern of band deposition and confirmed that band deposition took place daily from the day of hatching. Otolith marks in this species were always clear and belonged to Type Ia: an abrupt settlement mark with a sharp decrease in increment width across the settlement mark, completed within one increment[41,44]. The length of pelagic larval duration (PLD) for each individual was determined by counting the number of daily rings visible between the core and the settlement mark, and the age at sampling was determined by counting the total number of bands from core to margin. We analysed the otoliths along the longest radius from the centre to the edge and recorded otolith width during the pelagic larval period. We also recorded the size-at-hatching (radius from core to first increment) and size-at-settlement along the longest axis of the otolith. To minimize errors, all measurements were repeated three times. Trait differences were statistically evaluated in R (v. 3.3.1; RCoreTeam 2016)[45].

A total of 200 randomly selected otoliths from juveniles collected in the small-scale area were read. Total body length of the juveniles correlated very well with the determined age by otolith reading (Pearson’s R = 0.9151, p < 0.001, Supplementary Material Fig. S1), so the age of all 627 juveniles was estimated from body length with the regression equation. We inferred larval traits such as the hatching date for all juveniles by combining the age and the collection date of each individual. Inferred larval traits (size at hatching, size at settlement and pelagic larval duration) were compared between subsets of particular juvenile individuals (e.g. early vs late juveniles see below) with otolith readings by t-tests in R. To compare the beginning and the end of the reproductive season, juveniles were assigned to two groups: ‘early juveniles’ if born between April and mid-May (N = 164 within Blanes and N = 224 for the large-scale area) and ‘late juveniles’, born between the 15th of June and the end of the recruitment season, which is at the beginning of September (N = 121 in Blanes and N = 229 for the large-scale area). One additional comparison was made to be able to evaluate differences in possible selective mortality post-settlement. We could establish individuals that had just settled (<30 mm in size at collection and <30 days total age) named “recent settlers” and individuals that had survived a period of minimum of 30 days after settlement (>35 mm in size and >50 days total age) considered “survivors”.

Genetic analyses

All 1573 individuals were genotyped with 192 SNP markers developed for T. delaisi; details of the SNP development can be found in Schunter et al.[46]. Markers that were not in Hardy-Weinberg or linkage equilibrium were discarded and the remaining 178 SNP loci used in subsequent population genetic and parentage analyses, as described in Schunter et al.[20]. Heterozygosities for each marker and all markers combined were computed with GenAlEx v. 6.5[47]. Pairwise relatedness values (r-values) among adults as well as among juveniles were obtained by using the Queller and Goodnight equation[48] implemented in the SPAGeDi software[49]. To find patterns of high or low relatedness amongst thousands of pairwise values, we used a relatedness ‘ratio’ to identify individuals with large numbers of high relatedness values. This is not a direct measurement of how these individuals recruited, but an attempt to categorize individuals according to their relatedness values to other samples and find patterns within this well-connected population. For each individual, the number of r-values above 0.1 was counted and divided by the total number of comparisons. The top 25% of individuals with the highest ratios were termed ‘locals’, as these individuals are showing the highest relatedness values to a larger portion of other individuals. The samples with the lowest ratios (bottom 25%), which have an elevated number of highly negative relatedness values with the rest of the individuals, were referred to as ‘recent migrants’. A cutoff of 10% was also used for comparison, however, patterns did not differ although the number of samples was small (data not shown). Differences in spatial structure of these groups of individuals were evaluated using their geo-referenced locations of collection and with a Kolmogorov-Smirnov test in R.

Spatial and genetic clustering

To evaluate the spatial distribution of kinship we incorporated spatial structuring of samples into the analysis. Latitudes and longitudes for each individual’s sampling location were used to find spatial clusters within the small-scale (Blanes) area as well as in the large-scale sampling area and assigned individual geo-localized samples to the respective spatial clusters. To this end we used the R package EMCluster[35]. Based on the overall depiction of the positions, we assumed a Gaussian mixture model over the spatial distribution of individuals and performed several runs of the algorithm with 4, 5 and 6 initial Gaussian components (10 runs each) picking the model with lowest AIC[50] which resulted in 6 spatial clusters.

To take into account the spatial signature of relatedness estimates, we averaged individual genetic relatedness considering their spatial clustering across the study area. Specifically, we computed the genetic relatedness of each individual with each of its conspecifics in the same spatial cluster and with all other individuals (see Fig. 2). As pairwise relatedness values can be negative when using an unbiased estimator, especially in continuous large populations[51] for calculation needs we converted each single relatedness value into a measure between 0 and 1 with the formula (d = 1 − (r − a)/(b − a)), where a and b are the minimum and maximum value of relatedness (i.e., r), respectively. Thus, for each individual we computed two distinct values based on mean pairwise-relatedness: (i) the mean relatedness value of an individual in a given spatial cluster to all the other members of the same cluster (rIN), and (ii) the mean relatedness value of that individual to the rest of individuals outside of its spatial cluster (rOUT). Finally, we characterized each individual with these two measures and performed an Expectation-Maximization binary clustering with the R package EMbC[37]. The EMbC splits the range of each variable into high and low values based on a Gaussian mixture model. The thresholds of these binary partitions are determined by the frontier lines of equiprobability between the different Gaussian components. This analysis allowed us to partition all individuals according to four categories (Fig. 2): 1. MIXED-range dispersers: individuals with high mean relatedness (r) with respect to individuals in their own spatial cluster (rIN) and high mean r with respect to individuals in other spatial clusters (rOUT), which are most likely individuals with relatively common genotypes that confer medium relatedness values with many other individuals. 2. SHORT-range dispersers: individuals with high rIN and low rOUT, which therefore are the individuals that have high relatedness to neighboring individuals. 3. MEDIUM-range dispersers: individuals with low rIN and high rOUT, which are fish that came from the area but from a different spatial cluster. 4. LONG-range dispersers: individuals with low rIN and low rOUT, which means that these individuals have very low relatedness with most of the other samples and could be long-distance dispersers from a different genetic population.

Figure 2

Conceptual Plot showing the four categories of dispersers. Two spatial clusters are represented and the mean relatedness value (r) among the individuals within one cluster (rIN) indicated with a circular arrow and with individuals of other spatial clusters (rOUT) with a straight arrow. The thickness of the lines represents the magnitude: thick = high r, thin = low r.

Results

Genetic and phenotypic diversity

We were able to determine the hatch date of each individual juvenile and the number of days the larvae spent in the water column with a subset of 200 randomly selected otoliths. This allowed for identifying some life history traits of the collected juvenile fish. The amount of time a larva spent in the water column depended on when in the reproductive season it was born (Pearson’s r = −0.6057, p < 0.001; Fig. S2a). Individuals hatching in April can have a planktonic tenure of up to 10 days longer than those hatching in July. No significant differences were found between ‘early juveniles’ and ‘late juveniles’ for size-at-hatching (t = 1.73, p = 0.086) or size-at-settlement (t = −1.009, p = 0.32). Average PLD was 19.5 days for early juveniles and 16 days for late juveniles and were significantly different (t = 9.19, p = 2.596e-14). The data were found to be normally distributed and no significant differences were found between recent migrant and local juveniles, as categorized based on the relatedness ratio nor for any otolith characteristic (size-at-hatching t = 0.57, p = 0.56; size-at-settlement t = 0.25, p = 0.80 (Fig. S2b); PLD t = −0.08, p = 0.94).

To test for possible sweepstake effects from one recruitment year to the next, as well as from the beginning of the recruitment season to the end, we compared average heterozygosity between all juveniles and adults as well as early and late juveniles. We found no significant differences in heterozygosity neither between juveniles and adults (t = −0.53, p = 0.59) nor between early and late juveniles (t = −0.95, p = 0.34) for all loci combined (Fig. S3), nor between recent settlers and survivors (t = 0.93 p = 0.35). Small differences in heterozygosity between early and late juveniles were found when analysed locus by locus (Fig. S4). Relatedness did not differ between adults and juveniles (t = 1.301, p-value = 0.194; Fig. S5). Due to a longer PLD at the beginning of the season, individuals could potentially travel further distances and disperse further and we might expect differences in the relatedness of these individuals, given that T. delaisi larvae are able to travel 11.5 km in 10 days[20]. However, the average relatedness values of early juveniles (−0.00099 ± 0.015) and late juveniles (−0.00027 ± 0.015) were not significantly different (t = −0.55, p-value = 0.58) (Fig. S5). Of the 157 juvenile individuals categorized as ‘locals’ due to high relatedness ratios, 29 were early juveniles and 17 late juveniles, whereas of the 157 juvenile individuals considered ‘recent migrants’, with low relatedness ratios, 19 were early and 17 late juveniles. Interestingly, one of the individuals considered a ‘local’ is the sole sample from la Pilona (Fig. 1) in the southernmost part of the large-scale sampling area.

All collected fish were geo-localized allowing the analysis of spatial structuring. We did not find differences in spatial distribution patterns between adults and juveniles representing two different years of recruitment (Kolmogorov-Smirnov D = 0.438, p = 0.09) when only space was considered. Immigration of less-related individuals, possibly from more distant populations, into the small-scale area could potentially result in different spatial patterns for these ‘recent migrants’. Adult migrants were differently distributed when compared to the rest of adults (D = 0.72, p < 0.001) with some differences evident in the southernmost part of the distribution (Fig. S6). However, ‘recent-migrant’ juveniles were not distributed differently from the rest of the juveniles (D = 0.25, p = 0.6994), nor were spatial differences found between early and late ‘recent-migrant’ juveniles in the small-scale area (D = 0.125, p = 0.996; Fig. S7). Furthermore, no spatial change was detected between recently recruited individuals (recent settlers) and individuals that survived a period after settlement (survivors) (D = 0.3125, p = 0.42; Fig. S8).

Fine-scale genetic and spatial clustering analysis

With hundreds of pairwise relatedness comparisons per individual, it is difficult to detect fine-scale genetic structure, although there could be spatial clusters. To evaluate genetic clustering in space, we first found significant spatial clusters in the small-scale area for both adults and juveniles (Fig. 3). We found six spatial clusters for the adults and seven in the juvenile samples. The additional cluster found for juveniles (SplClust6; Fig. 3b) is in an area that could correspond to non-suitable habitat for T. delaisi, as we found no adults and only a few juveniles.

Figure 3

Spatial clusters within the small-scale area for (a) adults and (b) juveniles. Spatial cluster 6 only exists for the juveniles and few individuals were found. The map was created using ArcGIS 10.5 (www.arcgis.com).

The proportion of individuals in each cluster that came from the SHORT-range disperser category varied from southwest to northeast along the coastline (Table S1). For juveniles, the frequency of SHORT-range dispersers increased along the coast from the southwest to the northeast, where they constituted 81% of all fish sampled, whereas the pattern was reversed for adults; that is, the frequency of SHORT-range dispersers decreased from southwest to northeast (Fig. 4). In contrast, LONG-range dispersers, which could be considered immigrants, were less frequent than SHORT-range dispersers at both life stages/years of recruitment (adults likely recruited in 2009 and juveniles in 2010) and were relatively homogeneously distributed along the coast (Fig. 4).

Figure 4

Proportion of adults and juveniles in each spatial cluster (SplClust; Fig. 3) in the SHORT-range disperser category, which are individuals with high mean relatedness with respect to their own spatial cluster (rIN) and low mean relatedness with respect to other spatial clusters (rOUT) as well as LONG-range dispersers, which have low relatedness to all clusters. SplClust6 is not included because it was only found in juveniles and in very low frequency.

Juveniles born at the beginning of the season (early juveniles) and at the end of the season (late juveniles) have differences in pelagic larval duration (Fig. S2), and different patterns were also observed for specific spatial clusters (Table 1). Spatial Cluster 4 (Fig. 3) is a rocky reef patch called ‘La Roca’ that is separated from the coastline by sand and has the furthest distance from shore. In Spatial Cluster 4, at the beginning of the season (with longer PLDs), only LONG-range disperser juveniles were found (Fig. 5). These individuals could be immigrants dispersing from distant populations. Late juveniles, with shorter PLDs, on the other hand, had a higher proportion of SHORT-range dispersers compared to LONG-range dispersers at La Roca (SplClust4) indicating a more local origin. Most LONG-range dispersers in late juveniles were found in spatial clusters 2 and 5 (Fig. 5).

Table 1

Proportion of individuals that belong to each spatial cluster in the large-scale sampling area.

Juveniles large scale area	SplClust1	SplClust 2	SplClust 3	SplClust 4	SplClust 5
N	98	119	174	79	157
MIXED-range dispersers	0	0.62	0.3	0.07	0.03
SHORT-range dispersers	0.7	0.34	0.55	0.75	0.1
MEDIUM-range dispersers	0.3	0.04	0.16	0.08	0.84
LONG-range dispersers	0	0	0	0.1	0.03

The clusters correspond to Fig. 6.

Figure 5

Proportion of low dispersers (a) which are individuals with a high relatedness with respect to their own spatial cluster (rIN) and low mean relatedness with respect to other spatial clusters (rOUT) and high dispersers (b) which are individuals with low mean relatedness with respect to all clusters, in early and late juveniles in each spatial cluster (from Fig. 3).

Figure 6

Spatial clusters for juveniles within the whole large-scale area. The map was created using ArcGIS 10.5 (www.arcgis.com).

Large-scale genetic and spatial clustering analysis

Genetic spatial clustering for juveniles across the whole large-scale sampling area, spanning a distance of almost 50 km, identified five distinct clusters, four of them including individuals within the small-scale area (Fig. 6). Some of the clusters of the large-scale analysis combined clusters identified in the small-scale area. For instance, SplClust 2 of the large-scale analysis (Fig. 6) included all individuals of SplClust 2, 3 and 4 in the small-scale analysis (Fig. 3).

The proportion of LONG-range dispersers, and thus likely immigrants, was very low for the large-scale dispersal area (Table 1). As the large-scale area is almost 50 km of coastline, these individuals would have to come from even more distant locations. All clusters had a higher proportion of SHORT-range dispersers than of MEDIUM or LONG-range dispersers, except Spatial Cluster 5, which is the northeastern most sampled area, which had a higher proportion of MEDIUM-range dispersers (Table 1), meaning that individuals in this cluster had less genetic relatedness to each other and were more related to the other locations from the southwest.

MIXED-range dispersers: individuals with high mean relatedness (r) with respect to individuals in their own spatial cluster (rIN) and high r with respect to individuals in other spatial clusters (rOUT). SHORT-range dispersers: individuals with high rIN and low rOUT. MEDIUM-range dispersers: individuals with low rIN and high rOUT. LONG-range dispersers: individuals with low rIN and low rOUT. Values above 0.5 are in bold.

Discussion

Dispersal and recruitment patterns in the black-faced blenny population

In the present work, we document very fine-scale patterns of dispersal and recruitment in a population of the black-faced blenny (Tripterygion delaisi) on the northern Mediterranean coast of Spain, through the combined use of genetic relatedness indices and a novel approach of binary data clustering that incorporates spatial patterns. In such a well-connected population with large numbers of dispersing individuals, genetic homogeneity would be expected and no population structure or population-wide pairwise relatedness was found with traditional methods[20]. However, binary clustering of genetic relatedness indices mapped in space revealed a slight clustering of T. delaisi individuals with different dispersal patterns. More adult blennies (sampled in 2010 and recruited mostly in 2009) considered short distance dispersers, were found in the southwestern parts of the small-scale study area. In contrast, juvenile individuals that recruited and were sampled in the following year (2010), had more short distance dispersers in the northeastern part of the 2 km coastline. Previous population genetic studies of this species with microsatellite loci, in the same small-scale area and the surrounding areas (up to 50 km), identified only one genetic unit[4,38]. Sibship analysis of T. delaisi from the same location found that full siblings can recruit many km apart from each other, further indicating a genetically homogenous population[20]. Nonetheless, the binary clustering analysis across the large-scale area demonstrated that a large proportion of juveniles collected northeast of the small-scale sampling region were MEDIUM-range dispersers and thus more related to the juveniles of this sampling region. This could potentially be the result of a previously unknown directional movement of larvae from southwest to northeast. Although high site-fidelity is observed in this species[40] we cannot completely rule out selective movement behaviour of adults to the southwestern part after having recruited the year before. However, there was no selective movement between recent settlers and survivors, individuals that survived a period after settlement, within one recruitment year. The detection of fine-scale clusters at both small and large-scale study areas shows that there is a complex but not full genetic mixing and suggests oceanic/coastal circulation as the main driver of the fine-scale chaotic genetic patchiness within this otherwise homogeneous population.

Larval dispersal by oceanic/coastal current drives fine-scale chaotic patchiness

Chaotic genetic patchiness, a term coined decades ago[52], has become commonly used to describe patterns of genetic diversity and population structure[27,28,53]. Broquet et al.[54] concluded through simulations that chaotic genetic patchiness could occur by a combination of genetic drift due to small effective size and collective dispersal at the larval stage. Such collective dispersal would most likely occur through retention of local larvae near the natal site and close to shore, as suggested by hydrodynamic modelling of larval dispersal[55]. Even though chaotic genetic patchiness is generally associated with marked genetic distances between cohorts, our results seem to be consistent with the dispersal hypothesis as the main driver of our fine-scale genetic patchiness. In essence, the higher degree of relatedness in the southwestern portion of the study area could be a result of a physical retention process of ‘locals’ at settlement or a collective movement after settlement. Importantly, these aggregations in space varied across two recruitment seasons, an additional sign of dispersal-driven fine-scale chaotic genetic patchiness.

The reproductive biology of Tripterygion delaisi could suggest that other mechanisms may play a minor role in generating fine-scale patterns of chaotic genetic patchiness. For this species, the number of adults (i.e., potentially reproductively active) is very large, with more than 1000 individuals on a narrow 2 km stretch of coastline[20]. Reproductively active individuals include females and two types of males, the dominant nest-making male and the sneaker male[40]. So-called sweepstake reproduction has been shown to cause genetic patchiness in other species[27,56] and for T.delaisi it could be that a few successful (i.e., the dominant) males contribute to the majority of offspring. However, the testes of T. delaisi sneaker males are significantly larger than those of dominant males and gene expression analysis found genes involved in differentiation and development, suggesting that these males are reproductively active[57]. A decline in diversity would be expected if only a proportion of the population is contributing to the next generation[58]. However, in our analyses, we did not find any evidence of sweepstakes reproductive success as genetic diversity did not decrease between adults and juveniles. Furthermore, neither heterozygosity nor relatedness differed for juveniles recruiting at the beginning or at the end of the breeding season, or between recent settlers and survivors. Sweepstakes-driven patchiness thus seems unlikely, even though high variability in the number of contributing adults over time has been shown in other marine fishes[56,59].

Patchiness could also be due to recruitment differences related to the conditions encountered by larvae. Pelagic larval duration (PLD) varied significantly with the hatching date of the larvae, with PLD up to 10 days shorter at the end of the season than at the beginning. Larvae that spend more time in the water column have more opportunities to disperse further and might include more immigrants from other populations. Yet, in the whole sampling area a similar proportion of T. delaisi migrants arrived early and late in the season. Nonetheless, at ‘La Roca’ (Spatial Cluster 4), a submerged rocky habitat not continuous with the rocky shoreline habitat, all early juveniles collected (N = 29) were LONG-range dispersers and can be considered migrants. This could hint at collective dispersal, at least at the beginning of the season, which Broquet and coauthors[54] suggested as the second source for chaotic genetic patchiness and greater than expected levels of kinship as an indicator have been found in other marine organisms[60]. However, a previous study only found one full sibling pair of T. delaisi settling together in the same area[20] and the fine-scale patchiness pattern might be due largely to seasonal variation in current patterns.

One factor influencing the PLD in T.delaisi was the sea temperature. In our study area, the difference in sea surface temperature between the beginning (April) and the end of the reproductive period (July) is approximately 8 °C. This increase in temperature resulted in a decrease of PLD. Consequently, increased ocean temperatures, due to global climate change, could decrease connectivity between populations and increase self-recruitment due to a general decline in PLD at higher temperature[13,61-63]. The effect of temperature on larval connectivity has also been modelled and results suggested a considerable change in larval recruitment and connectivity[64]. Such changes in dispersal potential and connectivity in T. delaisi could have detrimental effects on population genetic diversity. Nonetheless, the seasonal increase is larger than the predicted temperature increase in the western Mediterranean Sea (ca. 0.25 °C per decade) towards the end of the present century due to climate change[65]. Mean PLD of the black-faced blenny clearly decreases with increased temperature, but in the small-scale area studied in the present work this had little effect on the overall connectivity of the population since the proportion of individuals arriving from distant localities was similar regardless of PLD or water temperature. Even with a possible reduction in PLD over time due to climate change, intra-annual variation in conditions encountered by larval T. delaisi may provide a buffer and allow maintenance of generally well-connected populations.

Conclusion

The species Tripterygion delaisi on the Mediterranean coast of Spain is a well-connected unit but appears to have fine-scale structure in both space and time. The combination of geo-localized individual locations, hatching date information, and genotypic SNP data allowed detection of patches of larvae with particular dispersal histories. Yet, very fine-scale genetic patchiness across space and time was only apparent when using binary clustering techniques, previously used primarily for animal movement trajectory segmentation[37]. The ability to detect such fine-scale dispersal patchiness will aid in our understanding of the underlying mechanisms of population structuring and chaotic patchiness in a wide range of species, including those with high potential dispersal abilities.

Supplementary Figures & Table