Gene flow, population growth and a novel substitution rate estimate in a subtidal rock specialist, the black‐faced blenny Tripterygion delaisi (Perciformes, Blennioidei, Tripterygiidae) from the Adriatic Sea

Abstract Population histories depend on the interplay between exogeneous and endogeneous factors. In marine species, phylogeographic and demographic patterns are often shaped by sea level fluctuations, water currents and dispersal ability. Using mitochondrial control region sequences (n = 120), we infer phylogeographic structure and historic population size changes of a common littoral fish species, the black‐faced blenny Tripterygion delaisi (Perciformes, Blennioidei, Tripterygiidae) from the north‐eastern Adriatic Sea. We find that Adriatic T. delaisi are differentiated from conspecific populations in the remaining Mediterranean, but display little phylogeographic structure within the Adriatic basin. The pattern is consistent with passive dispersal of planktonic larvae along cyclonic currents within the Adriatic Sea, but limited active dispersal of adults. Demographic reconstructions are consistent with recent population expansion, probably triggered by rising sea levels after the last glacial maximum (LGM). Placing the onset of population growth between the LGM and the warming of surface waters (18 000–13 000 years BP) and employing a novel expansion dating approach, we inferred a substitution rate of 2.61–3.61% per site per MY. Our study is one of only few existing investigations of the genetic structure of animals within the Adriatic basin and is the first to provide an estimate for mitochondrial control region substitution rates in blennioid fishes.


Introduction
Environmental factors such as geology, hydrology and climate can interact strongly with ecological, morphological and reproductive characteristics of species in shaping their phylogeographic structure. For example, the glacial cycles in the Quaternary caused dramatic environmental changes that severely influenced the distribution, phylogeographic structure and demography of both aquatic and terrestrial taxa (Hewitt 2000). In marine habitats, the sea level fluctuations associated with glacial cycles had a strong impact on shelf regions and peripheral seas, and have particularly affected the demography and phylogeographic patterns of shallow water taxa, especially of those with low active dispersal capacity (e.g. Palumbi 1994;Hellberg 2009). Another important determinant of phylogeographic structure in the marine environment are water currents, which provide for dispersal during planktonic life stages (e.g. Cowen and Sponaugle 2009;Shanks 2009). The Adriatic Sea, a part of the Mediterranean Sea, is one of the marine regions, in which both sea level fluctuations and water currents may have shaped current population genetic and phylogeographic patterns to a great extent. It is connected to the rest of the Mediterranean Sea by the narrow Otranto Strait (Fig. 1) and features characteristic salinity, temperature, depth, circulation and productivity patterns (Astraldi et al. 1999). Circulation patterns in the Adriatic are dominated by cyclonic gyres, which may interrupt planktonic dispersal both within the Adriatic basin and across the Otranto Strait. Large parts of the Adriatic Sea were dry during the last glacial maximum (LGM;~18 ka BP) when the sea level was 120 m below its present level (Siddall et al. 2003). Subse-quently, rising sea levels shifted the coastline northwards and expanded the distribution ranges of Adriatic taxa considerably (Lambeck and Purcell 2005). This rapid range expansion should have left detectable traces of population growth, especially in taxa with no or little gene flow to and from the rest of the Mediterranean Sea (Debes et al. 2008;Luttikhuizen et al. 2008). The black-faced blenny Tripterygion delaisi Cadenat & Blache 1969 (Perciformes, Blennioidei, Tripterygiidae) is a common littoral fish species in the Mediterranean Sea and the eastern Atlantic from Senegal north to the British Isles, including the Azores and Canary islands. It inhabits rocky habitats at depths of 3-40 m (Wirtz 1978) and can be observed close to the water surface only at very shadowy places (S. Koblm€ uller and K. M. Sefc, personal observation). Adults have a reduced swim bladder and are not able to traverse long distances across open water, which may turn, for example, larger stretches of sandy habitat into serious dispersal barriers. This low active dispersal capacity, together with high levels of territoriality in males and a well-developed homing behaviour over short distances (Heymer 1977;Wirtz 1978), is a precondition for pronounced geographic population structure. However, in contrast to the rather stationary adult life stage, dispersal may occur during the planktonic larval phase that lasts for 16-21 days (Ravent os and Macpherson 2001). While the duration of the planktonic stage would allow for long-distance dispersal, high levels of self-recruitment indicate that a significant proportion of larvae remain near their natal area (Carreras-Carbonell et al. 2007), and indeed, larvae are predominantly found in coastal areas (Borges et al. 2007). Nonetheless, larval drift via water currents likely represents the predominant mode of dispersal and enables T. delaisi to colonize new habitat.
Tripterygion delaisi from the eastern Atlantic islands (Azores, Canaries) are genetically divergent from T. delaisi in the European Atlantic coast and the Mediterranean, and population differentiation exists between distant locations in the latter areas (Carreras-Carbonell et al. 2006;Domingues et al. 2007). At a smaller geographic scale, microsatellite markers revealed significant differentiation and a pattern of isolation by distance among populations along the Spanish Mediterranean coast (Carreras-Carbonell et al. 2006). Except for two individuals from an unspeci-fied location in northern Croatia (Domingues et al. 2007), Adriatic T. delaisi have never been examined. Here, we use mitochondrial control region sequences of T. delaisi from the northeastern Adriatic Sea (Croatia) to analyse the phylogeographic structure within this area and compare them to published sequences to investigate the relationship to T. delaisi from other parts of the Mediterranean and from the Atlantic Sea. The sampled area belongs to one of the cyclonic gyres (the middle Adriatic gyre, Fig. 1), within which water currents may provide for substantial dispersal of planktonic larvae. Given the life history of the species and the geomorphological and hydrological features of the Adriatic Sea, we predicted that the sampled populations would show genetic signatures of recent population expansion, only weak phylogeographic structure, but genetic differentiation from their conspecifics in the rest of Mediterranean. Furthermore, we apply an expansion dating approach to provide a first substitution rate estimate of the mitochondrial control region for the suborder Blennioidei.

