Population and seascape genomics of a critically endangered benthic elasmobranch, the blue skate Dipturus batis

Abstract The blue skate (Dipturus batis) has a patchy distribution across the North‐East Atlantic Ocean, largely restricted to occidental seas around the British Isles following fisheries‐induced population declines and extirpations. The viability of remnant populations remains uncertain and could be impacted by continued fishing and by‐catch pressure, and the projected impacts of climate change. We genotyped 503 samples of D. batis, obtained opportunistically from the widest available geographic range, across 6 350 single nucleotide polymorphisms (SNPs) using a reduced‐representation sequencing approach. Genotypes were used to assess the species’ contemporary population structure, estimate effective population sizes and identify putative signals of selection in relation to environmental variables using a seascape genomics approach. We identified genetic discontinuities between inshore (British Isles) and offshore (Rockall and Faroe Island) populations, with differentiation most pronounced across the deep waters of the Rockall Trough. Effective population sizes were largest in the Celtic Sea and Rockall, but low enough to be of potential conservation concern among Scottish and Faroese sites. Among the 21 candidate SNPs under positive selection was one significantly correlated with environmental variables predicted to be affected by climate change, including bottom temperature, salinity and pH. The paucity of well‐annotated elasmobranch genomes precluded us from identifying a putative function for this SNP. Nevertheless, our findings suggest that climate change could inflict a strong selective force upon remnant populations of D. batis, further constraining its already‐restricted habitat. Furthermore, the results provide fundamental insights on the distribution, behaviour and evolutionary biology of D. batis in the North‐East Atlantic that will be useful for the establishment of conservation actions for this and other critically endangered elasmobranchs.


