A single genomic region involving a putative chromosome rearrangement in flat oyster (Ostrea edulis) is associated with differential host resilience to the parasite Bonamia ostreae

Abstract European flat oyster (Ostrea edulis) is an ecologically and economically important marine bivalve, that has been severely affected by the intracellular parasite Bonamia ostreae. In this study, a flat oyster SNP array (~14,000 SNPs) was used to validate previously reported outlier loci for divergent selection associated with B. ostreae exposure in the Northeast Atlantic Area. A total of 134 wild and hatchery individuals from the North Sea, collected in naïve (NV) and long‐term affected (LTA) areas, were analysed. Genetic diversity and differentiation were related to the sampling origin (wild vs. hatchery) when using neutral markers, and to bonamiosis status (NV vs. LTA) when using outlier loci for divergent selection. Two genetic clusters appeared intermingled in all sampling locations when using outlier loci, and their frequency was associated with their bonamiosis status. When both clusters were compared, outlier data sets showed high genetic divergence (F ST > 0.25) unlike neutral loci (F ST not ≠ 0). Moreover, the cluster associated with LTA samples showed much higher genetic diversity and significant heterozygote excess with outlier loci, but not with neutral data. Most outliers mapped on chromosome 8 (OE‐C8) of the flat oyster genome, supporting a main genomic region underlying resilience to bonamiosis. Furthermore, differentially expressed genes previously reported between NV and LTA strains showed higher mapping density on OE‐C8. A range of relevant immune functions were specifically enriched among genes annotated on OE‐C8, providing hypotheses for resilience mechanisms to an intracellular parasite. The results suggest that marker‐assisted selection could be applied to breed resilient strains of O. edulis to bonamiosis, if lower parasite load and/or higher viability of the LTA genetic cluster following B. ostreae infection is demonstrated.


| INTRODUC TI ON
The European flat oyster (Ostrea edulis L.) is a marine bivalve mollusc distributed from the Atlantic coast of Morocco to the Norwegian Sea and through the Mediterranean Sea up to the Black Sea (Perry et al., 2017). Due to their economic value, flat oyster populations exist outside their natural range in many different countries following human-mediated transfer (Colsoul et al., 2021). Ostrea edulis is a sessile filter-planktotrophic species (Ezgeta-Balić et al., 2020;Peharda et al., 2012) that can reach 20 cm and live up to 20 years, preferably in firm benthic habitats comprised of mud, rocks, hard silt or even artificial substrates created with broken shells (cultch), appearing at depths up to 80 m. The flat oyster is a protandrous hermaphrodite, maturing from 8 to 10 months, depending on environmental factors. Reproduction takes place within the females' pallial cavity, where ova are fertilized by sperm that pass through the gills as part of the normal feeding process. After an incubation period between eight and 10 days, free-swimming larvae are dispersed into the environment before metamorphosis and settlement (Bayne, 2017). Dispersal can last 10-30 days and involve distances longer than 10 km under favourable conditions (Berghahn & Ruth, 2005). Each generation tends to settle on top of the preceding one, such that reefs grow vertically and provide a habitat for other organisms (Chambers et al., 2017). This property is why oysters are one of a number of bivalve species that are considered 'ecosystem engineers' (Smaal et al., 2018).
Flat oyster has been considered a biomonitor to detect the impact of different contaminants in the wild (Medaković et al., 2006;Valbonesi et al., 2003) and further used as a biofilter in aquaculture for removal of dissolved and particulate nutrients from fishpond effluents (Cranford et al., 2011;Shpigel, 2005). Moreover, it is a highly appreciated gastronomic and nutritional resource, representing one of the most important aquaculture species of the 20th century.
However, in the past 40 years, the production of European flat oyster has declined drastically worldwide from the peak production of nearly 30,000 t in 1961 to just 2000 t in 2012 (FAO, 2020), due to a combination of overexploitation, habitat destruction, environmental change and diseases. European flat oyster aquaculture has been increasingly replaced by other oyster species such as Crassostrea gigas (Troost et al., 2008) and flat oyster now accounts for <0.2% of the total global production of all farmed oyster species (5 Mt in 2015; FAO, 2020).

