Combining population genomics with demographic analyses highlights habitat patchiness and larval dispersal as determinants of connectivity in coastal fish species.

Gene flow shapes spatial genetic structure and the potential for local adaptation. Among marine animals with nonmigratory adults, the presence or absence of a pelagic larval stage is thought to be a key determinant in shaping gene flow and the genetic structure of populations. In addition, the spatial distribution of suitable habitats is expected to influence the distribution of biological populations and their connectivity patterns. We used whole genome sequencing to study demographic history and reduced representation (double-digest restriction associated DNA) sequencing data to analyse spatial genetic structure in broadnosed pipefish (Syngnathus typhle). Its main habitat is eelgrass beds, which are patchily distributed along the study area in southern Norway. Demographic connectivity among populations was inferred from long-term (~30-year) population counts that uncovered a rapid decline in spatial correlations in abundance with distance as short as ~2 km. These findings were contrasted with data for two other fish species that have a pelagic larval stage (corkwing wrasse, Symphodus melops; black goby, Gobius niger). For these latter species, we found wider spatial scales of connectivity and weaker genetic isolation-by-distance patterns, except where both species experienced a strong barrier to gene flow, seemingly due to lack of suitable habitat. Our findings verify expectations that a fragmented habitat and absence of a pelagic larval stage promote genetic structure, while presence of a pelagic larvae stage increases demographic connectivity and gene flow, except perhaps over extensive habitat gaps.

INTRODUCTION

Gene flow influences evolutionary processes in multiple ways and may promote or constrain population structure and local adaptation (Lenormand, 2002; Slatkin, 1987). The life history and dispersal ability of species may determine the opportunities for gene flow in different environments. A key factor facilitating extensive gene flow is the presence of an early life‐history stage with high dispersal potential. In the marine environment, pelagic eggs and larvae represent such a potential with their ability to be transported by ocean currents, often in great numbers and over long distances. Many marine invertebrates and fishes have pelagic larvae and/or eggs, and this is widely regarded as an explanation for the generally observed low level of genetic structure among marine populations as compared to organisms inhabiting terrestrial or freshwater environments (Palumbi, 1994; Waples, 1998; Ward et al., 1994). Nevertheless, more recent studies, using new genetic and genomic tools, have uncovered genetic structuring in many marine species previously thought to be largely panmictic, including highly mobile organisms such as fishes (Hauser & Carvalho, 2008; Knutsen et al., 2003; Manel et al., 2020). While some of this genetic structuring has been related to environmental gradients (Catarino et al., 2015; Hellberg, 2009; Johannesson et al., 2020; Reid et al., 2016; Stanley et al., 2017; Wang & Bradburd, 2014), cryptic or hidden genetic heterogeneity is also observed in the apparent absence of environmental gradients (Knutsen et al., 2018; Selkoe et al., 2010), underlining our incomplete understanding of the mechanisms that restrict gene flow in coastal and oceanic systems.