| INTRODUC TI ON
Many elasmobranchs have experienced drastic population declines as a consequence of fishing pressure during the last century, representing a major conservation concern. Almost one-third of elasmobranch species globally are threatened with extinction, yet nearly half remain too data-deficient to be assessed (Dulvy et al., 2014(Dulvy et al., , 2017IUCN, 2021). Despite fishing restrictions, a large number are still caught as by-catch, particularly in unregulated coastal and continental waters (Dulvy et al., 2017) where they can be of significant socio-economic importance to local fisheries ICES, 2020). The K-selected life history that most elasmobranchs exhibit exacerbates the impacts of exploitation; their characteristically slow growth, late-onset maturity and relatively low reproductive output limit population recovery potential (Dulvy et al., 2017).
In addition, evidence is mounting on the consequences of climate change for elasmobranch fitness (Di Santo, 2016;Dziergwa et al., 2019;Pistevos et al., 2015). For many data-deficient elasmobranchs, instituting appropriate conservation actions requires a better understanding of their population structure and of their current and future realized niche in the face of environmental changes.
Elasmobranchs exhibit a range of life history traits that translate to different degrees of population structuring. Some species demonstrate high levels of gene flow across ocean basins (Lieber et al., 2020), while others are divided into smaller subpopulations with limited gene flow (Le Port & Lavery, 2012;Thorburn et al., 2018). A wide range of behaviours such as site fidelity and natal philopatry (Corrigan et al., 2018;Feutry et al., 2017;Pardini et al., 2001;Thorburn et al., 2018), long-distance migrations (Blower et al., 2012;Cameron et al., 2018;Corrigan et al., 2018) and aggregating behaviour among closely related individuals (Lieber et al., 2020;Thorburn et al., 2018) can shape patterns of elasmobranch population connectivity and genetic diversity. In addition, environmental discontinuities such as bathymetric barriers (Le Port & Lavery, 2012) and temperature gradients (Griffiths et al., 2010) can influence species distributions and population connectivity, especially for less vagile species. The diversity and complexity of elasmobranch life histories are likely underappreciated due to issues such as taxonomic confusion (Iglésias et al., 2010) and misreporting of catches (ICES, 2020). Consequently, current conservation strategies that include marine protected areas (MPAs) have been suggested by some as oversimplified and ineffective (Dulvy et al., 2017;Dureuil et al., 2018), requiring more comprehensive species-specific assessments.
Climate change represents a major threat to global biodiversity.
In particular, climatic extremes such as maximum temperatures may lead to higher probabilities of local extinctions for species that are unable to disperse or adapt to these conditions (Román-Palacios & Wiens, 2020). In the North-East Atlantic Ocean, coastal waters are projected to experience temperature rises and acidification, and decreasing dissolved oxygen and salinity levels by the end of the century, while extreme oceanographic events are expected to increase in frequency and magnitude (MCCIP, 2020;Penny Holliday et al., 2020). Rising temperatures have already been associated with poleward distribution shifts for many species (Barton et al., 2016;Brattegard, 2011;Chaudhary et al., 2021;Perry et al., 2005) and have been linked to decreases in individual growth rates and fitness through the pejus effect (Morrongiello & Thresher, 2015).
Behavioural changes have also been documented on a local scale, with some benthic elasmobranchs exploiting deeper thermal refugia (Perry et al., 2005). However, species that are more sedentary in nature may not be capable of undertaking spatial distribution shifts; in these cases, survival may depend upon physiological adaptation to a changing environment. For marine elasmobranchs, the projected environmental changes are likely to incur important physiological costs, particularly in relation to osmoregulation and acid-base regulation to maintain homeostasis. While some elasmobranchs have adapted strategies to cope with environmental extremes (Dziergwa et al., 2019;Heinrich et al., 2014), others are likely to suffer greater losses in individual fitness (Di Santo, 2016;Pistevos et al., 2015).
For nonmodel species that cannot be studied in situ or experimentally, novel molecular approaches in the era of next-generation sequencing (NGS) can provide insights into the structure and local adaptation of wild populations. Ideally, the assembly and annotation of full genomes would provide a functional basis for genomic investigations of a species. However, genome assembly remains prohibitively costly and resource heavy to address urgent conservation questions at the scale of populations, especially given that elasmobranch genomes can be large and complex (Hara et al., 2018).
Reduced-representation sequencing (RRS) methods provide an alternative approach, whereby thousands of genome-wide single nucleotide polymorphisms (SNPs) can be examined in the absence of a reference genome (Andrews et al., 2016). These generate highresolution data to estimate genetic differentiation even when sample sizes are small (Willing et al., 2012), which is an advantage for studies on rare species that rely on opportunistic sampling. Genomewide SNPs can also be used to estimate effective population size (N e ), a theoretical estimator of population size after accounting for genetic drift, which is often used in conservation genetics (reviewed in Allendorf et al., 2013). Furthermore, genotype-environment studies have taken a leap forward with the arrival of NGS methods.
Landscape (or seascape) genomics combines genomic and environmental data to investigate how genetic structuring may be driven by environmental variables, and can reveal candidate genes under selection in certain environmental conditions (Balkenhol et al., 2017;Riginos et al., 2016;Roffler et al., 2016).
In this study, we used a population and seascape genomics approach to investigate patterns of population structuring, abundance and local adaptation in a critically endangered elasmobranch, the blue skate Dipturus batis (Linnaeus, 1758). D. batis has only recently received species status after morphological and genetic investigations distinguished it from the parapatric flapper skate D. intermedius (Griffiths et al., 2010;Iglésias et al., 2010). Both species have become a conservation priority as a result of population declines and range restrictions, but are still caught as by-catch despite an EU landing ban (Ellis et al., 2016). In addition to the continual threat posed by fisheries-induced mortality, the ability of the skates to adapt to environmental changes has become a pertinent question; D. batis currently exploits a narrower thermal niche than D. intermedius (Frost et al., 2020) and consequently may respond differently to ocean warming. Using samples collected from the widest available extent of D. batis' range, obtained through a combination of research surveys and samples of opportunity, we applied a RRS approach (DArTseq TM ; Kilian et al., 2012) to (i) assess the level of gene flow and levels of genetic diversity among extant D. batis populations, (ii) estimate their effective population sizes and (iii) identify potential signals of selection in relation to environmental conditions. Assessment of D. batis' population structure and the potential consequences of environmental change on its ecological niche are necessary to define current and future spatial management units, and may help identify areas that may qualify for additional protection. The threat of continued fishing mortality and uncertain effects of environmental change, combined with the challenge presented by data deficiency, are characteristics that D. batis shares with many endangered elasmobranchs. It is hoped that our primarily molecular approach will effectively address important knowledge gaps for D. batis and that it may be extended to other elasmobranchs of conservation concern.
2 | MATERIAL S AND ME THODS