Sample collection and DNA sequencing
Eighty-four individuals of T. delaisi were caught with hand nets at 10 localities in Croatia between 2006 and 2012 ( Fig. 1, Appendix). Fin clips (small parts of the caudal fin) were taken and preserved in 96% ethanol, and fish were immediately released. Whole genomic DNA was extracted following a rapid Chelex protocol (Richlen and Barber 2005). The most variable part of the mitochondrial control region was amplified and sequenced according to the protocols described in Koblm€ uller et al. (2011) andDuftner et al. (2005), respectively. The primers used for PCR and chain termination sequencing were L-Pro-F_Tropheus (Koblm€ uller et al. 2011) and TDK-D (Lee et al. 1995). DNA fragments were purified with Sephadex TM G-50 (GE Healthcare, Vienna, Austria) and visualized on an ABI 3130xl capillary sequencer (Applied Biosystems, Vienna, Austria). Additional sequences from Portugal, Italy, Cyprus, Croatia and the Azores and Canary islands (Domingues et al. 2007) were obtained from GenBank ( Fig. 1, Appendix). Sequences were aligned by eye (no indels were present in the alignment) in MEGA 6.06 (Tamura et al. 2013). The length of the final alignment was 352 bp. Sequences are deposited in GenBank under the accession numbers KT267998-KT268081.