Various diseases impact both wild and cultured populations of
O. edulis (Sas et al., 2020). The two main disease-causing organisms-Martelia refringens (marteliosis) and Bonamia ostreae (bonamiosis)started to spread in the 1970s and 1980s, respectively, and despite changes in management practices and intensive restocking programs, production has remained steady and low in recent years. Bonamia ostreae is a destructive intracellular parasite transmitted horizontally and putatively also vertically, causing a heavy inflammatory reaction that leads to high mortality following the loss of normal architecture and function in multiple organs (da Silva et al., 2008). The parasite proliferates inside the haemocytes occurring mostly in the connective tissue but also in epithelia of numerous organs. Due to its small size (3-5 μm in diameter) and its intracellular location, B. ostreae infection is often difficult to detect visually, which is why standard diagnostic methods use cytology and histopathology together with qPCR to screen oyster tissues (Sas et al., 2020). Although adult individuals appear to be the most susceptible to mortality in the field, the parasite can also infect larvae, seeds and juvenile oysters, suggesting that larvae displacement with ocean currents may result in the movement of the parasite through Europe (Arzul et al., 2011).

However, the patchy distribution of the parasite in the Northeast
Atlantic Ocean suggests that infection in one area does not necessarily ensure the transference of the parasite to neighbouring areas.
Understanding the genetic basis of host response to parasites and host-parasite interactions during infection is essential to ascertain the best strategy to control diseases. Several studies have supported that both in vitro contact of B. ostreae with haemocytes and its inoculation into healthy oysters activates the immune response of the host, inducing up-regulation of immune-related genes and modification of cellular parameters (Pardo et al., 2016;Ronza et al., 2018).
Following an oligo-microarray approach, genes involved in resistance such as antimicrobial peptides and extracellular matrix components were identified by comparing gene expression profiles of naïve and long-term affected oyster strains (NV and LTA, respectively;Ronza et al., 2018). Breeding programs to produce oysters resistant to bonamiosis rely on the underlying genetic variation and have been successfully applied by selecting broodstock among survivors after long exposure in heavily affected areas (Ronza et al., 2018). Quantitative trait loci (QTL) and markers associated with resistance to bonamiosis have been reported, and this information could accelerate the selection process (Lallias et al., 2009). Combining functional and markerassociation genomic strategies, Harrang et al. (2015) also identified SNP (single nucleotide polymorphism)-associated markers linked to differential expression of candidate genes (eQTL).
Recently, a SNP array containing a total of 14,065 putative flat oyster SNPs (Gutiérrez et al., 2017) and 37 SNP-associated markers from candidate genes identified by Ronza et al. (2018) were analysed to identify signals of divergent selection to bonamiosis by comparing two sets of NV and LTA flat oyster populations (Vera et al., 2019).
Within a low differentiation genomic background in the Atlantic area (F ST = 0.0079; Vera et al., 2016), a set of 22 consistent and 87 suggestive outliers were associated with divergent selection to bonamiosis.
These markers discriminated NV vs. LTA populations, providing useful information for their putative application in breeding programs K E Y W O R D S bonamiosis, chromosome rearrangement, disease resilience, Ostrea edulis, SNP-chip to obtain parasite resistant strains. Interestingly, the strong linkage disequilibrium observed between most of these markers suggested that a major QTL could underlie resilience to the parasite.
In this study, we used the flat oyster SNP array to validate the genetic markers associated with resistance to B. ostreae reported by Vera et al. (2019) in a narrower sampling scenario, where populations of have been subjected to different bonamiosis pressures.
We used a new chromosome-level flat oyster genome assembly to map validated outliers and differentially expressed genes (DEG) previously reported by Ronza et al. (2018) for an integrative analysis.
The results confirmed the association of outliers with the bonamiosis status of samples in an independent scenario and demonstrated their location in a single genomic region containing genes enriched in relevant immune functions over the genomic background.