| Study species
The blue skate Dipturus batis is one of two rajids classified as Critically Endangered by the IUCN (Dulvy et al., 2006) formerly belonging to the common skate species complex (D. batis complex). The two common skate species, which co-occur in parts of their North-East Atlantic range (Frost et al., 2020), were recently differentiated into the smaller-bodied blue skate (D. batis) and the larger-bodied flapper skate (D. batis) based on morphological and genetic differences (Griffiths et al., 2010;Iglésias et al., 2010). At present, the confirmed geographic range of D. batis extends from the Celtic Sea to north of Orkney, with higher densities occurring in the Celtic Sea and Rockall (Frost et al., 2020;Griffiths et al., 2010), while their occurrence has recently been confirmed in Iceland (Bache-Jeffreys et al., 2021).

| Sample collection
Samples of D. batis were primarily obtained from fishery-dependent surveys conducted in the Celtic Sea during the autumn from 2011 to 2017 (Figure S1), and from fishery-independent surveys conducted along the Scottish coast and at Rockall. Faroese samples, obtained opportunistically from fisheries-independent surveys and the commercial fishing vessel 'Sandshavið', were donated by the Faroe Marine Institute. Individuals were identified as D. batis based on the morphological characteristics described by Iglésias et al. (2010), namely body size, eye colour, dorsal patterning and interdorsal space. In addition, 113 of the individuals included in this study (102 from the Celtic Sea, five from Northern Scotland and six from Rockall) were genetically validated as D. batis using diagnostic microsatellite markers by Frost et al. (2020). Sex and length (cm) data were collected from all samples except those from the Faroe Bank.
Fin or muscle tissue samples were collected and stored in 96% ethanol or RNAlater ® . A sample size of 4-6 individuals reportedly provides sufficient power for resolving population genomic structure when over 1000 biallelic markers are used (Willing et al., 2012). To ensure sufficient power in our study, we selected at least 10 individuals from different geographic areas where available, across the narrowest temporal range possible to minimize temporal structuring of our sample set. A total of 564 samples were selected for genomic analysis, from six locations: the Celtic Sea, the Scottish West Coast, Northern Scotland, Rockall, the Faroe Bank and the Faroe Shelf ( Figure 1, Table 1).

| DNA extraction and genotyping
Genomic DNA was extracted using a DNeasy ® Blood & Tissue Kit (Qiagen). DNA concentrations were quantified on a Qubit analytical pipelines. Twelve samples were identified as nontarget species (data not shown) and were removed from the data set before processing of raw sequences was repeated. This generated data for 17,620 sequences of ~69-bp length, each containing a SNP.