A range of factors have been suggested to affect population structure in the marine environment. Apart from historical patterns of recolonization (Le Moan et al., 2019; Mattingsdal et al., 2020) and ongoing selective forces (isolation by environment [IBE; Berg et al., 2016; Wang & Bradburd, 2014), persistent genetic structure requires low levels of gene flow to shape the genetic constitution of populations. Despite transport and advection of propagules by ocean currents (Jahnke et al., 2018; Knutsen et al., 2009), levels of gene flow can be reduced due to life‐history differences among ecotypes (Kirubakaran et al., 2016), sex‐dependent dispersal ability (Ashe et al., 2015), phenotypically dependent traits (Fobert et al., 2019), bathymetric boundaries and ocean currents (Catarino et al., 2015; Riginos et al., 2016), and habitat patchiness (Binks et al., 2019). The distance between suitable habitats may influence gene flow at both the drifting larval stage and active adult dispersal (Blanco et al., 2016; Seljestad et al., 2020).

To identify the mechanisms responsible for cryptic or hidden genetic stratification, we need first to correctly identify the geographical scales of relevant ecological and physical processes (Allen & Hoekstra, 1992). The use of multiple methods may aid in identifying the appropriate geographical scale of dispersal and gene flow, especially when such dispersal occurs mainly at the pelagic egg and larval stages (Cayuela et al., 2018). Directly tracking dispersal of larvae is challenging due to their small sizes. One commonly used indirect approach is biophysical modelling, to estimate dispersal probabilities over single or multiple generations and how the metapopulation is shaped by circulation patterns and dispersal barriers (Galindo et al., 2006). Genetic methods are also commonly employed that estimate realized gene flow and thereby successful dispersal (Beerli & Felsenstein, 2001; Hey & Nielsen, 2004; Rousset, 1997; Slatkin, 1985, 1994). Despite various shortcomings, genetic methods remain important (e.g., Allendorf, 2017) and can be powerful tools when combined with other approaches that provide complementary information on dispersal and population connectivity (Hawkins et al., 2016; Lowe & Allendorf, 2010; Marandel et al., 2018; Palumbi, 1994). As indirect genetic methods are influenced by rare long‐distance dispersal (Waples & Gaggiotti, 2006), complementary methods are needed to fully understand how connectivity operates in the ocean. Such complementary methods include biophysical modelling of dispersal (e.g., Selkoe et al., 2008), seascape genomics estimating the correlation between environmental variables and the frequency of particular alleles or genotypes (Gebremedhin et al., 2009; Selmoni et al., 2020), as well as estimates of spatial population synchrony (Östman et al., 2016; Rogers et al., 2014) and habitat availability (Buonomo et al., 2017; Selkoe et al., 2010).

Dispersal among populations not only has genetic consequences but may also change local demography (Hastings, 1993) and lead to spatial synchrony in population dynamics (Kendall et al., 2000). Spatial population synchrony, defined as the spatial covariation in fluctuations of population abundance (Bjørnstad et al., 1999), can be expected among population segments that exchange individuals through dispersal and experience similar environmental conditions affecting survival. Fluctuations in abundance typically become less synchronous as the distance between population segments (or sampling locations) increases, reflecting decreased connectivity through dispersal and (possibly) less similar environment conditions. This phenomenon opens another avenue to investigate dispersal in the marine environment, by analysing spatial abundance data and assessing the spatial scale of synchrony among sampling sites (Hameed et al., 2016; Östman et al., 2016; Rogers et al., 2017). Population connectivity inferred from demographic correlations can provide an independent indicator of dispersal that complements and may be contrasted with genetic patterns (Cayuela et al., 2018).

Here, we estimate and compare spatial genetic structure and demographic decorrelation scales (geographical distance at which correlation in annual fluctuations in abundance falls below a certain threshold value) in three coastal fish species. We provide new genetic and demographic data for broadnosed pipefish (Syngnathus typhle) and contrast it with two unrelated fish species, the corkwing wrasse (Symphodus melops) (Blanco et al., 2016) and the black goby (Gobius niger) (Catarino et al., 2022). All three species were sampled over ~900 km of coastline, typically from the exact same locations or nearly so, providing a unique opportunity for assessing drivers of population stratification in a consistent sampling design. The black goby and the corkwing wrasse both have larvae that remain pelagic for 14–28 days and are putatively subject to advection by ocean currents. In contrast, the broadnosed pipefish carry the eggs until they hatch and nurse the larvae until released as free‐swimming small juveniles.

Using a comparative approach, we test the hypothesis that species lacking a pelagic larval phase (pipefish) have reduced dispersal capability as compared to species (corkwing wrasse and black goby) that have pelagic larvae and thus a greater potential for dispersal with ocean currents. The hypothesis is tested by comparing the demographic decorrelation scale (Bjørnstad et al., 1999; Rogers et al., 2017) and the slope of the genetic pattern of isolation‐by‐distance (IBD) (Rousset, 1997). We search for genomic footprints of selection and reconstruct demographic history by analysing whole genome sequences, the latter with the aim of revealing possible deeper historical factors underlying genetic structure. We can thus detect if the present genetic structure is a result of historical contingencies (e.g., different postglacial colonization patterns) or, alternatively, caused by contemporary processes (gene flow, genetic drift and selection). The findings are discussed in a broader context of how dispersal capability, gene flow and local adaptation shape the genome‐wide genetic structure of a coastal fish species.

MATERIAL AND METHODS

The study was carried out along the coast of southern Norway, including the southern (Skagerrak) and western coasts (Figure 1). The target species for the present study, the broadnosed pipefish (Syngnathus typhle), is distributed in the Northeast Atlantic along coastal areas from Vardø in northern Norway to Morocco in the south (Dawson, 1986). This species prefers seagrasses or other macrophyte habitats down to about 20 m depth. Broadnosed pipefish are normally <30 cm in size in Norwegian waters. They are poor swimmers and are not expected to swim long distances. Rafting with algae may be possible but is still undocumented in northern waters (Bertola et al., 2020). The species spawns from late June to early July. Eggs from one to three females are deposited in the male's brood pouch, in batches of around 20 where they are fertilized. The fry hatch after about 4 weeks, then may linger around the male and retreat to the safety of the pouch should danger threaten. They are typically 1–3 cm in size and have some ability to swim, but are fragile (Foster & Vincent, 2004).

FIGURE 1

Map over all sampling localities (red dots) with the names of the location. Orange arrows denote Atlantic water masses coming towards Norway, while green arrows denote the main pathway for the Norwegian coastal current (redrawn from Sætre, 2007, p. 100). See Table 1 for details of the number of fish from each location

TABLE 1

Sampling sites of broadnosed pipefish for genetic analyses, arranged geographically (samples from Søgne and down are within Skagerrak)

Locality	No. of stations	ddrad sample size, n		Whole genome sample size, n	N total	F IS	H E
		2018	s
Ålesund	1	10	10	4	20	−0.024	0.0735
Bergen	2	4	21	4	25	−0.103	0.0796
Stavanger	1	23	0	4	23	−0.006	0.0720
Egersund	1	26	0	4	26	0.007	0.0626
Søgne	4	5	9	4	14	0.034	0.0671
Kristiansand	2	2	0		2	−0.052	0.0817
Lillesand	6	21	11	4	32	−0.014	0.0702
Tvedestrand	2	7	3	4	10	0.037	0.0682
Risør	8	17	13		30	−0.010	0.0699
Kragerø	6	28	14		42	−0.018	0.0717
Eidanger	5	14	8		22	−0.020	0.0706
Sandefjord	4	25	15		10	−0.020	0.0728
Færder	3	21	0		21	0.023	0.0668
Vrengen	4	17	0		17	0.031	0.0660
Oslo	5	24	0		24	−0.003	0.0681
Hvaler	3	28	0	4	28	−0.017	0.0702

No. of stations is the number of (closely located) beach seine sampling sites. F IS is the fixation index. H E is genetic diversity (expected heterozygosity under Hardy–Weinberg equilibrium).

Two other coastal fish species have been analysed in recent studies from the same area, allowing comparison with the broadnosed pipefish with respect to genetic structure and demographic scaling. For comparative purposes, we reiterate some findings of those studies herein and supply some new analyses (spatial decorrelation scale: below), while referring to the original publications for details. The corkwing wrasse (Symphodus melops: Knutsen et al., 2013; Blanco et al., 2016; Faust et al., 2021) and black goby (Gobius niger; Catarino et al., 2022) are common in coastal waters of southern Norway. Like the broadnosed pipefish, both prefer eelgrass but may also be found in habitats with macroalgae in sandy or muddy sheltered areas. Both species build nests and guard their eggs until hatching into larvae that are pelagic for 15–25 days for corkwing wrasse and slightly longer (mean about 28 days) for black goby (Darwall et al., 1992; Planes, 1998).

Sampling

We took advantage of a uniquely long and continuous time series of standardized juvenile fish collections in coastal nursery habitats in coastal Skagerrak (Stenseth et al., 1999). The abundance of broadnosed pipefish has been quantified in numbers (1989 to present), by recording the number of individuals in beach seine hauls conducted in September–October of each year. Biogenic habitat type (e.g., eelgrass or type of macroalgae) and coverage (scaled from 1 to 5 subjectively: 1 = bare substrate to 5 = full coverage) have been characterized at each haul location. The beach seine is 40 m long, has a stretched mesh size of 1.5 cm, and covers an area up to 700 m2 of nearshore (<15 m depth) habitat. The beach seine survey normally encompasses four to 10 stations in each fjord or segment of the coast. Samples from each sample location in the present study (cf. Table 1) normally included fish from more than one beach seine station to increase sample sizes for genetic analyses. Spatial decorrelation analyses (below), however, were done on individual stations to maximize spatial resolution. To extend the study area outside the long‐term beach seine monitoring programme in the Skagerrak, we performed additional beach seine sampling in 2018 and 2019 at four locations along the Norwegian west coast (Figure 1: localities Egersund, Stavanger, Bergen and Ålesund). Together, we sampled a total of 367 individual broadnosed pipefish from 16 localities for genetic analysis (Table 1). Individual fish were immediately killed upon catch to obtain sufficient tissue of both flesh and skin ensuring that we obtained high enough DNA quality for the genetic analysis. Fish were collected and transferred in the field to tubes with 96% ethanol. The tubes were stored at 4°C until DNA extraction.

Habitat associations

We used a binomial regression with a random intercept by station to test for effects of biogenic habitat type and coverage on presence or absence of broadnosed pipefish. Although recorded as a categorical variable with five categories (“bare substrate” to “full cover”), coverage was treated in the model as a continuous variable (ranging from 1 to 5) after preliminary model runs indicated an approximately linear relationship between presence and coverage. For this analysis, we included data from all stations surveyed for at least 15 years, including stations with few to no positive catches, in order to retain information on conditions leading to absence of pipefish. Models were fitted using the lme4 package (Bates et al., 2015) in the R statistical environment (R Core Team, 2020).

Whole genome sequencing

To reconstruct the demographic history among samples and to scan for potential candidate loci or groups of loci under selection, we carried out WGS genome sequencing (WGS) on a subset of samples. For this, a subset of four pipefish individuals were taken from eight of the 16 locations, resulting in 32 individuals for WGS (cf. Table 1). DNA from the pipefish samples were extracted using an E‐Z 96 Tissue DNA kit or E.Z.N.A, E‐Z 96 Plant DNA kit (Omega Bio‐tek). All samples were treated with RNAse A (#19101). DNA concentration (ng µl–1) was measured either on a Qubit or a Fluoroskan instrument (Thermo Fisher Scientific). Subsamples were tested for DNA degradation and purity by using a NanoDrop instrument (Thermo Fisher Scientific). Library construction was performed using the Illumina TruSeq DNA PCR Free protocol and checked on a Bioanalyser High sensitivity chip and Tapestation (both Agilent) followed by Kapa Biosystems’ qPCR (quantitative polymerase chain reaction) assay for quantification of Illumina libraries. Whole genome resequencing was conducted on the Illumina HiSeq platform at the Norwegian High‐Throughput Sequencing Centre (NSC: Oslo, Norway). The paired‐end sequencing reads from the 32 individuals were aligned to the broadnosed pipefish reference genome (Roth et al., 2020, GenBank accession no. GCA_901007915.1) (N50 = 485 kb) using a Burrows‐Wheeler Aligner BWA (bwa mem: Li, 2013). Amongst all 848,053,562 reads (mean 26,501,673 per individual), 95.2% mapped to the genome, of which 702,372,912 (85.5%) were properly paired. The mean depth of coverage across the genome was 10.9×. Sequence variants were detected by freebayes (Garrison & Marth, 2012), which reported 3,896,949 variations. Of these, 2,837,389 were bi‐allelic single polymorphism nucleotides (SNPs) with a quality score >40. We applied minimum and maximum read depth filters by setting calls with depth of ≥6 and maximum depth of ≤30. SNPs with >10% missing data were then removed, resulting in 431,910 SNPs, at an SNP density of one per727 bp.

Demographic history reconstruction

We extracted SNPs from whole genome data on contigs larger than N50 > 3,046,963 bp. These were 30 contigs representing half the genome of the species (158 Mb). We used the software smc++ (Terhorst et al., 2017), and calculated the composite likelihood across four individuals (as recommended by the manual) using a mutation rate of 1 × 10−8 per base per generation. Changes in effective population size through time were estimated separately for western and southern sites. To estimate variance and to ascertain convergence of the inferred historical profiles, we performed 10 iterations. The pipefish is sexually mature after 1 year and lives for 2–3 years, so we set the generation time (i.e., mean age of parents when their offspring are born) to 3 years (Dawson, 1986) in the analysis.

Scan for selection

We conducted a scan for candidate regions under evolutionary constraint for the full‐genome data (see Supporting Information for similar analysis and results from the double‐digest restriction associated DNA [ddRAD] data). Rare variants were excluded (minimum allele count, MAC > 2), and we allowed a maximum 10% missing data, resulting in 97,749 SNPs from the full‐genome sequences. We interrogated this data set for signals of selection using bayescan (version 2.1., Foll & Gaggiotti, 2008). bayescan estimates F ST values which are decomposed to population‐ and locus‐specific components, where we compared the Skagerrak samples and the Western samples. The script vcf2bayescan (Sánchez‐Ramírez et al., 2018) was used to prepare input files. bayescan was run with default settings (‐n 5000 ‐thin 10 ‐nbp 20 ‐pilot 5000 ‐burn 50000), except increasing the Bayes factor from 10 to 100 to adjust for the number of SNPs (>1000). Putative outliers were detected by plotting the F ST coefficient against the posterior odds and q‐values against alpha values. We used the commonly used thresholds for significance to q < 0.05 and log10 p > 2. To adjust for multiple testing, we set the false discovery rate (FDR) level to 5% using the R package q‐value (Storey et al., 2021).

Reduced genome sequencing (ddRAD) and SNP calling

DNA from all samples was normalized to the same concentration of 7 ng µl–1. ddRAD (Peterson et al., 2012) sequencing was used to detect SNPs. The restriction enzymes applied were Sbf1 and Sph1, and fragment size ranged from 250 to 600 bp. Each library contained 96 double‐indexed specimens and was sequenced by 150‐bp paired‐end sequencing on an Illumina Miseq V3 flow cell (Illumina) spiked with 5% PhiX.

SNP calling

Reads for individual specimens were retrieved with the process_radtags module in the stacks software package version 1.48 (Catchen et al., 2013). Retrieved reads were aligned to the reference genome sequence (link given above) with bwa (Li & Durbin, 2009). Initial SNP filtering was conducted with the populations module to only include SNPs present in at least 10% of the individuals from at least one location.

SNP filtering

We applied vcftools (versionb 0.1.16: Danecek et al., 2011) with the following set of parameters to filter the data further (starting with 14,850 loci, 302 individuals). First, we set the maximum missing data as 30% for loci (down to 10,887 loci, 302 individuals) and 40% for individuals (10,887 loci, 283 individuals). Then, we applied the MAC =2 filter to include only bi‐allelic loci, HWQ to 0.001 to remove strong deviations from Hardy–Weinberg equilibrium (HWE), and finally, we applied plink version 1.07 (Purcell et al., 2007) to remove loci exhibiting linkage disequilibrium (LD >0.5 level) to filter and retain only nonlinked SNPs in HWE, resulting in 1,461 SNPs and 283 individuals for downstream analyses.

Genetic population structure

Allelic frequencies, expected heterozygosity (H E) and fixation index (F IS) were estimated for the ddRAD SNPs using genepop software (version 4.7.5: Raymond & Rousset, 1995, Rousset, 2008). Genetic differentiation among all samples (overall F ST) and pairwise F ST between pairs of samples were estimated using Weir and Cockerham's (1984) estimator theta (Ꝋ). Tests for genetic heterogeneity among all samples and for pairs of samples were carried out with the allele frequency test option in genepop, using default numbers of iterations and batches. To identify spatial patterns in genetic structure we used the R package adegenet 2.1.1 (Jombart & Ahmed, 2011) to perform principal component analysis (PCA) on all 1,461 SNPs. Further spatial genetic analyses (below) focused on genetic structure as a function of geographical distance and habitat availability.

Drivers of population structure

Due to the lack of a pelagic larval stage, we hypothesized that dispersal and gene flow in broadnosed pipefish are limited by nearby habitat availability. To test potential drivers of genetic structure, we adopted the maximum likelihood population‐effects (MPLE) approach (Clarke et al., 2002), with sample sites (“populations”) as random effects and environmental factors (geographical distance and habitat gaps) as fixed effects. This approach accounts for nonindependence in pairwise data by including the correlation structure among sample localities. We used the cormple function in the R package with the same name (Pope, 2022). Geographical distances over water among all samples were calculated using the R package marmap (Pante & Simon‐Bouhet, 2020) with bathymetric data from the ETOPO1 database (Amante & Eakins, 2009) hosted by NOAA (https://www.ngdc.noaa.gov/mgg/global/relief/ETOPO1/docs/ETOPO1.pdf). Slight adjustments of sample positions were necessary to avoid samples on land due to limited resolution of the bathymetric data in coastal waters. Because our habitat analysis (above) showed a strong preference in this species for eelgrass, we took eelgrass beds as suitable habitats. The positions and sizes of eelgrass habitats were taken from a previous habitat mapping project, available at: https://kart.kystverket.no/share/9220e0e277e4. We used these positions to calculate the longest span without an eelgrass bed (“habitat gap”) between pairs of sample localities, as measured over the shortest water path connecting them. We included geographical distance (over water) and habitat gaps in a generalized least squares model with pairwise F ST as the response variable in the gls function in the R package nlme (Pinheiro et al., 2022): where pop1 and pop2 denote the two samples in each pair. Submodels without the interaction term between distance and habitat gaps, and with only one of the two factors alone, were constructed and the optimal model was chosen on the basis of the corrected Akaike Information Criterion (AICc), using the model.sel function in the R package mumin (Bartoń, 2020).

Spatial population synchrony and decorrelation scale

Dispersal among population segments is expected to reduce genetic differentiation and increase the spatial scale of population synchrony. This suggests that species with extended pelagic larval durations should remain correlated in their abundance fluctuations at longer distances than those with no or limited dispersal, reflecting differences in demographic connectivity. We compared the spatial population synchrony and decorrelation scales among the three study species to determine whether the data were consistent with this expectation.

To calculate spatial population synchrony, we calculated the pairwise correlation in natural log‐transformed catches (counts) among each pair of sampling stations. A constant of 1 was added to each count to avoid taking the logarithm of zero. Some stations were not sampled every year or had catches of zero in most years. These were excluded from the synchrony analysis by restricting stations to those sampled a for a minimum of 15 years and having at least 10 years with positive (nonzero) catches. This resulted in a total of 61, 123 and 132 stations for broadnosed pipefish, corkwing wrasse and black goby, respectively, distributed along the Skagerrak coast from Kristiansand in the west to Hvaler in the east (cf. Figure 1).

For each species, an exponential decay model (Bjørnstad et al., 1999) was fitted to pairwise correlations by distance, measured as distance over water between stations: where ρ(d) is the pairwise correlation at distance d, ρ r is the asymptotic correlation, or background regional correlation in recruitment, and ρ r + ρ 0 is the estimated correlation at zero distance, following Rogers et al. (2014). The parameter v (“decorrelation scale”) describes the distance at which the pairwise correlation between time series is reduced to ~37% (e−1) of that at zero distance, relative to the background regional level of correlation. Confidence intervals (90%) were estimated based on 10,000 bootstrap replicates, following Bjørnstad et al. (1999).

Species comparisons

Spatial population structure in the broadnosed pipefish was contrasted with published genetic data from the corkwing wrasse and the black goby (Blanco et al., 2016; Catarino et al., 2022) and archived abundance data from the long‐term beach seine survey, described above. For this, genetic structure along the coast was analysed by regressing pairwise F ST between the easternmost sample (at Hvaler or the closely situated Oslofjord) with the other samples. This approach retained geographical context to the genetic patterns and revealed a marked genetic discontinuity or “break” at the transition from Skagerrak to the west coast (between Egersund and Stavanger: Figure 1). We therefore applied a nonlinear (sigmoid) regression to analyse the relationship between genetic differentiation (F ST) and geographical distance (nls function in R). This was supplemented with linear regressions done separately for the two coastal sections on either side of the break (lm function in R). Statistical tests for significance employed the t‐statistics reported by the lm function. No tests were carried out for the west coast section because of few sample localities for two of the species.

The spatial scale of demographic influence was compared among species using the estimated spatial demographic decorrelation scales (above). These estimates were restricted to the Skagerrak (cf. Figure 1; Table 1) and used broadly overlapping sample sites with the genetic analyses.

RESULTS

Analysis of beach seine catches showed a strong and significant effect of biogenic habitat type on the presence or absence of broadnosed pipefish (likelihood ratio test; p < 0.001). The pipefish showed a strong affinity for eelgrass and was more often present and abundant in stations with eelgrass than in those without (Figure S1). In habitats with green algae (with or without brown macroalgae) and no eelgrass, pipefish were rare, although these habitats were infrequently observed in the beach seine survey (n = 33 out of 2741 total hauls). We also found a significant positive effect of coverage on abundance of broadnosed pipefish (z = 7.06, p < 0.001). Pipefish were rarely caught in habitats with bare substrate, and increasingly present as biogenic habitat cover increased (Figure S1).

Demographic history reconstruction based on whole genome sequences

The estimated demographic history, that is profile of effective population size over time (N: Figure 2), using whole genome data was made separately for the Western Norway and for the Skagerrak samples. Both groups were relatively consistent for the 10 iterations and suggested a shared history of post‐glacial expansion for both groups, ~10,000 years before present (years bp). Although estimating recent history is more uncertain by this method, the analysis suggested a decline in effective population size in the western sample with an onset about 1000 years bp, assuming a generation‐time of 3 years (cf. Figure 2).

FIGURE 2

The estimated demographic history using smc++ on whole genome sequencing data from individuals from two populations of pipefish in Norway using 10 iterations each. The samples included individuals (n = 4) from Ålesund (the western Norwegian coast) and Tvedestrand (the Skagerrak coast: Figure 1), respectively. The results suggest a shared history, with a decline in the west at 1000 years before present, possibly causing increasing genetic drift over this period

Scan for selection in whole genome data

We obtained few (112 out of ~100,000 loci, <1% after FDR), but some highly significant outliers from bayescan (Figure S2). Further, one genomic region was found to be monomorphic in the western samples (scaffold00010: 185,320–188,302). Although this region is modest in size (~3000 bp), it contains 23 SNPs. We extracted the sequence from the region and by using ncbi blast found high homology to the greater pipefish, Syngnathus acus, chr22: 10,218,727–10,220,596. This region contains the bend3 gene, involved in protein regulating transcription and in chromatin remodelling. This region obtains extreme p‐values in bayescan and is not included in Figure S2.

There was no sign of large blocks of elevated F ST indicative of chromosomal rearrangements in our data, and the Manhattan plot suggests that genomic divergence is evenly spread out over the genome (cf. Figure S3).

ddRAD analysis

For the ddRAD data, our sample of 367 individuals was reduced to 302 after filtering out those with many missing values (see Material and Methods). For these, we recovered 1,461 SNPs that met the specified requirements for missing values, minor allele frequency and deviations from HWE (above). Genetic variability (H E) in the different samples is given in Table 1 and ranged from 0.0626 (Egersund) to 0.0817 (Kristiansand, with only two individuals scored). F IS values representing an excess or deficit of heterozygotes were low, with 10 localities displaying a negative average F IS and five localities a positive F IS. SNPs had already been filtered for compliance to HWE, and consequently no further tests for HWE were performed.

Overall genetic differentiation among broadnosed pipefish samples was moderate, but highly significant (average among localities: F ST =0.0342, p ≪ 0.001). The genetic structure as depicted by the PCA (Figure 3) displayed two or three major clusters, represented by a dense southern (Skagerrak) cluster (at right in the figure), a more sparse western cluster (left) and with the Egersund sample somewhat of an outsider (top). This overall pattern was reflected in the pairwise F ST estimates (Table 2): among Skagerrak samples F ST ranged from −0.0005 to 0.0384 (mean 0.0133; p ≪ 0.001), whereas among western localities pairwise F ST was higher and ranged from 0.0267 to 0.1114 (mean 0.0646; p ≪ 0.001).

FIGURE 3

PCA plot of broadnosed pipefish from the study area, grouped by sample locality (Table 1, Figure 1), using the 1461 SNP ddRAD data set

TABLE 2

Pairwise F ST among all sample localities averaged over all 1461 SNPs. Bold numbers denote significant values (q < 0.05)

	Ålesund	Bergen	Stavanger	Egersund	Søgne	Kristiansand	Lillesand	Tvedestrand	Risør	Kragerø	Eidanger	Sandefjord	Færder	Vrengen	Oslo
Ålesund
Bergen	0.0390
Stavanger	0.0267	0.0305
Egersund	0.0748	0.0924	0.0500
Søgne	0.0663	0.0711	0.0514	0.0873
Kristiansand	0.0602	0.0837	0.0491	0.1114	0.0034
Lillesand	0.0725	0.0768	0.0447	0.0857	0.0166	0.0196
Tvedestrand	0.0541	0.0595	0.0314	0.0705	0.0118	0.0113	0.0149
Risør	0.0655	0.0705	0.0425	0.0790	0.0207	0.0345	0.0265	0.0067
Kragerø	0.0664	0.0752	0.0457	0.0817	0.0145	0.0188	0.0110	0.0095	0.0201
Eidanger	0.0734	0.0749	0.0518	0.0911	0.0236	0.0257	0.0218	0.0131	0.0343	0.0121
Sandefjord	0.0558	0.0648	0.0363	0.0786	0.0072	0.0125	0.0099	0.0013	0.0166	0.0048	0.0084
Færder	0.0631	0.0717	0.0383	0.0740	0.0116	0.0223	0.0106	0.0100	0.0194	0.0070	0.0135	0.0014
Vrengen	0.0599	0.0752	0.0422	0.0833	0.0104	0.0124	0.0084	0.0046	0.0203	0.0050	0.0170	0.0025	−0.0005
Oslo	0.0659	0.0740	0.0417	0.0807	0.0191	0.0384	0.0173	0.0149	0.0232	0.0118	0.0186	0.0058	0.0028	0.0025
Hvaler	0.0654	0.0728	0.0393	0.0752	0.0133	0.0303	0.0127	0.0047	0.0207	0.0074	0.0164	0.0015	−0.0007	0.0025	0.0078

The MLPE modelling of pairwise F ST against environmental drivers yielded highly significant effects of geographical distance and of habitat gaps between samples, both considered separately and in combination (Table 3). The best fit model included both factors without the interaction term, and demonstrated that geographical distance and habitat gaps both contributed to increased F ST between samples (Figure 4). This model described observed pairwise F ST reasonably well, except for sample pairs that included the Egersund sample (indicated with square symbols in Figure 4) and which displayed higher F ST values than predicted.

TABLE 3

Results of the maximum likelihood population‐effects (MLPE) models which were fitted to explain genetic divergence (F ST) between sample pairs related to two environmental drivers (geographical distance over water, and gaps in available eelgrass habitat)

Model	Estimated effects				df	AICc	ΔAIC
	Intercept	Distance (×1000)	Habitat gap (×10)	Interaction term (×1000)
Distance	0.0122	0.0860*	na	na	4	−740.4	29.43
Habitat gaps	0.0093	na	0.0120*	na	4	−743.9	25.97
Distance and gaps	0.0073	0.0502*	0.0067*	na	5	−769.8	0.00
Distance and gaps w/int.	0.0090	0.0407*	0.0058*	0.0003	6	−740.4	29.44

The best fit model according to AICc is highlighted in bold. Significant estimates are marked with an asterisk. na = not applicable.

FIGURE 4

Observed genetic divergence (F ST: dots) as a function of geographical distance and size of largest habitat gap between pairs of samples (both in km). The optimal MLPE regression model (cf. Table 3) is depicted as a plane, with observed F ST falling below the model prediction coloured grey and those lying above in black. Square symbols depict sample pairs with the Egersund sample

Contrasting genetic divergence patterns in the broadnosed pipefish with black goby and corkwing wrasse inhabiting the same coast revealed differences among species with respect to the shape and steepness of the regression of F ST with distance (Figure 5). For this analysis we compared F ST for each sample location with the easternmost sample for each species, which for all three was represented by the Hvaler sample or close by in the outer Oslofjord. Both the corkwing wrasse and the black goby displayed a marked genetic discontinuity at approximately the same geographical position, at the southwest tip of the coast, whereas the broadnosed pipefish changed somewhat more gradually (Figure 5a). Because of the geographical discontinuities within species, we compared species after splitting the coast into two sections for each species, corresponding to the Skagerrak and the west coast, respectively. Within the Skagerrak, genetic divergence among broadnosed pipefish samples increased with distance by 3.7 × 10−5 km–1 (p = 0.005) whereas this trend was an order of magnitude smaller for the two other species (black goby: slope −0.84 × 10−5 km–1, p =.29; corkwing wrasse: slope 0.19 × 10−5 km–1, p =.78: Figure 5b). Along the west coast all three species displayed an apparent increase in F ST with distance but the number of samples (three) was insufficient for a meaningful comparison.

FIGURE 5

Comparing patterns of genetic divergence (F ST) among three fish species inhabiting the same coast. For each species, pairwise F ST between the easternmost (Hvaler or Oslofjord, cf. Figure 1) and the other samples (symbols) is plotted against geographical distance over water. (a) Nonlinear regression (using the nls function with a sigmoid response curve in R) illustrating the abrupt shift or “break” in genetic divergence at around 400 km (between Egersund and Stavanger) for the corkwing wrasse and the black goby. (b) The same data as in (a), with separate linear regressions on either side of the break. Note the marked increase in F ST with distance for pipefish, while there were nonsignificant trends in IBD for black goby and corkwing wrasse on either side of the break

All three species showed a decrease in demographic correlation (i.e., density synchrony) with increased geographical distance (Figure 6). Regional background correlations in abundance were above zero for all three species, and highest for corkwing wrasse (ρ r = 0.20; 90% confidence interval [CI] = 0.17–0.22), and slightly lower for pipefish (ρ r = 0.08; 0.04–0.11) and black goby (ρ r = 0.10; 0.05–0.13). While fluctuations in abundance were more highly correlated at nearby stations for all species, the estimated decorrelation scale (v) was smallest for broadnosed pipefish (2.4 km; 90% CI = 1.3–21.5 km), followed by corkwing wrasse (13.5 km; 90% CI = 6.3–24.9 km) and black goby (27.9 km; 90% CI = 13.8–94.8 km) (Figure 7). Confidence intervals for the bootstrapped difference in estimated decorrelation scales between the broadnosed pipefish and black goby did not overlap zero, suggesting a significant difference between decorrelation scales for these two species. Using this same approach, the decorrelation scale for corkwing wrasse was not significantly different from the other two species. Results were robust to our choice of which stations to include based on a minimum number of years with positive (nonzero) catches. The shorter decorrelation scale in pipefish implies finer‐scale structuring in population dynamics of pipefish relative to the black goby, with corkwing wrasse falling between the two species.

FIGURE 6

Pairwise correlations among time‐series of catches at stations along the Norwegian Skagerrak coast (cf. Figure 1) (grey points here) for pipefish, black goby and corkwing wrasse. Solid coloured lines show a fitted exponential decay model along with bootstrapped 90% confidence intervals (shaded areas). Confidence intervals are based on 10,000 bootstrap replicates (Bjørnstad et al., 1999)

FIGURE 7

Estimated decorrelation scale (v) in kilometres (km) for three species based on the distance–decay of correlations in beach‐seine catches in the Skagerrak. Points show model estimates, thick bars show bootstrapped 25% and 75% quartiles and thin lines show bootstrapped 5% and 95% quantiles

DISCUSSION

Despite the tremendous increase over the last decade in amounts of genetic information generated, we still lack a firm understanding of the importance of the various mechanisms that restrict gene flow and demographic connectivity among marine fish populations. The present study improves our understanding of such mechanisms by contrasting the geographical scale of genetic structure and demographic decorrelation in three fish species within a common sampling design along the same coastal environment. The existence or absence of a pelagic larval stage (none of the three species have pelagic eggs), combined with geographical distance and habitat patchiness (“gaps”), influences the genetic structure and demographic independence of neighbouring populations in the broadnosed pipefish.

The main study area, and the coastal section from which abundance data are available, was the Skagerrak coast. Here, we found significant genetic heterogeneity in the broadnosed pipefish (F ST =0.0342, p ≪ 0.001), with a pattern of increased genetic divergence with distance (IBD slope =0.00010 km–1). A similar pattern has been documented for other pipefish species elsewhere (Graaf, 2006). Interestingly, eelgrass, the main habitat of broadnosed pipefish, shows a similar genetic structure in the Skagerrak and the Norwegian coast with significant IBD and fjord‐scale differentiation (Jahnke et al., 2018, 2020; Olsen et al., 2013). An alternative explanation to increased genetic divergence with distance could be that local adaptation and/or historical contingency underlie the genetic structure. There was, however, inconclusive evidence of selection, probably due to small sample sizes, suggesting that local adaptation is not a major factor in explaining the observed genetic structure (Figures S2 and S4). Further, demographic history analyses indicated that pipefish populations along the coastline have a shared history over a long timescale (Figure 2). A similar observation was made for corkwing wrasse over the same coastline (Mattingsdal et al., 2020), suggesting that historical events cannot explain the patterns of genetic structure in these two species (similar demographic history analyses are not presently available for the black goby). Instead, the observed genetic differences appear to be modulated by contemporary processes.

The availability of suitable habitat patches has been suggested to be important for successful dispersal and gene flow in other marine species (D’Aloia et al., 2015; Selkoe et al., 2010; Van Wynsberge et al., 2017). We found that both geographical distance and the length of gaps in suitable habitat (here: eelgrass beds) had significant, positive effects on genetic difference between sample localities (Table 3; Figure 4). This finding is in accordance with our hypothesis, based on their lack of pelagic eggs and larvae and on the limited swimming ability of adults: gene flow is limited in the broadnosed pipefish, and both geographical distance and lack of suitable habitat are limiting factors for gene flow.

The genetically outlying Egersund sample indicated particularly limited connectivity with neighbouring populations (Figure 4). This locality was sampled inside a confined bay, and there may be limited exchange with the surroundings. Eelgrass patches are not very large and the number of inhabitants per patch is likely to be restricted, giving ample opportunity for either genetic drift or bottlenecks promoting divergence. Indeed, this sample had the lowest heterozygosity (H E) of our samples (cf. Table 1), lending support to the notion of a semi‐isolated, small population. It is quite possible that such small populations go extinct from time to time and that the habitat patch subsequently becomes recolonized from a neighbouring population(s), thus constituting a metapopulation system.

The corkwing wrasse and the black goby displayed weaker genetic structure within Skagerrak as compared to the broadnosed pipefish, although still statistically significant (Blanco et al., 2016; Catarino et al., 2022; Knutsen et al., 2013) and either no or very weak IBD (cf. Figure 5), and indicating restriction in gene flow also in these two species. The notable genetic”break” between Skagerrak and western samples (cf. Figure 5) provides additional evidence for restrictions in gene flow in all three species. The mechanism responsible for this and other genetic break zones in the corkwing wrasse appears to be associated with areas that lack suitable hard bottom substrate that the species requires (Blanco et al., 2016; Faust et al., 2021; Mattingsdal et al., 2020), implying that lack of suitable habitats restricts gene flow and thereby shapes the genetic structure of this species. Limited dispersal abilities for the two species are corroborated with the estimated demographic decorrelation scales which, while larger than for the broadnosed pipefish, are still restricted to a few tens of kilometres (Figures 6 and 7), which may be insufficient to cross the genetic break zone. It may seem unexpected that the corresponding genetic break in the broadnosed pipefish is comparatively shallower and apparently shifted somewhat eastward in geographical position. A plausible explanation is that the overall greater genetic divergence among samples elsewhere in the broadnosed pipefish masks a genetic pattern that is more accentuated in the other two species because of their lower levels of divergence elsewhere. Also, the position of the outlying Egersund sample complicates the characterization of the genetic break in this species. All three species reach notable levels of genetic divergence when the whole 900‐km coastline is considered (cf. Figure 5), and the presence or absence of a distinct genetic break seems to have little or limited impact on the maximum levels of divergence reached.

We combined both population genomics and population demographic approaches to determine how gene flow and dispersal shape population structure in three fish species. While dispersal affects spatial genomic structure in the long term, it will have a spatial demographic fingerprint in the short term, as reflected in the scale of population synchrony and decorrelation. Alone, spatial decorrelation scales may be difficult to interpret for underlying mechanisms as there are typically multiple mechanisms at play: populations may be synchronized by dispersal, as well as by shared environmental or food‐web effects. However, by comparing decorrelation scales of multiple species in the same system, sampled at the same locations and subject to the same scales of environmental variability, we can infer that differences in decorrelation scale are related to the species life histories (e.g., pelagic larval duration) and not environmental influences. This approach relies on spatial abundance data collected over a number of years. While we had 30 years of data available for this study, a shorter duration would still be informative, especially for commonly occurring species. Broadnosed pipefish abundances were low in our study, which decreased certainty in decorrelation scale estimates due to relatively high observation noise associated with our fine‐scale sampling units. This will be less of a problem for more numerous species where larger sample sizes increase the signal to noise ratio.

In conclusion, comparing genetic structure and demographic decorrelation among species inhabiting the same coastline allowed us to test hypotheses on the probable mechanisms restricting or promoting population connectivity and gene flow in a marine environment. Clearly, a pelagic larval stage is associated with increased dispersal, but perhaps not to the extent that might have been expected. Instead, we find that population connectivity in all three species is restricted and that gaps in suitable habitat may impose barriers to effective dispersal in coastal species in general. The present study thus corroborates recent findings emphasizing that life‐history influences genetic patterns of marine coastal fish species (Barry et al., 2021) and illustrates the power of a comparative multispecies approach.

CONFLICT OF INTERESTS

We declare we have no competing interests.

AUTHOR CONTRIBUTIONS

HK, PEJ, LR, MS, KJ, MMH and CA designed the study; HK, PEJ, MS, MM and IM collected the data; HK, PEJ, DC, LR, MM, MS, SHE and MJ analysed the data; HK, PEJ, DC, LR, SJ and MS wrote the manuscript; MM, MJ, JAH, KJ, OR, MMH, CA and SH contributed to the writing of the manuscript and HK, PEJ, MM and LR prepared the figures.

Fig S1‐S4

Click here for additional data file.