Phylogeographic structure and population connectivity of a small benthic fish (Tripterygion tripteronotum) in the Adriatic Sea

Population connectivity of benthic marine organisms depends strongly on planktonic larval dispersal and is controlled by geographic distance and oceanographic structure. We examine isolation by distance versus resistance to barriers (ocean current boundaries) against a background of post‐glacial habitat expansion in a small benthic fish of the Adriatic Sea.


| INTRODUC TI ON
The apparent continuity of the marine environment suggests virtually unrestrained dispersal and large-scale connectivity within species, but this superficial view is contradicted by genetic studies (Hellberg, 2009;Palumbi, 2003). Although many large-bodied marine fish species are highly mobile, their generational dispersal seems to be restricted and does not preclude the evolution of strongly structured populations (e.g. Castro et al., 2007;Graves, 1998). Similarly, planktonic larval dispersal -the main if not only manner of dispersal for sessile and small benthic species -is not at all boundless. Numerous factors, including egg type, pelagic larval duration and swimming behaviour, have been found to influence the extent of planktonic larval dispersal (Nanninga & Manica, 2018;Selkoe & Toonen, 2011). Moreover, rather than following a simple pattern of isolation by distance, planktonic dispersers often show 'isolation by resistance', with population connectivity affected by oceanographic features (Abrahamson, Eubanks, Blair, & Whipple, 2001;Banks et al., 2007;Benestan et al., 2016;Galarza et al., 2009;Gilg & Hilbish, 2003;Riginos, Douglas, Jin, Shanahan, & Treml, 2011;Schunter et al., 2011;Thomas et al., 2015;Truelove et al., 2017;White et al., 2010;Xuereb et al., 2018). As a consequence, adjacent populations can become differentiated across short geographic distances while exchanging genes with more distant populations (Gilg & Hilbish, 2003;Xuereb et al., 2018). By identifying spatial patterns of connectivity and isolation, genetic data provide valuable information for the design and improvement of conservation and management strategies (Giakoumi et al., 2013;Jahnke et al., 2017;Rilov et al., 2019).
The Adriatic Sea is a biogeographically distinct part of the Mediterranean Sea (Bianchi, 2007;Spalding et al., 2007;examples: Koblmüller, Steinwender, Weiß, & Sefc, 2015;Maltagliati, Di Giuseppe, Barbieri, Castelli, & Dini, 2010;Souche et al., 2015;Wagner et al., 2019), to which it is connected by the narrow Otranto Strait. The 800-km long body of water is enclosed by a regular, sandy and gently sloping coast to the west and a steep, rocky coast along with numerous coastal islands to the east. It is divided into three sub-basins -the northern, central and southern Adriatic (Russo & Artegiani, 1996; Figure 1), which also correspond to recognized biogeographic sectors (Bianchi et al., 2012) with variation in depth, temperature and the distribution of boreal and thermophilic fauna (Bianchi & Morri, 2000;Jardas, 1996). In addition to the basin-wide surface circulation, which runs along the eastern and western Adriatic coast (the East and West Adriatic currents, respectively), recirculation gyres form in each of the three sub-basins. The northern transversal edge of the Southern Adriatic gyre approximately matches the boundary between the central and southern Adriatic sub-basins, the northern transversal edge of the Central Adriatic gyre approximately matches the boundary between the northern and central Adriatic sub-basins and a fourth cyclonic gyre exists in the very northern part of the Adriatic between the Istrian peninsula and the Po river delta (Poulain, 2001; Figure 1). Surface circulation is influenced by winds and is seasonally variable (Poulain, 2001;Russo & Artegiani, 1996;Ursella, Poulain, & Signell, 2007). Consequently, passive transport in the water currents underlies seasonal variation and population connectivity in species with planktonic larval dispersal will to some degree depend on the time of spawning (Melià et al., 2016). In spring and summer, for instance, the East Adriatic Current, which runs northwards along the Eastern coast, is weak, whereas the recirculation cells in the northern and central sub-basins as well as the cyclonic gyre along Istria are well developed (Poulain, 2001;Ursella et al., 2007).
Apart from current seascape structure, a major impact on population structure of aquatic biota in the Adriatic was produced by sea level changes during the glacial cycles. During the Last Glacial Maximum (24,000-18,000 years BP), the sea level was more than 120 m lower than today (Lambeck, Rouby, Purcell, Sun, & Sambridge, 2014), leaving the northern Adriatic desiccated. With rising sea levels, this area became available for re-colonization by marine organisms and the present coastline was established about 7,000 years BP (Fleming et al., 1998;Lambeck, Antonioli, Purcell, & Silenzi, 2004). Population structure in the northern Adriatic is therefore of recent origin and possibly affected by bottlenecks, founder events and population expansion (Koblmüller et al., 2015;Ledoux et al., 2018;Maltagliati et al., 2010).
Despite the surging interest in seascape genetics in general and the oceanographic appeal of the Adriatic Sea, relatively few studies have addressed population structure in the Adriatic.
Among those that did, some found that genetic clusters within species reflected the Adriatic sub-basins (Carreras et al., 2017;Ledoux et al., 2018;Schiavina, Marino, Zane, & Melià, 2014), whereas in others, there was only weak or no differentiation by between realized and potential connectivity illustrate the value of integrating different data sources to understand population structure and inform conservation planning.
In addition to the division into large sub-basins, the ragged shoreline and wealth of near-shore islands along the eastern Adriatic coast provide a richly structured environment, which could further shape population differentiation of littoral and benthic species on a more local scale (Ledoux et al., 2018). The analysis of patterns of genetic differentiation within the Adriatic sub-basins requires a much denser sampling regime than was performed in most previous studies.
In the present study, we performed dense sampling of the triplefin blenny Tripterygion tripteronotum (Tripterygiidae; Figure 1) along the eastern Adriatic shoreline and investigated population structure in relation to the Adriatic sub-basins and recirculation gyres. The small (<7 cm) fish are restricted to the shallow (0-6 m) rocky littoral (de Jonge & Videler, 1989;Wirtz, 1978).
Their distribution therefore follows a narrow band along the rocky coast of the eastern Adriatic and is interrupted between islands and along stretches of sandy coast, which rarely exceed a few hundred metres in length along this part of the Adriatic.
As adult triplefins are non-migratory, highly territorial, epibenthic to cryptobenthic fish with relatively sedentary lifestyle (Schunter, Pascual, Garza, Raventós, & Macpherson, 2014;Wirtz, 1978; own observations), dispersal over any significant distance would necessarily occur during the larval phase. Tripterygion tripteronotum attach their eggs to algal aufwuchs in the males' territories (de Jonge & Videler, 1989). After hatching, larvae disperse into the inshore waters for 2-3 weeks (Macpherson & Zika, 1999;Raventos & Macpherson, 2001), but apparently remain close to the shore as no larvae were found in hauls >2.5 km from the shore (Sabates, Zabala, & Garcia-Rubies, 2003). Spawning occurs from early spring to late summer (de Jonge & Videler, 1989;Macpherson & Zika, 1999;Wirtz, 1978), and accordingly, Tripterygion larvae were abundant in plankton samples collected at depths of 2-3 m near the Spanish coast from April to September (Sabates et al., 2003). In the related species T. delaisi, which has a similar planktonic larval duration and distribution as T. tripteronotum, the estimated percentage of self-recruitment to a study population in north-eastern Spain was as low as 6.5% (Schunter et al., 2014), which, however, did not prevent population differentiation to evolve on a geographic scale of >100 km in the same area (Carreras et al., 2017;Carreras-Carbonell, Macpherson, & Pascual, 2006, 2007. We hypothesized that planktonic dispersal of T. tripteronotum larvae may likewise be restricted, and its scale and direction determined by larval behaviour, oceanic currents and the topography of the eastern Adriatic coast. Differentiating between effects of geographic distance and effects of potential dispersal barriers called for a dense sampling of locations both within and across recirculation gyres. To assess the role of passive larval dispersal following the Adriatic currents in shaping population connectivity, we compared genetic population structure with potential connectivity estimated from simulations of particle drift.