| Sampling
The analysis was performed on 134 flat oyster individuals collected in 2019 from locations with different bonamiosis prevalence: (i) three oyster beds from the North Sea (OS, GBR, NO); (ii) two hatchery batches (OSH, WZH) from the same area (FAO subarea 27.4 https://www.fao.org/fishe ry/en/area/27/en; Figure 1) ( Table 1). Hatchery spat was produced from wild individuals of the same area that were conditioned in the hatchery until spawning.
Hatchery samples were included in the study considering their interest for restocking bonamiosis-affected areas using Bonamia-resilient strains from marker-assisted selection programs. These samples provided additional information for the association analysis and validation of markers, and specifically in one case, the Wadden Sea, the only representative sample was from the hatchery. Furthermore, hatchery data allowed an estimation of the contribution of parents for checking the reproduction protocols in batches at hatcheries to maintain genetic diversity in future restoration programs. The sampling locations were sorted according to bonamiosis status from previous information as: (i) naïve (NV), supposedly to have never been in contact with the parasite and (ii) long-term affected (LTA), where parasite presence was first reported over 30 years ago and monitored over time to the present (Engelsma et al., 2010).

| DNA extraction and SNP genotyping
Genomic DNA was extracted from gill tissue using the E.Z.N.A.® Mollusc DNA kit (OMEGA Bio-Tek) following the manufacturer's recommendations. Briefly, between 30 and 50 mg of fresh gill was treated with a solution containing 20 μl of proteinase K and 700 μl of lysis buffer (CSPL) overnight while gently shaking to mix thoroughly at 65°C. Samples were subjected to RNAse A treatment and to several buffer washes by centrifugation. DNA quantity and purity was measured with a Nanodrop spectrophotometer and DNA concentrations were normalized for genotyping.
Samples were shipped to Identigen (Dublin, Ireland) for genotyping on a medium-density Affymetrix Axiom SNP array (14,950 SNPs; 1 SNP per ~63 kb of the flat oyster genome; 935.6 Mb; see below). After filtering SNPs with <90% call rate, 11,641 were retained, which supports the robustness of the SNP chip platform (Gutiérrez et al., 2017). A total of 611 SNPs were monomorphic across all sampling locations and discarded for structure analyses, as they would not influence the results. Three different SNP data sets were considered according to the analyses performed

| Genetic diversity
Genetic diversity per sample was estimated using expected (H e ) and observed (H o ) heterozygosity, and allelic richness (Ar), computed using the rarefaction method. These analyses were performed with the DiveRsity R package v. 1.9 using the 'basicStats' function (Keenan et al., 2013). Non-parametric Mann-Whitney and Kruskal-Wallis tests were performed to check for genetic diversity differences between pairs or groups of samples, respectively.
Conformance to Hardy-Weinberg expectations was evaluated with an exact test implemented in the R package Genepop v 1.1.7 (Rousset, 2008). This test was performed using the complete enumeration method to estimate p-values (Louis & Dempster, 1987).
Global p-values per sample were obtained with Fisher's method (Fisher, 1933). Deviations from panmixia were estimated for each locus and averaged across loci for each sampling location using the intrapopulation fixation index (F IS ; Wright, 1949). The significance of F IS was tested with the 95% bias corrected confidence intervals (CI) obtained after 1000 bootstraps iterations with 'divBasic' function (DiveRsity R package). A full-pedigree reconstruction method implemented in Colony v2.0.6.6 (Wang, 2004) was applied to estimate parentage, and the number of breeders in all locations/hatchery batches. We used 1000 highly polymorphic SNPs (Minimum Allele Frequency [MAF] > 0.35), homogeneously distributed across the genome, among the 11,641 SNPs available due to the computational limits of the program working with thousands of markers. The maximum likelihood method implemented in Colony was also used to estimate effective population size assuming random mating.

| Genetic structure
Pairwise interpopulation fixation index (F ST ) between sampling locations was calculated using the StaMPP R package v 1.6.2 with the 'stamppFst' function (Pembleton et al., 2013); 10,000 bootstrap replicates across loci were used to generate 95% confidence intervals and p-values regarding the null hypothesis (F ST = 0). Global F ST was calculated using the R package Genepop with the 'Fst' function (Raymond & Rousset, 1995).
The number of different genetic clusters in our sampling collection was evaluated with STRUCTURE v 2.3.4 (Pritchard et al., 2000).
This program looks for the number of different genetic clusters (K value) in the samples compatible with the data by analysing linkage disequilibrium and Hardy-Weinberg equilibrium following a Bayesian clustering approach. Tests without a priori information were performed using the R package ParallelStructure v 1.0 (Besnier & Glover, 2013) with a burn-in of 100,000 iterations and 200,000 Markov Chain Monte-Carlo steps. Ten independent replicate runs were used to increase precision. The ADMIXTURE ancestry model was performed with K1 − K6 (N samples +1). To identify the most likely number of clusters, two different K estimators were used: del-taK (Evanno et al., 2005) and Mean LnP(K) (Pritchard et al., 2000).
StructureSelector Web-based software (Li & Liu, 2018) was used to obtain K estimators and CLUMPAK graphical outputs (Kopelman et al., 2015). Moreover, an analysis with a priori population information was performed (POPINFO = 1) on our samples with K = 2 using consistent and suggestive outliers taking as references the LTA Discriminant analysis of principal components (DAPC), a multivariant method to infer the number of clusters in a sample of genetically related individuals, was applied as a complementary method to infer structure in the samples studied using the Adegenet package function 'dapc' in RStudio (Jombart & Ahmed, 2011). A principal component analysis (PCA) from the matrix of genotypes was performed and then, a selected number of principal components (PCs) used as input for linear discriminant analysis (LDA). The selection of the optimal number of PCs to be further used in the LDA was done via cross-validation, and those associated with the lowest root mean square error (RMSE) were retained. In addition, DAPCs retaining at least 90% of the cumulative variation of the data were evaluated.

| Linkage disequilibrium
Pairwise linkage disequilibrium (LD) between outliers mapping in the same chromosome was estimated using the square of the correlation TA B L E 1 Characteristics of Ostrea edulis samples analysed in this study coefficient r 2 as implemented in the R package gaston v 1.5.7 with the 'LD' function (Perdry & Dandine-Roulland, 2020). The results were represented in a LD heatmap plot using the same R package.
A p-value for estimated deviations from the null hypothesis (r 2 = 0) was obtained using an exact test for genotyping disequilibrium with Genepop 4.7 (Rousset, 2008). LD between consistent outliers was also estimated using the previous data set from Vera et al. (2019), considering that they represented the original reference and also the broader sampling area studied. 93.6% of total genome assembly) corresponding to the haploid chromosome number of flat oyster (2n = 20; Leitao et al., 2002). All SNPs, including the consistent and suggestive outliers detected for resistance to bonamiosis, and the DEGs between NV and LTA oysters were mapped against the OEROSLIN genome to check for co-localization.

