The impact of aquaculture on the genetics and distribution of the onuphid annelid Diopatra biscayensis

Abstract Aim Evolutionary history of natural populations can be confounded by human intervention such as the case of decorator worm species Diopatra (Onuphidae), which have a history of being transported through anthropogenic activities. Because they build tubes and act as ecosystem engineers, they can have a large impact on the overall ecosystem in which they occur. One conspicuous member, Diopatra biscayensis, which was only described in 2012, has a fragmented distribution that includes the Bay of Biscay and the Normanno‐Breton Gulf in the English Channel. This study explores the origin of these worms in the Normanno‐Breton region, which has been debated to either be the result of a historic range contraction from a relic continuous population or a more recent introduction. Location Northeastern Atlantic, the Bay of Biscay, and the Normanno‐Breton Gulf. Methods We utilized a RAD‐tag‐based SNP approach to create a reduced genomic data set to recover fine‐scale population structure and infer which hypothesis best describes the D. biscayensis biogeographic distribution. The reduced genomic data set was used to calculate standard genetic diversities and genetic differentiation statistics, and utilized various clustering analyses, including PCAs, DAPC, and admixture. Results Clustering analyses were consistent with D. biscayensis as a single population spanning the Bay of Biscay to the Normanno‐Breton Gulf in the English Channel, although unexpected genetic substructure was recovered from Arcachon Bay, in the middle of its geographic range. Consistent with a hypothesized introduction, the isolated Sainte‐Anne locality in the Normanno‐Breton Gulf was recovered to be a subset of the diversity found in the rest of the Bay of Biscay. Main conclusions These results are congruent with previous simulations that did not support connectivity from the Bay of Biscay to the Normanno‐Breton Gulf by natural dispersal. These genomic findings, with support from previous climatic studies, further support the hypothesis that D. biscayensis phylogeographic connectivity is the result of introductions, likely through the regions’ rich shellfish aquaculture, and not of a historically held range contraction.


| INTRODUC TI ON
Anthropogenic activities can make assessing historic organismal distributions challenging. Within short time intervals, humans can transport individuals across broad geographic regions, cause extinction of populations, or facilitate range shifts (Strauss et al., 2006) impacting the ability to assess population-level dynamics. Distinguishing between changes over time due to natural causes or due to humanmediated activities is important to the evaluation of anthropogenic environmental impacts and to aid management and conservation practices. One such example where human activities appear to have obscured biogeographic history includes decorator worm species belonging to Diopatra (Audouin & Milne Edwards, 1833), Onuphidae.
Large species of Diopatra have present-day phylogeographic patterns that are particularly difficult to interpret as they have been moved for sale as fishing bait, in association with shellfish aquaculture, and have undergone natural shifts in ranges partly in response to climate change (Arias et al., 2016;Berke et al., 2010;Saito et al., 2014;van der Have et al., 2015;Wethey et al., 2016;Woodin et al., 2014). The oldest historical records of D. biscayensis are from collections from Spanish sites in the early 1900s. These were misidentified as D. neapolitana at the time of collection (Arias & Paxton, 2015). Up until at least the 1930s, and probably into the 1960s, neither species was known north of Arcachon, France (44.66°N, 1.1417°W) (Berke et al., 2010;Fauchald et al., 2012;Faure, 1969). Diopatra biscayensis produces short-lived gelatinous egg masses attached to its tube, which disintegrate after ~ 2 days, releasing lecithotrophic larvae that rapidly settle (Arias & Paxton, 2015). Based on larval transport simulations, the larvae are unlikely to disperse more than 50 km and probably much less (Woodin et al., 2014). However, most of the disjunct localities within the Normanno-Breton Gulf in the English Channel (2.4°W-1.6°W, 48.6°N-48.7°N) are more than 450 km away from other known localities and are associated with extensive imports of aquaculture material (Goulletquer & Le Moine, 2002;Woodin et al., 2014). Further, the 450 km of coastline between disjunct localities, including the Brittany Peninsula, was shown in previous simulations to not be hospitable to Diopatra until likely the late 21st century due to cold water temperatures (Wethey et al., 2016).
Important source localities for aquaculture are in the region between the Loire and Gironde rivers on the central French Biscay coast (Héral, 1989;Muehlbauer et al., 2014), areas with very large populations of D. biscayensis (Woodin et al., 2014).
Whereas Woodin et al. (2014) hypothesized human-assisted transport seemed likely given the lack of evidence of stepping stone dispersal and limited larval transport probability, Arias and Paxton (2015)  Currently available genetic and morphological evidence (Arias et al., 2016;Berke et al., 2010;Pires et al., 2010) is consistent with the human-assisted transport hypothesis of Woodin et al. (2014) and Wethey et al. (2016) as is the common finding of juvenile Diopatra that have settled on clam boxes, ropes for seed collection, and living oyster and mussel shells (unpub.data, A. Arias and S.A. Woodin). However, available genetic data are based solely on mitochondrial, single-locus markers, which have limited resolution over the timescale of interest. Thus, to distinguish between the human-mediated transport hypothesis and the range reduction hypothesis (Arias & Paxton, 2015), we employed a reduced representation genome approach, specifically the 2b-RAD method (Wang et al., 2012), to assess evolutionary patterns of population genetic variation. Notably, the two competing hypotheses differ in their predictions of isolation by distance and genetic similarity when comparing northern localities to those further south. If the disjunct northern worms represent a relic of a formerly continuous population, one would expect the genetic signature to reflect that isolation given the multiple generations necessary to traverse over 450 km of coastline given a very short larval period. Alternatively, if the disjunct northern locality is the result of a recent introduction via human-assisted transport, then the expectation is of a genetic signature congruent with, and a subset of, the haplotype diversity recovered from the contiguous population of the Bay of Biscay. Additionally, genetic diversity indices observed heterozygosity (H o ), expected heterozygosity (H e ), inbreeding coefficient (F is ), and the lower limit (ll) and upper limit (ul) confidence intervals calculated from 1,000 bootstrap replicates are also included.

