Distribution and population structure in the naked goby Gobiosoma bosc (Perciformes: Gobiidae) along a salinity gradient in two western Atlantic estuaries

Many species of fish produce larvae that undergo a prolonged dispersal phase. However, evidence from a number of recent studies on demersal fishes suggests that the dispersal of propagules may not be strongly correlated with gene flow. Instead, other factors like larval behavior and the availability of preferred settlement habitat may be more important to maintaining population structure. We used an ecologically important benthic fish species, Gobiosoma bosc (naked goby), to investigate local and regional scale population structure and gene flow along a salinity gradient (∼3 ppt to ∼18 ppt) in two North Carolina estuaries. G. bosc is an abundant and geographically widespread species that requires complex but patchy microhabitat (e.g. oyster reefs, rubble, woody debris) for reproduction and refuge. We sequenced 155 fish from 10 sites, using a common barcoding gene (COI). We also included recent sequence data from GenBank to determine how North Carolina populations fit into the larger biogeographic understanding of this species. In North Carolina, we found a significant amount of gene flow within and between estuaries. Our analysis also showed high predicted genetic diversity based upon a large number of rare haplotypes found within many of our sampled populations. Moreover, we detected a number of new haplotypes in North Carolina that had not yet been observed in prior work. Sampling along a salinity gradient did not reveal any significant positive or negative correlations between salinity and genetic diversity, nor the proportion of singleton haplotypes, with the exception of a positive correlation between salinity standard deviation and genetic diversity. We also found evidence that an introduced European population of naked gobies may have originated from an Atlantic source population. Altogether, this system offers a compelling way to evaluate whether factors other than dispersal per se mediate recruitment in an estuarine-dependent species of fish with a larval dispersal phase. It also demonstrates the importance of exploring both smaller and larger scale population structure in marine organisms to better understand local and regional patterns of population connectivity and gene flow.


Manuscript to be reviewed
To quantify population structure in these estuaries, we used the common mitochondrial 150 barcoding gene Cytochrome Oxidase I, which has been used as a tool for population genetics 151 studies across multiple diverse taxa over the past 20+ years (COI 2014; Barcode 2018), and has 152 proven effective at sequencing the genetic diversity of marine fishes (e.g. Ward et al. 2005). 153 Furthermore, the use of this marker allowed us to explore our data in combination with other  Two replicate crates were deployed at each sample location, making for a total of ten collecting 196 units on the Neuse River and eight on the Pamlico River. Crates were zip-tied to 0.75 m wooden 197 stakes secured in the nearshore subtidal zone, or deployed from fixed or floating docks using 198 rope. Crates on both rivers were checked every 6-8 weeks, and the contents sorted using a large 199 sieve (55.9 x 55.9 x 12.7 cm) with 2 mm mesh. A maximum of ten naked gobies were collected 200 from each sample site during each sampling event. In order to minimize selection bias, fish from 201 both crates were pooled together, and ten individuals were randomly selected from a grid divided 202 into four quadrants. Only sexually mature adults were used in this study.    In other words, this popset contained no information on haplotype frequencies per population.

247
On the other hand, our North Carolina dataset included frequency data that were explored   We also used ARLEQUIN to obtain genetic diversity values for each population to investigate 266 whether there was any effect of salinity on haplotype diversity in the North Carolina populations.

267
For this salinity analysis, we included sites that were sampled from the two rivers (n=9 sites: 4 in 268 the Pamlico River and 5 in the Neuse River). We used JMP Pro 13 (SAS Institute Inc., Cary, (calculated from the haplotype analysis described above). This latter approach was used to 271 determine whether there was any influence of salinity on the proportion of rare haplotypes in a 272 population. Additionally, given the variability in salinity within these systems, we also explored 273 whether there was any relationship between the standard deviation of salinity and genetic 274 diversity.

276
In addition, Primer 7.0.13 was used to construct rarefaction and extrapolation curves of

289
For our second analysis focused on goby biogeography, we performed an AMOVA to explore 290 differences at the subregional level and also at the larger regional level, specifically between the 291 Atlantic and Gulf of Mexico regions. We again used PopArt to graphically create haplotype  in naked gobies using a major estuarine system as our focal region. Our study also contributes 384 previously undocumented sequence data from North Carolina populations of G. bosc, which will 385 help in further resolving questions related to gene flow in this species at both the local and 386 biogeographic scales. In addition, we incorporated a salinity gradient into our study design, as 387 individuals would need to be sampled to produce an asymptote; in turn, greater sampling would 408 produce significantly more haplotypes in each estuary than were initially detected.

410
In addition to demonstrating high diversity and a large number of singleton haplotypes (just a   (Fig. 3), Pamlico sites are closer together spatially than Neuse sites. The location of our 435 sample sites along each river, as well as the topography of the two rivers themselves, could be 436 influential in differentially affecting gene flow in these rives. Most of our sites were located on 437 creeks that were tributary to the main stem of either river. Although these sites were all 438 positioned within 1 km of the river, some creeks would be subject to more flushing than others, 439 while more hydrodynamically-isolated areas would favor greater larval retention. In addition, the 440 orientation of the Pamlico itself is straighter (Fig. 1) and thus may be more favorable to wind-  Though COI is generally treated as a neutral marker (i.e., it is widely used as a "barcoding"      Table 1 and Table S1 for more information). Base map layers