|
Functional enrichment of genes located at genomic regions associated with bonamiosis resistance was obtained with R package GOfuncR version 1.14.0 (Grote, 2021). Enriched Gene Ontology (GO) terms related to 'Biological Process', 'Molecular Function' and 'Cellular Component' and the associated genes were ascertained against the background of all genes in the flat oyster genome that received GO annotations.

| RE SULTS
The present study aimed to validate previous Bonamia-resistance associated markers reported by Vera et al. (2019) in an independent narrower scenario following a population genomics approach and exploiting a new flat oyster genome assembly for a more comprehensive interpretation of results. Different SNP panels were used according to the goal of the analyses performed ( Table 2). The whole SNP data set, including polymorphic and monomorphic markers, was used to estimate genetic diversity, which made possible comparisons with previous studies in flat oyster, while the neutral and outlier data sets were used to investigate genetic structure with respect to the differential bonamiosis pressure operating on flat oyster samples (i.e. NV vs. LTA). The number of polymorphic markers available for the different classes (i.e. neutral and outlier) was lower than reported by Vera et al. (2019) because of the smaller geographical range analysed.
In particular, the number of suggestive and consistent SNP outliers dropped from 87 to 71 and from 21 to 16, respectively ( Table 2).

F I G U R E 2 STRUCTURE analysis in
Ostrea edulis from the North Sea with: (a) neutral for K = 3 and (b) 16 outlier, and (c) 71 outlier SNP data sets for K = 2. Each vertical bar represents a single individual and its colour proportion the posterior probability of its assignment to the defined STRUCTURE identified The genetic structure of flat oyster samples was also tested with DAPC analyses (Figure 4). For the neutral panel, 101 PCs that explained 90.4% of the variance were used, and the higher differentiation component was found between sampling locations, with the wild NO being the most divergent (Figure 4a). For the consistent outlier panel (16 SNPs), we used 6 PCs that explained 90.5% of the variance (Figure 4b), and in this case, a great similarity was found between sampling locations, despite a much higher dispersion of individuals regarding the population centroids than observed with the neutral data set. In fact, the more distant individuals within each location mostly corresponded to the 'orange' samples detected in the STRUCTURE analysis, a pattern matching the barplot with K = 2. A very similar representation was observed with the suggestive outlier data set (data not shown).
From the STRUCTURE information, we split the whole sample into two subgroups pertaining to the two clusters observed with both outlier data sets: 27 'orange' individuals associated to the LTA cluster and 107 'blue' individuals to the NV one. For this, we applied the individual membership coefficient to each cluster (q value) using a threshold of q "orange" > 0.5 and q "blue" ≥ 0.5 to classify the few admixed individuals observed (Figure 2b,c). Genetic diversity (H e ) was almost twofold higher in the 'orange' than in the 'blue' cluster with both outlier data sets, while the F IS was negative and highly significant in the 'orange' cluster, indicating a heterozygote excess, but not in the 'blue' one, which conformed to Hardy-Weinberg expectations ( Table 5). The deviation was nearly double for the consistent outlier data set (−0.376) than for the suggestive one (−0.196), while the 'blue' group did not depart from F IS = 0 in both cases. Furthermore, we also performed a DAPC analysis with the orange and blue groups, K = 1 being obtained when using the neutral panel, while K = 2 being the most likely explanation when using both outlier data sets.