| Sample collection
Only two large species of Diopatra are known from the Atlantic coast of Europe. Diopatra biscayensis is an onuphid polychaete that builds a very conspicuous tube that extends above the sediment surface (i.e., "tube-cap"), and it is heavily decorated with debris such as shell fragments, and the aperture typically points downward (Wethey et al., 2016). Its congener D. neapolitana, with a similar but more southerly range, has a tube with a more limited extension above the sediment surface, minimal decoration, much thicker interior tube lining, and upward-pointing tube aperture (Arias et al., 2016).
Approximately 20-30 individuals of D. biscayensis were taken by shovel from each collection locality by the authors, except those from Santander, which were obtained by AA from bait sellers who had collected them that morning ( Figure 1, Table 1  HiSeq using 50 bp single-end chemistry.

| Genomic data
Raw sequence reads were demultiplexed and quality-filtered to remove any reads with less than 90% of nucleotides having a Phred score above 20. Additionally, reads were filtered for the presence of the AlfI recognition site (scripts from https://github.com/Eli-Meyer/ 2brad_utili ties). Retained filtered sequences were then assembled using the Stacks v 2.1 (Catchen et al., 2011) denovo_map.pl pipeline with most parameters left in the default setting (except for distances allowed between stacks, M = 3, and distance allowed between catalog loci, n = 3, based on recommendations by Paris et al. (2017)).
Retained loci were then processed in the populations program of the Stacks package with a minimum minor allele frequency of 0.05, a maximum observed heterozygosity of 0.5, and retaining only one SNP per RAD locus with the write_single_snp command. To reduce the amount of missing data, a standard filtering approach (Benestan et al., 2015;Galaska et al., 2017b;Paris et al., 2017) was applied.
A given locus had to be present in at least 70% of individuals at a given sampling locality and in five of the seven localities to be retained. This filtering scheme was selected to remove any loci that did not have adequate sequencing coverage across the sample set and to mitigate the effects that allele dropout can have on estimating variation within and between populations (Gautier et al., 2013).
The resulting SNP matrix was then exported in "Structure" format for downstream analyses based on the allelic variation. Analyses that required different data formats were converted through the program PGDSpider v2.0.8.1 (Lischer & Excoffier, 2012).

| Genetic diversity, population structure, and gene flow
The resulting SNP matrix was imported into the R v3.4.2 statistical environment (R Core Team, 2017) for further filtering and population genetic inference. BayeScan v.2.01 (Foll & Gaggiotti, 2008) was used with a 50,000 interval burn-in, 100,000 intervals, and priors of 1,000 and 10,000 to identify any loci potentially under selection. Identification of putative loci under possible selection was attempted using a blastn analysis against the NCBI nucleotide database.
Tests of standard genetic distances and general nucleotide diversity indices were also calculated in the R statistical environment v 3.4.2 (R Core Team, 2017) with the package hierfstat v0.04-22 (Goudet, 2005). The statistic "Dch" or the Cavalli-Sforza and Edwards Chord distance (Cavalli-Sforza & Edwards, 1967) was selected for its increased accuracy in estimating relationships among samples, in comparison with other measures of genetic distance (Takezaki & Nei, 1996). A pairwise F st estimate using Weir and Cockerham (1984) statistic was also calculated and then bootstrapped for 1,000 replicates to provide upper and lower limit confidence intervals. To calculate signals of isolation by distance (IBD), the Mantel tests were performed in the ade4 v1.7-13 (Dray & Dufour, 2007) software package using the Cavalli-Sforza and Edwards Chord genetic distance values calculated between localities. Analyses of molecular variance (AMOVAs) were calculated in the poppr v2.8.1 software package (Kamvar et al., 2015).
All analyses were done after filtering of loci and individuals.
The number of individuals in Table 1 reflects the 209 available for analysis of the 229 individuals collected. Several specimens were lost due to degraded DNA before the RADseq prep, a few due to low sequencing and several due to low coverage of loci postfiltering.
All such individuals were removed from the entire analysis pipeline so that individuals with missing data did not impact downstream filtering or further analyses.

| Genomic data set generation and filtering
The data set consisted of 209 individuals coded for 4,078 polymorphic SNP loci with an average sequencing depth of 51.4x, after filtering all RAD loci and samples for quality and missing data. BayeScan analyses using priors of 1,000 and 10,000 initially recovered seven putative loci under selection and after removing those loci (which were found throughout sample localities), a secondary run at both priors recovered no additional loci potentially under selection ( Figure   S1). Further, these loci were blasted against NCBI's database but no significant matches were found, likely due to the short length of

| Genetic diversity, population structure, and gene flow
Observed and expected heterozygosity, along with inbreeding coefficients are reported in   (Table 3) recovered the largest variance from individual samples within a locality as expected for panmictic populations, while differences among localities represented the least amount of variance, though both were highly significant. These tests were performed on the entire data set and also with Sainte-Anne removed, as Sainte-Anne is hypothesized to be a subset from the Bay of Biscay, but results only differed slightly between the two analyses. In both cases, differences among localities explained the least amount of the variance, 1.3% in both instances, consistent with little differentiation across the range.
Both DAPC and Admixture analyses, which try to estimate the number of populations (K), found that over multiple simulations K = 1 had the highest likelihood, with K = 2 having the second highest support. In conjunction with the analyses of genetic distances, DAPC and Admixture analyses were influenced by the genetically

Estimates of migration rates and directionality calculations with
Gst done in divMigrate further support these findings (Figure 3).
With a moderate Gst cutoff of 0.7, we recovered no meaningful connections of Arcachon to its surrounding localities. Migration analyses were subsequently run without Arcachon to mitigate any impact on the inference between the remaining localities, but results did not differ. According to divMigrate, Le Bile, the furthest north locality in the Bay of Biscay, only showed connection with San Vicente, and Sainte-Anne was most strongly connected with Tharon Plage.
However, given that DAPC analysis, Admixture analysis, and PCA all show limited genetic structure between localities, we are reserved about the biological meaning of the divMigrate results. Calculating migration between genetically similar localities, with potential human-mediated transport, is problematic.
All tests of IBD were found to not be significant, with most recovering a negative correlation. Given that Arcachon was a geographically central but genetically distinct locality, tests of IBD were performed with and without these individuals. The recovered IBD value was still negative at −0.075 and found to not be significant, with a p-value of 0.504 (Figure 4).

F I G U R E 4 Mantel test of isolation by distance based on Cavalli-Sforza and Edwards
Chord distances between all seven sampling localities of Diopatra biscayensis. Tests were found to not be significant, suggesting that geographic distance over the sampling range is not a good predictor of genetic differentiation Gulf is many times greater than the planktonic larval duration of this worm (Wethey et al., 2016). Sainte-Anne also had the most limited diversity and represents a subset of the diversity recovered from the rest of the geographic range of D. biscayensis (PCA convex hull analyses), indicative of an introduction rather than an historic range contraction (Figure 2). An introduction of few founders followed by reproduction among those individuals would yield this result, similar to that seen in bays with limited water exchange (i.e., long water time residence), and thus facilitated autorecruitment (Ayata et al., 2009;Plus et al., 2009).
The two Spanish collection sites (Santander and San Vicente de la Barquera) retained the largest window space in the convex hull analysis (Figure 2b), indicating the highest genetic diversity among all studied localities. Interestingly, these localities correspond to the first records of D. biscayensis in Europe, dating back to the early 1900s (Arias & Paxton, 2015). Both are also locations with natural oyster banks that supported a considerable industry at that time (Torres & Corral, 2001). Given the long history of translocation of shellfish among regions in Europe (e.g., Héral, 1989;Muehlbauer et al., 2014), and given the existence of an oyster fishery in San Vicente de la Barquera at least since the 1790s (Graells, 1870), an artisanal fishery in Santander in the 1860s (Graells, 1870), and a commercial fishery there in the 1890s (Dean, 1891), the development of commercial oyster culture in the Basque Country in the 1860s (Graells, 1870), and the use of oyster culturing methods from Arcachon in Spain (Balaguer y Primo, 1878), D. biscayensis plausibly dispersed from northern Spain with oyster transfers.
Migration analyses (Figure 3) suggested the geographically isolated Sainte-Anne locality is more strongly connected to the Tharon Plage locality than to other sampling sites, but note that all connection values are >0.7 with the exception of Arcachon. As this relationship is unlikely through natural stepping stone dispersal across 450 km of inhospitably cold waters for D. biscayensis' reproduction (Wethey et al., 2016;Woodin et al., 2014), we considered possible anthropogenic vectors for this connectivity. Biofouling on ships is a common dispersal vector for many marine species (Floerl & Coutts, 2009), but the aquaculture industry remains a primary candidate for nontarget species dispersal through hitchhiking. Examples of polychaetes that overcame biogeographic barriers by humanmediated transport in mollusk aquaculture are common (Naylor et al., 2001;Simon et al., 2009;Williams et al., 2016). Mussel seed transport to Northern Brittany and the Bay of Mont-Saint-Michel originated from the Bay of Biscay in 1954 and 1965, respectively, and continues to this day (Goulletquer & Le Moine, 2002). Important potential source populations are in the region between the Loire and Gironde rivers on the central French Biscay coast (Héral, 1989;Muehlbauer et al., 2014), which include very large populations of D. biscayensis. Ropes are put out into the intertidal for collection of seed and then transported by truck to grow-out locations in the Bay of Mont-Saint-Michel and further north within the Bay of Biscay. Transport of oyster seed and adults follows similar pathways (Buestel et al., 2009;Goulletquer & Heral, 1997;Goulletquer, 1998;Goulletquer & Le Moine, 2002). Aquaculture materials are the likely transport mechanism given the finding of juvenile Diopatra on live mussel and oyster shells, ropes, etc. The long history of mussel and oyster aquaculture in France associated with these aquaculture practices makes this the likely vector for how D. biscayensis colonized the northern region of the species geographic range as suggested by Woodin et al. (2014). In the same way, another largesized errant polychaete Marphysa victori, originally described from Arcachon Bay in 2017 (Lavesque et al., 2017), is now considered to be an alien species in France (Lavesque et al., 2020). Based on morphological and molecular evidence, these authors have proposed an Asiatic origin of the species, considering that it was introduced from China or Japan into Arcachon Bay in the 1970s with the nonnative oyster Magallana (Crassostrea) gigas. Consistent with this notion, there is other abundant evidence for the association between aquaculture and the introduction of exotic species into Europe and secondary introductions among areas with aquaculture (Goulletquer & Le Moine, 2002;Mineur et al., 2007).
Although our interests originated with the connectivity of Sainte-Anne locality to the rest of the Bay of Biscay, the genetic distinctness of Arcachon necessitates further investigation. This bay, while centrally located, shows signs of genetic differentiation from the rest of the geographic range of D. biscayensis and migration analyses did not uncover even moderate connectivity to the surrounding localities. The recovered substructure of Arcachon within the singular population is unusual given its central placement in the Bay of Biscay. Arcachon is a semi-closed bay with aquaculture efforts focused on oyster farming and is a popular tourist location, increasing the potential for anthropogenic transport by the live bait industry. One possible explanation for the limited connectivity with surrounding localities is the unique hydrodynamics of Arcachon Bay. Specifically, river flows and wind-driven currents have only minor impacts on the hydrology, instead tidal flows account for the majority of water mass movement, which are somewhat restricted by a sill near the mouth at ~ 20 m depth; most of the water masses re-enter the bay after moving seaward on ebb tides (Plus et al., 2009). Plus et al. (2009) found that under ideal conditions with strong northerly and westerly winds, the flushing time for the bay ranged from 13.3 to 15.9 days; given that D. biscayensis larvae are in the water column for 4-5 days, this likely limits the species' capabilities to disperse in or out of the bay. Strong currents have been shown to influence the recovered phylogeography with multiple other marine species (Collins et al., 2018;Galaska et al., 2017a;Xuereb et al., 2018). Some larval distribution data support the idea of domination by tidal advection in Arcachon Bay (Marcano & Cazaux, 1994;Mathivat-Lallier & Cazaux, 1990). If currents are driving the genetic structure recovered at Arcachon, other species in the bay that have limited dispersal capabilities would also share similar phylogeographic patterns.
The single population of D. biscayensis recovered from the Bay of Biscay up to the Normanno-Breton Gulf serves as another example of how aquaculture can affect natural populations and influence ecosystems. Introduced species often can have major impacts on the local ecosystem, causing competition for local resources, interrupting food webs, and leading to an economic impact (Leung et al., 2002). Even estimating the number of introductions can become problematic once the species has become established (Resh et al., 2018). Diopatra biscayensis is an ecosystem engineer that influences the overall habitat which it occupies and can have a major impact on the biodiversity within the region (Berke et al., 2010). Although habitat modification by Diopatra species can increase overall biodiversity for a locality, it also can result in displacement of endemic fauna via competition (Berke et al., 2010). Interestingly, the biology of D. biscayensis is strongly influenced by temperature (Berke et al., 2010;Wethey et al., 2011;Wethey et al., 2016) and the presumably introduced population in Sainte-Anne is likely to expand given the continued progression of climate change.

ACK N OWLED G M ENTS
We would like to thank Viktoria Bogantes for secondarily identifying

CO N FLI C T O F I NTE R E S T
None declared.

Specimen vouchers have been deposited at the Auburn University
Museum of Natural History (Accession Numbers AUMNH 45517-45683). Raw sequences have been submitted to NCBI's SRA database as BioProject PRJNA706051, Accession Nos. SAMN18115220-SAMN18115428. The final SNP data set has been deposited at Dryad under project code https://doi.org/10.5061/ dryad.4j0zp c8b4.