| Sampling and sequencing of the mitochondrial control region
Tripterygion tripteronotum were sampled from 25 locations along the Croatian coast of the eastern Adriatic Sea (n = 550; Figure 1), and from four Eastern Mediterranean locations (n = 39). Table 1  Fish were captured at depths of 0-2 m with hand nets while snorkelling. Fin clips were taken from the caudal fins and preserved in 96% ethanol. Fish were released immediately and were frequently observed to return to their territories. Whole genomic DNA was extracted following a rapid Chelex protocol (Richlen & Barber, 2005), and a part of the mitochondrial control region was PCR amplified with primers L-Pro-F_Tropheus (Koblmüller et al., 2011) and TDK-D (Lee, Conroy, Huntting Howell, & Kocher, 1995) as described previously (Koblmüller et al., 2011(Koblmüller et al., , 2015. We expected the mtDNA locus to be a suitable marker for population structure, as mito-

| Geographic patterns of differentiation and diversity
We conducted distance-based redundancy analysis (dbRDA; Legendre and Anderson (1999)) using the R package 'vegan' (Oksanen et al., 2015) to separate the effects of ocean circulation and geographic distance on genetic structure. Two separate analyses were performed using the function 'capscale' with either pairwise θ ST or φ ST values between the Adriatic population samples as dependent variables. In each of the two analyses, we included a categorical predictor variable representing assignment to one of the four recirculation gyres and a continuous predictor variable representing geographic distance between locations. The assignment of sampling locations to recirculation gyres is shown in Table 1 and was based on maps of surface flow in the Adriatic basin (Poulain (2001); for a schematic representation, see Figure 1) and Lagrangian simulations of particle drift (see below). Geographic distances between sampling locations were calculated from their longitudinal and latitudinal co-  (Hedrick, 2005). Significance of predictor effects was assessed with 9,999 permutations using the function 'ANOVA.cca'.
We used non-metric multidimensional scaling of the popula-

| Coalescent modelling of divergence, migration and demographic history
To estimate rates of gene flow between the Istrian and Northern Adriatic gyre, we applied the isolation-with-migration coalescent model implemented in IMa2 (Hey, 2010). The model employed the HKY model of sequence evolution (Hasegawa, Kishino, & Yano, 1985 (Koblmüller et al., 2015), a species closely related to T. tripteronotum (Carreras-Carbonell, Macpherson, & Pascual, 2005).
Past population size changes were inferred separately for the two phylogeographic groups (Istrian and Northern Adriatic gyre; Central and Southern Adriatic gyre) by means of a Bayesian coalescent approach (Bayesian skyline tree prior) as implemented in BEAST 1.10.4 (Suchard et al., 2018). We employed the same model of molecular evolution as for ML-tree inference (see above) and assumed a strict molecular clock and alternative substitution rates of 1.11%, 2.61%, 3.61% and 6.77% per site per MY (Koblmüller et al., 2015). For both datasets, two independent runs of 100 million generations each were conducted, sampling every 1,000th step, with a burn-in of the first 10% of sampled generations. Tracer 1.5 (Rambaut & Drummond, 2009) was used to check for run convergence and sufficiently large effective sample sizes (ESS > 200; Kuhner, 2009). Runs were combined using LogCombiner (part of the BEAST package) and Bayesian skyline plots (BSPs) were visualized in Tracer.

| Lagrangian modelling of passive dispersal
Lagrangian simulations were performed using the Parcels tool- To assess potential patterns of passive larval dispersal in the study area and connectivity between the sampled populations, Lagrangian particles were advected in the AREG surface flow fields.
Particles were released at the study locations (see Table 1

| A phylogenetic break at the boundary between the central and the northern sub-basin
The ML-tree-based haplotype network (

| Population structure is associated with the recirculation gyres
We conducted distance-based redundancy analysis (dbRDA) to evaluate isolation by gyres and isolation by geographic distance, and found a significant association of genetic structure with gyres We used hierarchical AMOVAs to partition genetic variance between and within gyres. Accounting for genetic distances between haplotypes (Φ ST analysis), the largest proportion of variance was assigned between gyres (57% of total variance; Φ CT = 0.57, p < 0.001), with much less variance among populations within gyres (1.93% of total variance; Φ SC = 0.04, p < 0.001). Hierarchical AMOVA based on haplotype frequencies (θ ST analysis) likewise assigned a larger proportion of variance to between-gyre structure (21% of total variance; θ CT = 0.21, p < 0.001) than within-gyre population structure (3.75% of total variance; θ SC = 0.05, p < 0.001).

| Population differentiation within recirculation gyres
As shown above, a small but significant proportion of genetic variance was assigned to within-gyre structure by hierarchical ANOVA. Within the northern (11 populations) and the central

| Temporal stability
As Adriatic populations were sampled from 2006 to 2017 (Table 1) The analogous test for the phylogenetic clade comprising the central and southern sub-basins is not meaningful because sampling years are correlated with geographic regions (Table 1). Overall, there is no evidence in our data that the protracted sampling period might have an effect on our analyses.

| Lagrangian simulations of passive dispersal
We used Lagrangian simulations of particle drift to examine whether genetic differentiation across recirculation gyres could be explained by retention of passively drifting larvae within each of the gyres ( Figure 5, Appendix S4). Simulated particle drift reflected the recirculation gyres, but did not adhere strictly to gyre boundaries.
Northward drift between gyres was observed at a higher rate than southward drift, and as expected, drift patterns showed seasonal variation across the simulated period (February to September; Appendix S4).

| No northwards decline in genetic diversity
Estimates of genetic diversity within each population are shown in Table 1. Four to 18 different haplotypes (median value = 7) were detected in the Adriatic populations, and haplotype diversity (H E ) ranged from 0.38 to 0.93 (mean ± SD = 0.69 ± 0.15). Haplotype diversity was independent of sample size, but the number of haplotypes was positively correlated with sample size. Therefore, we estimated haplotype richness (HR) based on the minimum sample size of n = 13 (range: 3.58-8.83 haplotypes, mean ± SD = 5.75 ± 1.41) to allow comparisons among populations. Contrary to expectations for a step-wise northwards range expansion, neither haplotype richness nor haplotype diversity co-varied with latitude (

| Recent population expansion in the Adriatic Sea
The haplotype network of the Adriatic samples comprises a small number of highly frequent haplotypes and numerous low-frequency haplotypes, a pattern that is indicative of recent population expansion ( Figure 2). Coalescence models (Bayesian skyline plots) supported recent population expansion in both phylogeographic groups ( Figure 6). The onset of population growth in the northern sub-basin clade was dated to >10,000 to <75,000 years ago, dependent on the F I G U R E 5 Particle dispersal in Lagrangian drift simulations. Each dot represents the arrival location of a particle after 20 days of passive drift. Particles were released at the T. tripteronotum sampling locations within each of the four Adriatic gyres every 12 hr in the period from February to September throughout the years 2004-2018 (see Table 1) [Colour figure can be viewed at wileyonlinelibrary. com] assumed substitution rate. Consistent with the hydrologic history, population growth reconstructed for the central/southern clade was less pronounced and preceded the northward expansion.

| Relationship to Mediterranean populations
The ML-tree-based network (Figure 2)

| Population structure in relation to the Adriatic circulation patterns
The present study employed dense geographic sampling along the eastern Adriatic coast to differentiate between the effects of distance versus oceanic circulation patterns on marine population connectivity. We found that genetic population differentiation of T.
tripteronotum corresponded to boundaries between the Adriatic recirculation gyres. This result blends in with findings in other marine species and regions, where circulation systems control gene flow and generate patterns of 'isolation by resistance' (Banks et al., 2007;Benestan et al., 2016;Thomas et al., 2015;White et al., 2010;Xuereb et al., 2018). Hence, the observed population structure in the northern Adriatic sub-basin can be explained by recent divergence in the wake of the filling of the northern basin after the last glacial maximum (Lambeck et al., 2014), and incomplete isolation between the two gyres with asymmetric gene flow mediated by northward drifting larvae.
In contrast, drift of particles from the Central Adriatic into the Northern Adriatic gyre predicted that gene flow between these two gyres, if driven by passive drift, should occur at higher rates than indicated by the observed phylogeographic break. Our simulation results concur with Lagrangian analyses by Rossi, Ser-Giacomi, Lopez, and Hernandez-Garcia (2014), who likewise identified a leaking hydrodynamical province boundary in the same region. As in our study, incongruity between potential connectivity inferred from drift simulations and the connectivity realized by presumably passive dispersal was observed in the seagrass Posidonia oceanica across the Adriatic basin (Jahnke et al., 2017). Suggested explanations for this discrepancy included clonal reproduction, low settlement success and small-scale hydrodynamics (Jahnke et al., 2017). In T. tripteronotum, mechanisms such as larval swimming behaviour or differential settlement success (Riginos et al., 2011;Selkoe & Toonen, 2011;Shanks, 2009) may contribute to the restriction of gene flow between the Central and the Northern Adriatic gyre, but it is presently unclear which oceanographic features confine the dispersal of T. tripteronotum larvae in that region.
It is noteworthy that the phylogeographic break coincides not only with the sub-basin and gyre boundaries but also with a gap in the conglomerate of nearshore islands and with a stretch of unusually steeply sloping shoreline (Figure 1). Possibly, if larvae are retained close to their F I G U R E 6 Bayesian skyline plots depicting estimated changes in population size of Adriatic T. tripteronotum through time. Thick lines and broken lines represent the median and 95% HPD intervals respectively. The y-axis represents the population size parameter (female effective population size times the substitution rate). The x-axes are scaled by alternative substitution rates to illustrate the range of time estimates natal sites (Brandl et al., 2019), genes may move in a wide front along the island-rich coastline, whereas gene flow may be curbed along the narrow belt of suitable habitat in the region where we observe the phylogeographic break. Furthermore, selection could contribute to the observed break either in that fish achieve very little reproductive success after crossing between gyres or by selecting against mitochondrial introgression in the face of nuclear gene flow (Morales, Pavlova, Joseph, & Sunnucks, 2015). In both cases, an abrupt change in the selection regime between the central and the northern sub-basin would be required to counter the gene flow potential that was suggested by the drift simulations.
Based on available data, it is difficult to assess whether the phylogeographic break between approx. 43.6° and 44° northern latitude on the eastern Adriatic coast is generally shared across species.
In the densely sampled octocoral Paramuricea clavata, a distinct genetic discontinuity was detected between the Kvarner area (represented in our sample by locations 5-12, Figure 1) and the Kornati area (Ledoux et al., 2018). The Kornati sites of Ledoux et al. (2018) are located at ~43.8° northern latitude, i.e. within the latitudinal range of the phylogeographic break in T. tripteronotum, such that the genetic divergences in T. tripteronotum and P. clavata may indeed coincide geographically. Other studies in the Adriatic Sea lack the necessary geographic sampling to locate genetic breaks (Carreras et al., 2017;Jahnke et al., 2017;Koblmüller et al., 2015;López-Márquez et al., 2019;Schiavina et al., 2014).
In contrast to the situation in T. tripteronotum, no genetic divergence between populations from different gyres seems to exist in another triplefin species, T. delaisi, as specimens sampled from one location in the Central Adriatic gyre grouped among populations of the Northern Adriatic gyre (Koblmüller et al., 2015). Adult Tripterygion delaisi occur in the same regions of the Eastern Adriatic as adult T.
tripteronotum, but in a much wider range of depth (typically 3-40 m (Wirtz, 1978), except for occasional, more shallow occurrence at very shadowy places (Koblmüller et al., 2015)). In New Zealand, triplefin species inhabiting the very shallow benthic zone are also more strongly structured than deeper-dwelling species (Hickey, Lavery, Hannan, Baker, & Clements, 2009). Hickey et al. (2009) concluded that 'forces structuring populations, whatever these may be, decline with increasing depth'. It will be useful to complement the present study with investigations of more species, for instance, representatives of blenniid fishes, in order to identify determinants of (crypto) benthic population structure along the Adriatic coast.

| Connectivity among populations within gyres
With benthic, territorial adults and inshore-dispersing planktonic larvae (de Jonge & Videler, 1989;Sabates et al., 2003;Wirtz, 1978), gene flow among T. tripteronotum populations might be limited by distance, as observed, for instance, in the octocoral Paramuricea clavata in the same region (Ledoux et al., 2018) Macpherson, & Pascual, 2006). In contrast to this prediction, population structure within gyres was low and we detected no effect of geographic distance on genetic population structure after accounting for structure across gyre boundaries.
We also observed numerous disagreements between θ ST -and Φ ST -based inferences of population differentiation within gyres, which are likely due to the different sensitivities of the two statistics to drift and mutations at low levels of differentiation (Sefc, Payne, & Sorenson, 2007). Both test statistics indicate weak differentiation of populations situated in an enclosed bay on the island of Pag (locations 11 and 12) and in the Split Channel (locations 17-19) from several other populations within the respective recirculation gyre. We conclude that connectivity within gyres is generally high, but may be reduced in the case of populations at secluded locations.

| Post-glacial recolonization of the northern Adriatic and relationship between Adriatic and Eastern Mediterranean haplotypes
Genetic signatures of population expansion were consistent with a rapid recolonization of the northern sub-basin. A comparable signal was previously observed in the related triplefin T. delaisi (Koblmüller et al., 2015), whereas no evidence for rapid population expansion in the northern Adriatic was detected, for instance, in the octocoral P.
clavata (Ledoux et al., 2018). In contrast, populations of P. clavata displayed the expected diversity gradient along the recolonization axis (Ledoux et al., 2018), which was absent in the present study. In

| CON CLUS IONS
The dense sampling employed in the present study allowed us to separate correlations of population structure with geographic distance (isolation by distance) from correlations with oceanographic features (isolation by resistance) in the Adriatic Sea. While we document a strong relationship between population structure and the Adriatic sub-basin and gyre structure, we also demonstrate that population differentiation across one of the investigated gyre boundaries is stronger than predicted by the simulated effect of gyres on passive larval dispersal. Our data clearly support a role for ocean currents in shaping population structure and connectivity in our study species, but additional mechanisms seem to be at work to reduce gene flow between the central and the northern Adriatic sub-basin.
Biophysical modelling of potential connectivity has been recognized as a valuable tool for marine conservation planning (Coleman et al., 2017;Melià et al., 2016;Rossi et al., 2014). Our results corroborate the insight that species-specific differences in realized connectivity (Paterno et al., 2017) and discrepancies between realized and potential connectivity (Jahnke et al., 2017) call for integrating oceanographic simulations with biological data collection, preferably across a range of species, in order to inform efficient conservation planning.

ACK N OWLED G EM ENTS
We thank Simon Brandl for insightful discussion and comments on the manuscript, and Holger Zimmermann for assistance with sam-