| Data filtering
Ten samples were not reported by DArT Pty. Ltd. due to poor sample quality, and an additional sample was removed due to suspected genotyping error, as determined by visually scanning the raw data.
SNPs were further filtered based on a call rate of 80%, and when duplicate loci were present, only the locus with the highest call rate was retained. After this step, the proportion of scored loci per sample was assessed, and all samples were considered to have a sufficiently high score rate to be retained (>88%). Monomorphic loci and loci with low minor allele frequencies (MAF < 0.05) across all samples were identified using adegenet (v 2.1.2, Jombart, 2008;Jombart & Ahmed, 2011), as implemented in R (v 3.6.2, R Core Team, 2019), and subsequently removed. Because human error could lead to sampling an individual multiple times or to contamination during molecular laboratory work, we looked for duplicate samples based on a threshold of 700 mismatching loci (roughly 10% of remaining loci) using the R package CKMRsim (Anderson, https://doi.org/10.5281/ zenodo.820162). Where duplicates were found (i.e. >90% genetically identical), the sample with the highest score rate was retained.
Next, we tested for conformation of loci to the Hardy-Weinberg proportions using the R package pegas (v 0.12, Paradis, 2010), performing an exact test based on Monte Carlo permutation of alleles (Guo & Thompson, 1992) with 1000 replicates for each of the six sampling locations and for the entire dataset (Table S2). After applying the false discovery rate (FDR) correction method of Benjamini and Hochberg (1995), loci were removed if they significantly deviated from Hardy-Weinberg proportions (at a significance threshold of α = 0.05) in at least two sampling locations. We then tested for linkage disequilibrium among loci using the R package snpStats (v 1.36.0; Clayton, 2019) and removed one locus from each pair for which R 2 > 0.80. Following these filtering steps (summarized in Table   S1), the resulting data set contained 503 individuals genotyped at 6,350 loci (Table 1, Figure 1).

| Finding related individuals
Because related pairs of individuals may introduce a bias in population genomic analyses, particularly when sample sizes are small, we looked for first-order (e.g. parent-offspring, full-sibling) and secondorder (e.g. half-sibling) relatives in our data set and removed one in-  (Table 1).

| Population structure
Population-and locus-wide summary statistics were obtained using GenAlEx (v 6.5, Peakall & Smouse, 2006 with the exception of allelic richness, which was estimated using the R package PopGenReport (Adamack & Gruber, 2014). Spatial population structure was assessed using three approaches. First, we employed a Bayesian clustering algorithm using STRUCTURE (v 2.3.4, Pritchard et al., 2000). We used an admixture model with correlated allele frequencies, a burn-in length of 300,000 (more than enough to reach convergence) followed by 500,000 MCMC, and performed five iterations for each prior subpopulation number K (ranging from K = 1 to K = 6). In order to avoid impractically long computation times and potential biases resulting from imbalanced sample sizes (Wang, 2017), we randomly subsampled 10 individuals from the Celtic Sea,

West Coast Scotland and Rockall, and included all samples from
North Scotland (N = 9), the Faroe Shelf (N = 10) and the Faroe Bank (N = 8), such that the total sample size for STRUCTURE analysis was 57. The most likely value of K was estimated using the delta-K method of Evanno et al. (2005) in STRUCTURE Harvester (Earl & VonHoldt, 2012), and summary plots for each K were produced using CLUMPAK (Kopelman et al., 2015). In order to justify the pooling of samples from different years within sample sites, and to ensure that no fine-scale structure would go undetected by subsampling the data for the STRUCTURE analysis, we tested for popula- 1.2% of loci among all samples) were replaced with the mean allele frequencies across all samples. Cluster identification was performed using the find.clusters function, with the optimal number of clusters evaluated using the Bayesian information criterion (BIC). Two DAPC plots were produced: one in which the prior grouping of individuals was based on the evaluated number of clusters, and one with prior groupings based on the six predefined geographic locations.
Third, we performed a principal component analysis (PCA) with the R function prcomp. We used the prefiltered data set of 503 individuals, as the first two principal components from the PCA would later be used as covariates for each individual in our seascape genomics analysis, in which we included all 503 individuals. Since the function does not allow for any missing data, we utilized 3 540 loci with a call rate of 100%.
Overall F-statistics (F IS , F ST ) and pairwise F ST (Weir & Cockerham, 1984) between sampling locations were estimated with the R imple- Rousset, 2008). GenePop was also used to perform overall and pairwise tests of genic differentiation, testing the null hypothesis that all alleles are drawn from the same distribution in all populations. Here, we applied an exact G test (Fisher's method), using 1 000 dememorizations, 100 batches and 1 000 iterations per batch.

| Effective population size
The effective population size (N e ) is a theoretical estimator of population size after accounting for genetic drift that is useful in conservation genetics as it reflects the additive genetic variation, or evolutionary potential, of wild populations (reviewed in Allendorf et al., 2013). We estimated contemporary N e for each sampling location and for each putative population (inferred from the preceding analyses) using the linkage disequilibrium (LD) estimator (Hill, 1981;Waples, 2006;Waples & Do, 2010) in NeEstimator (v 2.1, Do et al., 2014). The estimate assumed random mating and was performed at critical values (i.e. MAF at which alleles should be excluded) of 0.05, 0.02 and 0.01. Confidence intervals were obtained using the jackknife-over-individuals method.

| Seascape genomics and candidate loci under selection
We tested for associations between allele frequencies and en- burn-in length of 50,000) and applied a q-value threshold of 0.05.
Outlier loci from both approaches could therefore be compared.
The sequences containing the SNPs detected in BayeScan and in our seascape genomics analysis were BLASTed on NCBI (NCBI, 1988) to identify functional sequences. The blastn function was used to allow for comparisons across species, given the poorly annotated nature of elasmobranch genomes. We reported the top five BLAST hits with an E-value <0.01.

| Sampled individuals
Sample sizes reflected sampling effort and corresponded to anecdotal reports of D. batis abundance. Most samples were obtained from the Celtic Sea and Rockall during scientific surveys (Table 1).
Males and females were collected at all sites, and a higher number

