Reconciling seascape genetics and fisheries science in three co-distributed flatfishes

Uncertainty hampers innovative mixed-fisheries management by the scales at which connectivity dynamics are relevant to management objectives. The spatial scale of sustainable stock management is species-specific and depends on ecology, life history and population connectivity. One valuable approach to understand these spatial scales is to determine to what extent population genetic structure correlates with the oceanographic environment. Here we compare the level of genetic connectivity in three co-distributed and commercially exploited demersal flatfish species living in the North East Atlantic Ocean. Population genetic structure was analysed based on 14, 14 and 10 neutral DNA microsatellite markers for turbot, brill and sole respectively. We then used redundancy analysis (RDA) to attribute the genetic variation to spatial (geographic location), temporal (sampling year) and oceanographic (water column characteristics) components.Genotypes were analysed in order to compare levels of genetic variation and genetic differentiation between the three species. Multi-locus genotypes were tested for deviations from Hardy-Weinberg equilibrium using the pegas package in the R software (Paradis, 2010; R Core Team, 2020). Linkage disequilibrium was evaluated using Fisher's exact test implemented in the genepop package in R (Rousset, 2008). R package hierfstat was used to test for significance of FIS (reflecting heterozygote deficiency/excess) using a randomization test (Goudet & Jombart, 2015). Subsequently the level of genetic variation for each sample was estimated as number of alleles (allelic richness), observed (HO) and expected (HE) heterozygosity. We evaluated the proportional importance of geographical location (SPACE), sampling year (TIME) and water-column dynamics (ENV) in explaining genetic connectivity patterns. To do so, the genotype matrix of each species was first converted into allele counts, where each row is an individual and each column indicates the count (0, 1 or 2) per allele. Redundancy analysis (RDA) identified the spatio-temporal and environmental features explaining genetic (i.e. allelic) variation among the individuals of each species. The significance of RDA models was determined using 1000 random permutations with the vegan package v2.5.6 in R (Oksanen et al., 2019). The RDA Significant models were then subjected to forward selection using the ordiR2step function implemented in the vegan package, including a threshold of α = 0.05 and given the adjusted R² parameter of the RDA with all variables included to obtain an unbiased selection (Blanchet, Legendre, & Borcard, 2008). Forward selection corrects for highly inflated type I errors and overestimated amounts of explained variation. The reduced set of explanatory variables based on forward selection, was then used to recalculate the explained proportion of genetic variation. Geographic variables (SPACE) were represented by Moran's Eigenvector Maps (MEMs), along with longitude and latitude. The MEMs were calculated for each individual species, based on a distance based matrix of the shortest geographic waterway distance between sampling locations (Borcard & Legendre, 2002). Temporal variables (TIME) were represented by dummy variables from sampling year indicators. The year a sample was obtained was scored with the value 1; other years were marked with the value 0. Lastly, water-column variables (ENV) for the greater North Sea area (including English Channel and Skagerrak) were downloaded from the ICES WGOOFE website (groupsites.ices.dk/sites/wgoofe). For each of the nine water-column parameters, the monthly and yearly average were available for the time period 1980-2004. The yearly standard deviation of each variable across the 12 months was calculated to capture the intra-annual variation. Seascape genetic analyses were conducted with two sets of environmental variables (ENV1 and ENV2). The first set (ENV1) consisted of 18 variables including the yearly average and standard deviation of each of the nine parameters. This set covers the broad environmental variation and captures the relevant biological information for a comprehensive analysis across the three species in an identical dynamic system. For the second set (ENV2), we selected specific month averages for each species. Specifically, we selected April, May, June, and September for turbot (4 months; 36 variables), March, May, June, and September for brill (4 months; 36 variables), and February, April, and September for sole (3 months; 27 variables). This set allowed us to test whether specific seasonal variation in reproduction time (start and peak spawning time) affects the genetic variation among individuals