| Genomic mapping of DEGs and outliers: linkage disequilibrium
Using the flat oyster genome, we mapped the outliers associated with bonamiosis resistance along with the differentially expressed genes (DEGs) detected by Ronza et al. (2018) when comparing the response of NV and LTA strains to a bonamiosis challenge for an integrative analysis (Table S1A). An important proportion of the 71 suggestive and nearly all the 16 consistent outliers associated with the bonamiosis status of sampling locations mapped at OE-C8 (37 and 14, respectively) (  Figure S3).
Interestingly, most SNPs located in that region showed highly significant F ST values between 'orange' vs. "blue" groups and highly significant negative F IS values within groups as compared to the whole genome background using average F ST and F IS values across sliding windows of 50 SNPs ( Figure 6). This enabled us to refine the differential region between both groups up to 33.83 Mb (Table S1A) and suggested that several other SNPs in that region should likely be outliers between LTA and NV groups ( Figure 6).
Moreover, among the 711 differentially DEG reported by Ronza et al. (2018) mapping on the flat oyster genome, 65 were located on OE-C8, representing the highest gene density (DEGs/ TA B L E 5 Genetic diversity for the "orange" and "blue" groups identified in Ostrea edulis from the North Sea with STRUCTURE for the 16 (A) and 71 (B) outlier panels

Note:
In parentheses two outliers with the same matching score in two different chromosomes counted twice. total genes; Table 6), although nonsignificant when tested across  (Table S1A). Three additional methyltransferase genes were detected within the 29.44 Mb region along with genes with functions linked to endo-and exocytosis and cell adhesion.
We also investigated the overrepresented GO terms on genes from OE-C8 in comparison to the flat oyster genome background.
Enriched biological processes included those with immune functions, such as interleukin-4-mediated signalling pathway, negative regulation of natural killer cell activation, regulation of peptidoglycan recognition protein signalling pathway, negative regulation F I G U R E 5 Linkage disequilibrium (r 2 ) between all pairs of suggestive outliers at OE-C8 of Ostrea edulis. LD intensity (r 2 from 0 to 1) is shown using a range from whitish to reddish colours, respectively, as shown in the scale of the figure F I G U R E 6 F ST values between 'orange' vs. 'blue"'groups and F IS values within groups with all the 9277 SNPs mapping in the 10 Ostrea edulis chromosomes using average values of sliding windows of 50 SNPs across all chromosomes (C) of response to interferon-gamma and chromosome segregation (Table S2). Terms associated with interleukin-4 signalling and response to interferon-gamma were explained by genes encoding histone-lysine N-methyltransferase setd1b and multiple copies of Poly (ADP-ribose) polymerase member 14. Several genes encoding peptidoglycan recognition proteins 1 and 2, along with lectin galactoside-binding soluble protein 4 were responsible for the enrichment of terms associated with natural killer cell activation and peptidoglycan recognition signalling. Finally, the enriched term 'chromosome segregation' was explained by a large set of genes, including those encoding anaphase-promoting complex 23, structural maintenance of chromosome protein 4, ATP-dependent DNA helicase MER3, cyclin-B3, meiosis sister chromatid cohesion complex protein, mitotic spindle checkpoint protein Bub3, MMS19 excision repair mechanism protein, dynactin subunit 1 and retinoblastomaassociated protein.