Phylogeographic analyses
Phylogenetic relationships among haplotypes were inferred by means of a statistical parsimony network (Templeton et al. 1992) and a maximumlikelihood (ML) tree in TCS 1.2.1 (Clement et al. 2000) and PHYML 3.0 (Guindon et al. 2010), respectively. For ML tree search, identical sequences were collapsed into haplotypes using DNACollapser implemented in FABOX (Villesen 2007). The HKY+I (Hasegawa et al. 1985) model was employed as best fitting model of sequence evolution selected by the Bayesian information criterion (BIC) in JMODELTEST 0.1 (Posada 2008), and statistical support was assessed from 1000 bootstrap replicates.
Haplotype (H d ) and nucleotide diversity (p) were calculated in DNASP 5.10 (Librado and Rozas 2009) for the entire Mediterranean + Atlantic coast samples, as well as for the Croatian samples only. A genetic landscape shape interpolation analysis from the Croatian samples was performed in ALLELES IN SPACE 1.0 (AIS: Miller 2005). The genetic landscape Fig. 1. Sampling localities of Tripterygion delaisi used in the present study, a representative picture of a male individual T. delaisi and genetic landscape shape visualization. (a) Sampling locations in the Mediterranean. Sequences from Portugal, Italy, Cyprus as well as two sequences from Croatia (exact sampling locality unknown) and the out-group sequences from the Azores and Canary islands (not shown in the map) are from Domingues et al. (2007). (b) A typical territorial male T. delaisi. (c) Bathymetric map of the Adriatic Sea with major currents (left) and sampling sites along the Croatian coast and islands (right), overlaid by the genetic landscape as inferred from the genetic landscape shape interpolation analysis (the degree of shading reflects genetic distance) of Adriatic T. delaisi. EAC, eastern Adriatic current; WAC, western Adriatic current; DWOC, deepwater outflow current; NAG, north Adriatic gyre; MAG, middle Adriatic gyre; SAG, south Adriatic gyre. For island populations, the island names are given in parentheses. Numbers to the right of the location name indicate sample size.
analysis employed a Delaunay triangulation connectivity network and residual genetic distances derived from a linear regression of genetic distances on geographic distances (as recommended for data sets with substantial variation in geographic distances between sampling sites connected in the Delaunay triangulation network; Manni et al. 2004). Grid size was set to 0.05 9 0.05 latitude and longitude degree, respectively, and a distance weighting parameter a = 1 was used (see Miller et al. 2006 for more details on the method). Qualitatively similar results were obtained applying different grid sizes and a range of distance weighting parameters (a = 0.3-3; not shown).