| Related individuals
The   Figure S2). We identified 10 related pairs, including one full-sibling pair and seven half-sibling pairs from the Celtic Sea, one half-sibling pair from Rockall and one half-sibling pair from the Faroe Bank (Table 2). All related pairs were therefore found at the same locality, and in four of these cases, the pair was collected in the same haul (i.e. at the same time and place).

| Genetic diversity
Patterns of genetic diversity were generally comparable across sites; however, samples from Rockall exhibited lower genetic diversity (

| Spatial population structure
Overall Results of the Bayesian clustering analysis implemented in STRUCTURE suggested that the most likely number of clusters was 2, as per the delta-K method of Evanno et al. (2005). The clustering clearly separated Rockall skates from the rest (Figure 2). When visualizing the output of K = 3, which had only a slightly lower mean log probability than K = 2 ( Figure S7), a third cluster consisting of Faroese skates was identified, with some proportion of assignment to the British cluster for a few samples from the Faroe Shelf ( Figure 2). Similar results were obtained in fastSTRUCTURE ( Figures S3 and S4). In order to assess the influence of loci under selection and the level of potentially adaptive population structure, DAPC was repeated using only the 21 outlier loci identified in BayeScan, and with a neutral data set excluding these 21 loci. Despite weaker clustering that could be expected with the small number of loci, the DAPC still revealed genomic differentiation across the putatively adaptive loci ( Figure S10), while the removal of these 21 loci from the total dataset did not affect the results ( Figure S11). The PCA also showed a clear separation between samples from Rockall and the remaining sites; however, the variation explained by each principal component was low (≤ 1.2%, Figures S12 and S13).

| Effective population sizes
Effective population sizes (N e ) were estimated for each sample site and after grouping samples into the three putative populations.
For the former, sample sizes were too small and N e could not be estimated for some sites (Table S4). For the latter, the British Shelf population had the highest N e (ca. 21,000,

| Seascape genomics and candidate loci under selection
The PCA showed a clear differentiation among sites based on all 34 environmental variables. The first two principal components (PCs) explained 72.7% of the environmental variation among sites ( Figure 4), with all variables contributing to at least one of these PCs ( Figure S14). Overall, there was a clear difference between southern (Celtic Sea) and northern (all other sites) sites across PC1 and PC2, whereas a distinction could also be made between 'offshore' (Celtic Sea and Rockall) and 'inshore' (Scottish and Faroese, also including Faroe Bank) sites across PC2. The variation among sites was complex, but a general pattern could be seen where the Celtic Sea was warmer, more saline and more acidic, while Rockall was the deepest of the sites (up to 385 m deep; Table 6). There was a strong correlation between primary productivity and chlorophyll concentration and between nitrate and phosphate concentrations (Spearman's correlation coefficients >0.9). Removing one variable from each collinear pair did not influence the results of our multiple logistic regressions, and actually reduced the environmental variation explained in the PCA. Therefore, we report results of logistic regression using all 34 environmental variables.
Testing for associations between allele frequencies (two alleles for each of 6,350 loci) and 36 environmental variables (including latitude and longitude as covariates) generated a total of 457,200 tests in Samβada. After the Bonferroni correction, one allele (alternate allele at locus 100069553) was significantly associated with seven environmental variables (p < 0.01 for both G-score and Wald score).
These were as follows: latitude, mean current velocity, maximum pH, minimum bottom temperature, and mean, minimum and maximum salinity. This locus was also detected as one of 21 outliers under putative positive selection in BayeScan (Table 7, Figure S16). On closer  Danio rerio strain Nadia (NA) genome assembly, chromosome: 3 2e-12 97.9 Zebrafish DNA sequence from clone DKEY-106C17 in linkage group 3, complete sequence 2e-12 97.9 Danio rerio strain Nadia (NA) genome assembly, chromosome: 7 5e-12 95.9 Zebrafish DNA sequence from clone CH211-72D16 in linkage group 17, complete sequence Island populations of benthic and benthopelagic fishes (Gonzalez et al., 2014;Johansen et al., 2020;Mattiangeli et al., 2002;Régnier et al., 2017;Saha et al., 2015) and squid (Shaw et al., 1999) to these bathymetric barriers. Unlike these animals, skates are almost exclusively epibenthic throughout their life cycle, laying their egg cases (mermaid's purses) on the sea bed (Last et al., 2016) and spending the majority of their life near the sea floor (Wearmouth & Sims, 2009).
Therefore, the bathymetric barriers are unlikely to be readily overcome unless skates swim long distances high in the water column.
The patterns of population structure observed in D. batis are also likely influenced by the species' highly resident behaviour. The proximity of close relatives identified in this study is indicative of site-attached behaviour, and supports the initial results of tagging in the Celtic Sea, which demonstrates that D. batis almost exclusively remain within relatively confined shallow areas (<200 m) of the continental shelf  and have been recaptured up to 900 km from their release sites (Little, 1998 (~2,300 individuals). The minimum N e required to maintain the evolutionary potential of a population has been a subject of debate, with proposed conservation thresholds set between 500 and 5,000 individuals (reviewed in Allendorf et al., 2013), below which populations may be prone to the loss of genetic diversity, higher inbreeding and the accumulation of deleterious mutations. Whereas the Celtic Sea had relatively high N e , the estimates for the Scottish, Rockall and Faroese populations could be construed as being of conservation concern. We note that N e estimation from high-throughput sequencing data can be subject to bias due to violations to the assumptions of unlinked loci that are difficult to correct for without knowledge of the species' genomic architecture (Waples et al., 2016). Its interpretation also depends on the life history and demography of the species and should ideally be compared with census population size (N c ) (Lieber et al., 2020;Waples et al., 2018). Although these factors may affect the reliability of our N e estimates, they provide an important first estimate of D. batis abundance across its distribution range and serve as initial indicators of their conservation status. The identification of close relatives among our samples suggests that N c could be estimated using close-kin mark-recapture, provided an appropriate survey design .
The seascape genomics analysis revealed a significant association between one locus under putative positive selection (which was identified as an outlier in BayeScan; Table 7) and a number of abiotic variables expected to be influenced by climate change. Previous work has demonstrated that D. batis may be tolerant to a narrower thermal range than its parapatric congener, D. intermedius, where their ranges overlap (Frost et al., 2020); this may lead to differential shifts in fitness under projected climate change. Not only did our findings identify an association between a putatively adaptive locus and minimum bottom temperature, but also associated with pH, salinity and current velocity. A dominance of homozygotes for the reference allele was found in the Celtic Sea (96% of individuals, com- Leucoraja erinacea from the western North Atlantic Ocean, in which individuals from northern parts of the range (colder and less thermally variable) were more sensitive to acidification and warming, as indicated by a lower aerobic capacity and slower recovery following anaerobic activity (Di Santo, 2016). Other studies have demonstrated that warming and acidification could impact elasmobranch hunting ability through impaired olfaction and increased energetic demands (Pistevos et al., 2015), influence juvenile skeletal development (Di Santo, 2019) and lead to denticle corrosion (Dziergwa et al., 2019). The association suggested with mean current velocity could relate to gas exchange rates during respiration when resting on the sea floor, or to the energetic costs of swimming. Due to a lack of reference genomes for D. batis and its relatives, BLAST searches for this and the 20 other putatively adaptive loci did not produce any unambiguous hits. Whereas this locus could be part of a functional sequence, it could also be physically linked to a gene under selection.
Future investigations will benefit from more targeted approaches as genomic resources for elasmobranchs improve.
Our results provide an updated assessment of D. batis' population structure, regional abundance and movement patterns that will have immediate relevance for the conservation of this critically en- demonstrates the efficacy of using increasingly affordable genomewide approaches to address fundamental knowledge gaps for one of many threatened and data-deficient elasmobranchs, leading to better spatial management.

ACK N OWLED G EM ENTS
We are grateful to the fishermen collaborating with CEFAS in the Celtic Sea, Marine Scotland Science and the masters and crew of the MRV Scotia, the crew aboard the Sandshavið and the Faroe Marine Institute (Havstovan) for donating skate tissue samples. We thank Martina Kopp for her assistance in the laboratory, and Andrzej Kilian and his team at Diversity Arrays Technology (DArT Pty. Ltd., Canberra, Australia) for performing the genotyping work. The study was funded by the Nord University.

CO N FLI C T O F I NTE R E S T
The authors declare no conflicts of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available in the Dryad Digital Repository http://doi.org/10.5061/dryad.98sf7 m0kc O RCI D Aurélien Delaval https://orcid.org/0000-0002-2223-2932