| DISCUSS ION
European flat oyster production is severely threatened in the Northeast Atlantic by the intracellular protozoan B. ostreae, (Engelsma et al., 2014). Within a very low, but significant genetic differentiation scenario in the Atlantic area (F ST = 0.0079; Vera et al., 2016), flat oyster populations that were never in contact with the parasite (naïve: NV) coexist with populations affected for many generations (long-term affected: LTA), dating back to the first report of the parasite in Europe Pichot et al., 1980). This observation suggests that while larvae can be infected by the parasite (Arzul et al., 2011), its transmission to other areas is not en- group showed twice the genetic diversity of the 'blue', coupled with a highly significant heterozygote excess, while the 'blue' group did not depart from Hardy-Weinberg expectations. This observation might suggest an overdominance mechanism underlying resistance to bonamiosis, but further work should be done for its confirmation. Tolerance or resistance to stress or pathologies has been associated with overdominance in other studies (Di et al., 2015;Gallaga-Maldonado et al., 2020;Maynard et al., 2016), sometimes involving specific immune genes, such as the histocompatibility complex (MHC; Kekäläinen et al., 2009;Penn et al., 2002). It is important to note that the term resistance applied until now in this and previous studies related to bonamiosis infection is rather inaccurate, as the response to the parasite may be related both to resistance against infection and the ability to survive once infected (tolerance) (Holbrook et al., 2021). A more appropriate term encompassing both components is resilience, and we suggest its use until more refined  (Fontdevila et al., 1983;Maroso et al., 2018;Mérot et al., 2021).
Flat oyster responses to B. ostreae have been associated with preventing the entrance of the parasite into haemocytes, and with apoptotic processes arresting proliferation of the parasite once it surpasses the external barrier (Gervais et al., 2016;Morga et al., 2012;Ronza et al., 2018). In our study, a set of 65 DEGs between naïve and long-term affected strains (Ronza et al., 2018) were located on OE-C8, representing the highest proportion among the 10 flat oyster chromosomes. To note, a block constituted by five genes close to the 33.83 Mb refined region at OE-C8 included several collagen-related genes involved in the extracellular matrix (external barrier of haemocytes) and a programmed cell death protein (avoidance of parasite proliferation). Interestingly, a SNP associated to the collagen-IV DEG (Ronza et al., 2018) was the only outlier associated with divergent selection among the 37 candidate genes evaluated by Vera et al. (2019). Also, several methyltransferases involved in epigenetic patterns were detected at OE-C8 and could explain the faster response of LTA oysters to bonamiosis reported by Ronza et al. (2018), considering the putative transmission of epigenetic patterns across generations (Eirín-López & Putnam, 2019). Furthermore, the functional enrichment analysis performed on OE-C8 also identified GO terms such as interleukin-4 regulation signalling pathway involving five genes related to epigenetic marking and a set of 26 genes related to chromosome segregation, many of them associated with apoptotic defence mechanisms.

| CON CLUS ION
Our data adds new evidence to the divergent selection outliers detected previously by Vera et al. (2019) associated with the bonamiosis status of samples on an independent sample collection. Also, the suggested genomic region associated with bonamiosis resistance by the same authors was supported and most outliers were located on a region of OE-C8. These outliers displayed a strong LD over a large chromosomal region, which supports a chromosome reorganization at OE-C8 that could be maintained by overdominance related to B. ostreae resilience. This genomic region include genes encoding proteins related to apoptosis and extracellular matrix components, considered key in the response to an intracellular infection. We also detected significant enrichment GO terms involving genes controlling epigenetic marks that could be related to the quick response to bonamiosis observed in LTA oysters. The identified markers at OE-C8 associated with B. ostreae resilience could aid for recovering its production and ecosystem services throughout the European coast.
However, further work is needed to test this association at individual level, to disentangle the components underlying bonamiosis resilience and to confirm the chromosome reorganization suggested at OE-C8 supporting an overdominance mechanism. Finally, the putative involvement of epigenetic mechanisms on the transgenerational immune memory for a quick response to the parasite is a suggestive hypothesis that deserves further attention.

CO N FLI C T O F I NTE R E S T
The authors declare that they have no competing interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
All genotyping information of this manuscript including neutral, outlier and bonamiosis-associated markers is available at the Dryad