Rare, long-distance dispersal underpins genetic connectivity in the pink sea fan, Eunicella verrucosa

Characterising patterns of genetic connectivity in marine species is of critical importance given the anthropogenic pressures placed on the marine environment. For sessile species, population connectivity can be shaped by many processes, such as pelagic larval duration, oceanographic boundaries, and currents. This study combines restriction-site associated DNA sequencing (RADseq) and passive particle dispersal modelling to delineate patterns of population connectivity in the pink sea fan, Eunicella verrucosa, a temperate octocoral. Individuals were sampled from 20 sites covering most of the species’ northeast Atlantic range, and a site in the northwest Mediterranean Sea to inform on connectivity across the Atlantic-Mediterranean transition. Using 7,510 neutral SNPs, a geographic cline of genetic clusters was detected, partitioning into: Ireland, Britain, France, Spain (Atlantic), and Portugal and Spain (Mediterranean). Evidence of significant inbreeding was detected at all sites, a finding not detected in a previous study of this species based on microsatellite loci. Genetic connectivity was characterised by an isolation by distance pattern (IBD) (r2 = 0.78, p<0.001), which persisted across the Mediterranean-Atlantic boundary. In contrast, exploration of ancestral population assignment using the program ADMIXTURE indicated genetic partitioning across the Bay of Biscay, which we suggest represents a natural break in the species’ range, possibly linked to a lack of suitable habitat. As the pelagic larval duration (PLD) is unknown, passive particle dispersal simulations were run for 14 and 21 days. For both modelled PLDs, inter-annual variations in particle trajectories suggested that in a long-lived, sessile species, range-wide IBD is driven by rare, longer dispersal events which act to maintain gene flow. These results suggest that oceanographic patterns may facilitate range-wide stepping-stone genetic connectivity in E. verrucosa, and highlight that both oceanography and natural breaks in a species’ range should be considered in the designation of ecologically coherent MPA networks.

Methods
Pink sea fan tissue samples (n = 285) were collected via SCUBA (under sampling license L/2019/00143 granted by the Marine Management Organisation) from depths ranging between 5 and 35 m between 2007-2019. In brief, apical 5 cm cuttings were removed from each sea fan colony underwater and preserved in >95% ethanol immediately after surfacing to prevent DNA degradation. Samples were stored long-term at 4 oC. Genomic DNA was extracted from polyp tissue using a salting-out protocol optimised for gorgonin-protein tissue (see Supplementary Material associated with the manuscript for full protocol).

Sequencing and SNP isolation were performed through nextRAD, a reduced representation sequencing method. Genomic DNA was first fragmented with Nextera reagent (Illumina, Inc), which also ligates short adapter sequences to the ends of the fragments. The Nextera reaction was scaled for fragmenting 5 ng of genomic DNA, although 7.5 ng of genomic DNA was used for input to compensate for degraded DNA in the samples and to increase fragment sizes. Fragmented DNA was then amplified for 26 cycles at 73 oC, with one of the primers matching the adapter and extending 9 nucleotides into the genomic DNA with the selective sequence GTGTAGAGG. Thus, only fragments starting with a sequence that can be hybridized by the selective sequence of the primer will be efficiently amplified. The nextRAD libraries were sequenced on an Illumina HiSeq 4000 with one lane of 2x150 bp reads (University of Oregon)

Raw reads were cleaned, de-multiplexed and trimmed to 140 bp using the process_radtags pipeline in Stacks v2.54 (Rochette and Catchen, 2017) and aligned to the pink sea fan reference genome (Macleod et al., 2023) using Bwa-mem2 v2.2.1 (Vasimuddin et al., 2019). Loci were built using gstacks from Stacks with a minimum stack depth of 3 (Paris et al., 2017).

SNPs were filtered using the populations pipeline of Stacks was used to retain SNPs with a minimum stack depth of 9, present in at least 16 sampling locations, minor allele frequency of 0.05, and maximum loci heterozygosity of 0.6. Further filtering using VCFtools (Danecek et al., 2011) removed genotypes containing less than 90% of SNPs, which were identified with the missingno function in R (R Core Team, 2023) using the package Poppr v2.8.3 (Kamvar et al., 2014), and a depth coverage <9. Loci in linkage disequilibrium were filtered in VCFtools (Danecek et al., 2011) using a correlation (r2) threshold of >0.7 between all marker pairs.  Finally, the function mlg from the R package Poppr was used to identify potential clones. 

After initial filtering in populations, 23,111 SNPs were identified in 285 individuals. After filtering out individuals and loci with greater than 10% missing data and monomorphic SNPs in VCFtools, 8906 SNPs in 246 individuals were retained. Additionally, 1,256 SNPs were removed from each pair of SNPs identified in linkage disequilibrium. No duplicated genotypes were detected. The final dataset consists 7,650 SNPs in 246 individuals.

In brief, for all 20 sampling sites neutral SNPs were identified using outlier detection programs OutFLANK and pcadapt, only neutral SNPs common to both programs were retained. Genetic diversity was assessed via observed/expected heterozygosity and inbreeding coefficient (Fis). Genetic differentiation was computed using Pairwise FST (Weir and Cockerham) and global FST. Genetic population structure was explored using three approaches: a DAPC (discriminant analysis of principal components), ADMIXTURE v1.3, and snapclust.

To compare oceanographic particle dispersal with patterns of genetic connectivity, passive oceanographic modelling was conducted in R using an Atlantic-Iberian Biscay Irish-Ocean Physics analysis and a Forecast oceanographic model from CMEMS (see Supplementary Material associated with manuscipt for full model details). For each site, raw modelled particle trajectories are contained within a geopackage file (.gpkg). 

Distance-based Moran's eigenvector maps (dbMEMs), asymmetric eigenvector maps (AEMs) were produced on marine least-cost distances and on the dispersal estimates between sampling sites, respectively. Environmental variables describing sea surface and bottom temperature, surface and bottom salinity, surface current velocity, were obtained from Copernicus (see Supplementary Material associated with manuscipt for full details). Redundancy analysis was performed to explore the effect of environmental variables, dbMEMs, and AEMs on modelled trajectories.