Substitution rate estimation by expansion dating
To test for signals of past population expansion, we calculated Fu's F S (Fu 1997) and Tajima's D (Tajima 1989) in DNASP 5.10 and employed a Bayesian coalescent approach [Bayesian skyline plot (BSP)], implemented in BEAST 1.8.0 (Drummond and Rambaut 2007). For the BSP analysis, we applied the HKY+I model of molecular evolution and a strict molecular clock model. Two independent MCMC runs of ten million generations each were conducted, sampling every 1000th step with a burn-in of the first 25% of sampled generations. Verification of effective sample sizes (ESS > 200 for all parameters), trace of MCMC runs and visualization of demographic changes were done in TRACER 1.5 (available from http://beast.bio.ed.ac.uk/tracer.). The logged parameter values and trees from the two replicate runs were combined using LOGCOMBINER 1.8.0 (part of the BEAST package), and TRACER was used to create a BSP.
We are not aware of published substitution rate estimates for the mitochondrial control region in Tripterygiidae, or even the Blennioidei. Rate estimates for closely related perciform families range from 1.8 to >33% per MY (reviewed in Bowen et al. 2006), rendering the selection of a rate for Tripterygion into guesswork. Therefore, we used the data produced in this study to derive a substitution rate based on expansion dating. During the LGM (~18 ka BP), the sea level was~120 m below its present level (Siddall et al. 2003) and large parts of the Adriatic Sea were dry (Trincardi et al. 1994;Lambeck and Purcell 2005). The following deglaciation was associated with a global rise of the sea level, which started gradually and reached a maximum rate at~14.5 ka BP (Crandall et al. 2012). Shortly thereafter, at around 13 ka BP, surface water temperatures in the Adriatic increased to temperate conditions (Zonneveld 1996). The rising sea levels shifted the Adriatic coastline northwards by~250 km (Lambeck and Purcell 2005) and greatly expanded the available habitat of T. delaisi. Provided that the habitat expansion gave rise to population growth, the palaeohydrologic and palaeoclimatic data suggest that the onset of population growth should have occurred between 18 and 13 ka BP.
The BSP of the Adriatic T. delaisi sample indicated a sudden transition from an extended period of constant population size to rapid growth in the more recent past. Therefore, we adopted the expansion dating approach of Crandall et al. (2012) and calibrated the tripterygiid molecular clock based on the assumption that this transition coincided with either the rising sea level or the warming of the surface waters. We applied a two-epoch coalescent model (Shapiro et al. 2004), as implemented in BEAST, that simulated a two-parameter exponential growth (female effective population size scaled by mutation rate, fN e l1; and intrinsic growth rate, r) preceded by a one-parameter constant size model (fN e l0), with a parameter for the transition time between the two epochs (t_transition). We used the same molecular evolution model as for the BSP. 1/x priors were used for all population size parameters, and simple uniform priors were employed for r and t_transition. Upper and lower limits of the prior distribution for each parameter in the model were set using the 95% highest posterior density (HPD) intervals from the BSP as guideline. Run length was 20 million generations, with sampling every 1000th step and a burn-in of the first 25% of sampled generations (ESS > 200). The mutation rate was calculated as l = t_transition/c, where c is the calibration point. Assuming that population growth started either with rising sea levels after the LGM at around 18 ka BP or with surface water warming at 13 ka BP, we estimated l = t_transition/ 18 000 to t_transition/13 000.
To explicitly test whether the two-epoch model indeed fits the data better than a simple constant population size model, we employed a Bayes factors (BFs) approach for model selection (Suchard et al. 2001). Marginal likelihoods for the BF calculation were estimated under both the two-epoch and the constant population size model by path sampling (PS; Lartillot and Philippe 2006) and stepping stone sampling (SS; Xie et al. 2011), two approaches that were shown to provide accurate and reliable marginal-likelihood estimates for model comparisons (Baele et al. 2012). PS and SS were conducted in BEAST (20 million generations, 100 path steps, following a burn-in of 200 000 generations), adopting the XML codes for PS and SS from Baele et al. (2012Baele et al. ( , 2013.

Phylogeographic pattern
In total, 24 haplotypes were detected in 117 Mediterranean + Atlantic coast samples, which were clearly divergent from the three Atlantic island haplotypes (also see Domingues et al. 2007). The Croatian sample contained 16 different haplotypes, none of which were shared with T. delaisi from other parts of the Mediterranean or from the Atlantic (Fig. 2). Genetic landscape shape visualizations revealed no major phylogeographic breaks within Croatia with little genetic distance between sampling localities (Fig. 1c). The darker shading in the area between the islands Cres and Pag suggests that genetic structure is higher among the northern island locations than between these locations and the more southern populations at Ugljan and Hvar.
The majority of the Croatian haplotypes belong to a starshaped haplogroup, which is indicative of recent population expansion. Additionally, the Croatian sample contains several more divergent haplotypes at low frequencies. In part owing to these divergent haplotypes, the genetic diversity within the Croatian sample is nearly as high as in the entire Mediterranean + Atlantic coast sample (Croatia, H d = 0.544, p = 0.00448; Mediterranean + Atlantic coast, H d = 0.734, p = 0.00659). The apparently derived position of the divergent low-frequency haplotypes in the ML tree (Fig. 2) is most likely an artefact resulting from the low in-group divergence, the low statistical support for in-group nodes and from the comparatively large divergence from the out-group, all of which compromise the accurate positioning of the root (see e.g. Wheeler 1990; Kirchberger et al. 2014).

Substitution rate estimation by expansion dating
Tajima's D and Fu's F S statistics calculated with the sequences of the full Croatian sample were significantly negative (D = À1.6261, p = 0.0257; F S = À8.2318, p = 0.0012), which supports the scenario of recent population growth suggested by the star-shaped part of the haplotype network in Fig. 2. The BSP analysis based on all Croatian haplotypes clearly indicated population growth in the recent past (Fig. 3).
Model selection analysis based on BFs favoured the two-epoch model with constant population size followed by recent exponential population growth, which simulates the inferred BSP patterns, over a pure constant population size model (Table 1). The 2 ln BF estimates of 6.42 (SS) and 6.6 (PS) indicate strong support for the two-epoch model (Kass and Raftery 1995). Mean values for the transition time between constant population size and exponential growth were similar to the infliction point depicted in the BSP model (Fig. 3a). Posterior distributions for t_transition, fN e l0, fN e l1 and r were all unimodal. We used the transition time between the constant population size period and the recent phase of exponential growth for calibrating the substitution rate of the analysed control region fragment. The mean parameter estimate for t_transition was 4.6974 9 10 À4 . Assuming that the population expansion in the Adriatic Sea was associated with rising water levels after the LGM about 18 ka BP or with surface water warming at 13 ka BP, this translates into mean estimates of l of 2.61% (calibrated with water level rise) to 3.61% (calibrated with warming) per site per MY (95% HPD, 1.11-6.77%).

Phylogeographic pattern
Consistent with expectations based on the life history of T. delaisi and the geomorphological and hydrological features of the Adriatic Sea, mitochondrial sequences indicated recent population expansion, extensive haplotype sharing among locations connected by the mid-Adriatic cyclonic gyre and differentiation from conspecifics in the Mediterranean and the Atlantic. The majority of the Adriatic samples were joined in one single 'Adriatic' haplogroup featuring the star-like shape indicative of population expansion. Additionally, several individuals from various Adriatic locations carried different divergent haplotypes some of which were similar to those sampled in other parts of the Mediterranean (Cyprus, Italy). These divergent haplotypes might either represent polymorphisms dating from the original colonization of the Adriatic, or gene flow from another Adriatic gyre system or other parts of the Mediterranean. While a more thorough sampling from outside the current study area is required to distinguish between these alternatives, gene flow into our study area after its original colonization is certainly a realistic assumption. Moreover, if several divergent haplotypes had been present in the Adriatic prior to the inferred population expansion, it is difficult to explain why only one of them gave rise to a star-shaped haplogroup, while the others remained at low frequencies.
The phylogeographic pattern in our data relates to the characteristic circulation systems of the Adriatic basin (Astraldi et al. 1999). The Adriatic circulation systems are strongly influenced by wind conditions, Po River discharge and weather conditions in general, which cause seasonal patterns of variation in the flow of water currents (Ursella et al. 2007). Mean currents parallel to the shoreline (eastern Adriatic current, EAC; western Adriatic current, WAC) are weak or almost absent in spring and summer, while cyclonic gyres are present throughout the year in the middle and southern Adriatic, but are particularly well-developed in spring and summer (Russo and Artegiani 1996;Ursella et al. 2007; see Fig. 1). Reproduction in T. delaisi takes place in spring and summer (Wirtz 1978). Consequently, planktonic larvae are present in periods with strong cyclonic circulations in the Adriatic Sea. Along the Portuguese Atlantic coast, larvae of T. delaisi appeared to be constrained to inshore regions, and no larvae were found offshore (Borges et al. 2007). Consistently, high levels of self-recruitment in populations of the Spanish Mediterranean suggested that larvae remain near their natal areas   Carreras-Carbonell et al. 2007). With their inshore distribution, larvae of T. delaisi in the Adriatic may largely escape the offshore water currents, but will still be caught by local, typically variable and reversible currents running near the coast and between the islands (Shanks 2009). The lack of phylogeographic structure in the Adriatic sample certainly indicates that gene flow occurs throughout the investigated area, probably helped by both the middle Adriatic gyre and small local currents. Additionally, the genetic landscape shape analysis identified a slightly increased genetic structure among the northern populations. A possible explanation for this may lie in the sampling of some populations from small and sheltered bays (e.g. Uvala Pastura on Pag), which may experience less gene flow than populations that are more exposed to water currents.
While local currents promote dispersal of T. delaisi on smaller geographic scales, any planktonic larvae drifting offshore will probably be trapped in their cyclonic gyre and prevented from dispersal across gyres and beyond the Adriatic basin. Investigations of fine-scale population structure of organisms within the Adriatic basin are scarce (Maltagliati et al. 2010), but differentiation between populations in the Adriatic Sea and the rest of the Mediterranean, as observed in T. delaisi, has also been reported for several other organisms (e.g. planktonic chaetognaths, Sagitta setosa: Peijnenburg et al. 2004;common shrimp, Crangon crangon: Luttikhuizen et al. 2008; edible sea urchin, Paracentrotus lividus: Maltagliati et al. 2010;red mullet, Mullus barbatus: Maggio et al. 2009; Mediterranean rainbow wrasse, Coris julis: Fruciano et al. 2011;marbled goby, Pomatoschistus minutus: Mejri et al. 2011;small-spotted catshark, Scyliorhinus canicula: Gubili et al. 2014). Notably, all these taxa are poor active dispersers and/or rely on passive dis-

Expansion dating
Population expansion in the wake of rising sea levels following the lowstand during the LGM has been reported for numerous marine species (e.g. Marko et al. 2010;Naro-Maciel et al. 2011;Crandall et al. 2012;Ho et al. 2014), including some species from the Adriatic Sea (e.g. Debes et al. 2008;Luttikhuizen et al. 2008). Given that Adriatic T. delaisi also display genetic signatures of population expansion, we took advantage of the likely concomitance of population growth with the postglacial habitat expansion and applied an expansion dating approach to estimate the substitution rate of the mitochondrial control region in Trypterigiidae. Using the period between the LGM and the warming of surface waters as probable onsets of population expansion (18- Genner et al. 2010). In contrast, the range of substitution rate estimates previously used in studies of blennioids, including Tripterygiidae, is comparatively wide and ranges from 1.8 to 10% per site per MY (e.g. Domingues et al. 2007;Hickey et al. 2009;Francisco et al. 2011). These estimates were derived from the trans-Isthmian divergence of two closely related pomacentrid fishes (Domingues et al. 2005) in one case, and adopted from data of galaxiid fish in New Zealand (Waters and Burridge 1999;Burridge et al. 2008) in another. The internally calibrated estimates obtained in the present study are likely an improvement over values borrowed from distant species. Due to the time dependency of molecular rates (Burridge et al. 2008;Peterson and Masel 2009;Ho et al. 2011; but see Emerson and Hickerson 2015), the very recent calibration points used in our dating restrict the applicability of our substitution rate estimates to analyses of recent evolution, in the scales of tens and perhaps hundreds of thousands of years.

Concluding remarks
The present study is one of only a few to address the genetic structure of animals within the Adriatic basin. High rates of gene flow, likely mediated by planktonic dispersal, were inferred among locations within a cyclonic gyre. Further sampling across gyres will be required to test whether for passive dispersers, the cyclonic gyres partition the Adriatic basin into discrete phylogeographic regions, or whether boundaries between gyres can be overcome with the help of local inshore currents. The phylogeographic separation of the Adriatic Sea from the rest of the Mediterranean, which has been observed in a range of taxa with limited active dispersal abilities, is also evident in the data produced and used in the current study. Finally, we report a novel substitution rate for our study species, which will facilitate future phylogenetic and population genetic work within and beyond Tripterygiidae.