The 48 families of M. edulis used in this study are detailly described in (Ajithkumar et al., 2024b). Briefly, three wild mussel populations were sampled, two from the Atlantic coast (OLE-PON and YEU_001), and a third from the North of France (WIM), and transferred to the Ifremer hatchery in La Tremblade in the fall of 2016. Each mussel population was cleaned and placed in separate tanks containing unheated UV-treated, and filtered seawater (400 L per hour). To favor gametogenesis, mussels were fed a cultured phytoplankton diet (Isochrysis galbana, Tetraselmis suecica, and Skeletonema costatum). Two sets of crosses were performed in January 2017 (set 1) and in February 2017 (set 2). For each population, 100 mussels were individually placed in 400 mL beakers, and spawning was triggered by alternating cold (10°C) and warm seawater (20°C). Depending on the ripeness of the mussels and the sex ratio, four males for OLE-PON, and 11 males for YEU_001 were used in set 1, while 9 males were used for WIM in set 2. Within population, each male was mated with two females, producing in total 24 half-sib families, each containing two full-sib families. These 48 full-sib families consist of 8, 22 and 18 families from OLE-PON, YEU_001 and WIM population, respectively, all should be considered as outbred. Each full-sib family was grown separately in 30 L tanks filled with filtered and UV-treated seawater at 20°C until the pediveliger stage. Then, downwelling system were used until mussels reached 1 cm. At that size, they were transferred to our nursery in Bouin in April and May 2017 for set 1 and set 2, respectively. For each family, 1,000 spat were maintained in 15 L SEAPA^© baskets, and all families were raised in a 20 m³ concrete raceway until the start of the experiment, which occurred in July 2018. More detailed on the larval and nursery culture are provided in Ajithkumar et al. (2024b).

Detailed step-by-step protocol of the experimental infection is given in Ajithkumar et al. (2024b). Briefly, two experimental infections (EI_1 and EI_2) were conducted in July 2018. Among the 48 families (mean individual total weight of approximately 5 g), 24 families were randomly sampled and tested in El_1, and the others were tested in El_2. Additionally, a third experimental infection (EI_3) was performed, involving the same 24 families from El_1, to increase the genotype sample size due to some samples from El_1 were lost. Finally, three families tested in El_2 were also tested in El_1, leading to a total of 27 families tested in El_1, and one of those three was also tested in El_3 indicating that only one family was tested in each of the three experimental infections (see Supplementary Table S1 for details). For each experimental infection and family, a highly pathogenic strain of V. splendidus (strain 14/053 2T1) isolated during AMM outbreak in 2014 was injected in 30 mussels to investigate their resistance to this pathogen. This particular strain was selected for its high virulence, the highest among the 50 strains tested, and it belongs to a specific pathogenic group. The lethal dose-50 (LD₅₀) was determined prior to the experimental infection by injecting varying concentrations of Vibrio splendidus strain into a large group of mussels. The selected dose was close to the LD₅₀, ensuring adequate variance in the response for effective assessment of resistance. First, mussels were anesthetized using MgCl₂ (50 g per L), and 50 μL of bacterial solution (10⁹ bacteria/mL) was injected into the muscle. Then, ten injected mussels per family, for all the 24 families of one set were hold in one 120 L tank containing UV-filtered seawater. Three replicate tanks were used and, in each tank, water recirculation was maintained using a TECO®pump (Ravenna, Italy), which also maintained the seawater temperature at 17°C. Dead mussels were counted and sampled daily up to 72 h post-injection. The adductor muscle/gills of the dead mussels during the experiment and the surviving mussels at the end of the experiment were collected using scalpels disinfected with 70% ethanol and stored in 1.5 mL sterile tubes at room temperature.

Among the dead and alive animal sampled, around 14–15 mussels were genotyped per family (ranging from 13 to 16) (Table 1; Supplementary Table S1). The number of dead and alive per family was calculated in proportion to the mortality observed within each family in El_1 and in El_2. As the number of animals sampled in El_1 was not sufficient, 158 dead mussels and 35 live mussels were sampled from EI_3 (Supplementary Table S1). In total, 768 individuals were genotyped; 348 dead, 348 alive, and the remaining 72 were their parents (48 dams and 24 sires). DNA extraction and genotyping were performed by the Gentyane INRAE Platform (Clermont-Ferrand, France) using the multi species medium-density 60K SNP-array, Axiom_Myt_v1_r1 (Thermo Fisher Scientific, Waltham, Massachusetts, United States), which comprises 23,252 markers for M. edulis (Nascimento-Schulze et al., 2023). Quality controls on the 60K SNPs from the SNP array and genotyped individuals were performed as described in D'Ambrosio et al. (2019). Firstly, genotypes of all individuals were analyzed using the Axiom Analysis Suite software (AxAS; v.4.0.3.3) with the default best practice workflow suggested by the manufacturer, with few threshold modifications, which includes individual quality control and SNP quality control analysis (Dish Quality Control ≥0.20; Quality Control call rate ≥90; percent of passing samples ≥98; average call rate for passing samples ≥92%; call rate cutoff ≥95; Fisher’s Linear Discriminant ≥2.6). Consequently, 7,476 polymorphic SNPs were retained for further analysis. Subsequently, final quality control was performed using PLINK v1.9 software (Chang et al., 2015). Two individuals with an identity-by-descent value over 0.90 were considered as duplicated and both individuals were removed from the analysis. Only SNPs with a minor allele frequency (MAF) higher than 0.01 and those passing the Hardy-Weinberg equilibrium test (p-value <0.0000001) in the genotyped mussels were retained. After the quality control, data comprised of a total of 766 genotyped individuals for 3,406 SNPs.

Parentage assignment was performed in the R package APIS (Griot et al., 2020) with a mismatch number set to 5%. The best 1,471 SNPs, selected with call rate greater than 90% and MAF value greater than 0.1 were used. Parentage assignment allowed the reconstruction of the pedigree of 647 offspring with assignment rates reaching 93.2% of the mussels having both parents assigned, while the remaining 47 mussels potentially from outside the cross-mating design, which were excluded from following analyses.

To evaluate potential genetic sub-structuring of populations and any associated biases, a principal component analysis (PCA) was performed using PLINK 1.9 (Chang et al., 2015) and the genetic structure was visualized using in RStudio (Team, 2024). Three individuals were identified as outliers beyond the population structure and were subsequently excluded from further analysis. Genetic differentiation between populations was measured through pairwise fixation index (F_ST) estimates using PLINK 1.9 (Chang et al., 2015).

All markers of the array along with their flanking regions were blasted using a BLASTn^® procedure on the reference genome (Mytilus edulis genome assembly, xbMytEdul2, GenBank accession number: GCA_963676595.2). To map SNPs, considering the high polymorphism in the mussel genome, four mismatches were allowed over a length of around 71 base pairs. Only SNPs mapping to a unique position on the reference genome were retained for the subsequent stage of quality control as mentioned in previous section. Out of the 3,406 SNPs, only 2,204 matched our mapping criteria and were successfully positioned on the reference genome (Supplementary Table S2).

The pairwise linkage disequilibrium (LD) analysis was performed between all SNPs and adjacent markers for each linkage group and population to determine LD decay within the genome of M. edulis using Plink 1.9 (Chang et al., 2015).

To identify SNPs associated with resistance to V. splendidus, a genome wide association study (GWAS) was performed using a mixed linear model association through ssGBLUP analysis. The postGSF90 module (Misztal et al., 2014) from the BLUPF90 package was used to estimate the effects of the SNPs ( a ^ i ) based on the genomic breeding values g ^ i predicted for the genotyped animals. The SNP effects were estimated according to the following equation: a ^ i = d Z ′ Z d Z ′ − 1 g ^ i where d is the vector of weights associated with the SNP effects and Z is the incidence matrix relating SNP effects to genomic breeding values.

Estimates of SNP effects ( a ^ i ) can be used to estimate the proportion of additive genetic variance of each SNP effect: % V a = 2 p 1 − p   a ^ i 2 σ a 2 * 100 with σ a 2 the total genetic variance estimated using the linear mixed model with postGSF90 and p the minor allele frequency of the target SNP.

Prediction accuracy for the PBLUP, GBLUP, and ssGBLUP models was assessed using the ‘leave-one-out’ method, implemented with a linear model in BLUPF90 (Resende et al., 2012; Kristensen et al., 2018). In this approach, each observation is systematically excluded one at a time. The model is then trained on the remaining data, and the (G)EBV for the excluded individual is predicted by masking its phenotype.

The accuracy (r) of prediction was computed as the correlation between the (G) EBVs and the corrected phenotype ( y ^ ) of the mussel divided by the square root of the heritability, using the formula: r = G EBV , y ^ h 2

The heritability value h 2 used in this analysis was calculated using the variance components ( σ a 2 and σ e 2 from the ssGBLUP model.

The cumulative mortality rate 72 h post-injection was 47%. At endpoint, mortality rates were 63% for EI_1, 41% for EI_2, and 37% for EI_3. Among mussel populations, the WIM population (54%) showed higher susceptibility to V. splendidus compared to the YEU_001 (45%) and OLE-PON (37%) populations (Supplementary Table S4). Mortality rates varied significantly among families upon exposure to V. splendidus, ranging from 17% to 83%. The mean mortality rates for all families are depicted in Figure 1.

Figure 2 illustrates the results of the principal component analysis (PCA), revealing the population structure of the mussel population. The first two PCA axes collectively account for over 15% of the total genetic variation. The YEU_001 and OLE-PON populations were generally homogeneous, although there were few stratifications within YEU_001. The WIM population have quite different genetic background from others, with two families whose offspring showed even greater isolation. Nevertheless, F_ST analysis revealed low genetic differentiation between populations. The mean genetic distances between populations are shown in Table 2, with F_ST values ranging from 0.02 to 0.03, suggesting genetic similarity across all three populations (Figure 3).

2,204 SNPs were positioned on the reference genome, resulting a loss of 1,202 SNPs. The positions of markers on the chromosomes are illustrated in Supplementary Figure S1. The average SNP density per megabase (Mb) ranges from 0.57 to 2.37, varying among chromosomes and within chromosome (Supplementary Table S5). Approximately, only 9% of all 1 Mb segments contain more than 5 SNPs. SNP density exhibits non-uniformity throughout the genome, with each chromosome demonstrating varying densities. The lower marker density results in greater mean average distances between adjacent SNPs, ranged from 421 kb to 1739 kb depending on the chromosome.

Figure 4 illustrates that linkage disequilibrium (LD) decreases sharply as the distance between pairs of SNPs increases, with the most rapid decline occurring within the first 100 kb. Beyond this range, LD continues to decline and becomes more variable. The OLE-PON population consistently shows higher LD throughout the genome compared to other populations. On average, the LD values (r ²) for SNPs less than 15 kb apart are 0.12 for OLE-PON, 0.10 for WIM, and 0.06 for YEU. Linkage disequilibrium values are generally low between adjacent SNPs for all the chromosomes, where distances between adjacent SNPs are larger.

The estimates of heritability using the linear and Gibbs sampling models are summarized in Table 3. Pedigree-based heritability estimates for resistance to V. splendidus in M. edulis ranged from 0.22 to 0.31. Genomic heritability was slightly higher, varying between 0.33 and 0.36. The ssGBLUP based estimated heritability ranging from 0.28 to 0.33, which combines genomic and pedigree information.

GWAS for resistance to V. splendidus suggest that this trait is likely to be impacted by multiple genomic regions. However, there were 20 SNPs that explained >0.5% of additive genetic variance on chr 1, chr 2, chr 3, chr 4, chr 5, chr 6, chr 8, chr 9, chr 13, and chr 14 (Figure 5). However, none of these markers explained more than 1.06% of the additive genetic variance (Figure 5; Table 4).

Accuracy with all data are 0.36, 0.43, 0.43 for PBLUP, GBLUP and ssGBLUP, respectively. Genomic selection (GBLUP and ssGBLUP) is better than PBLUP by 19%. Overall, prediction accuracy for genomic selection increased with the density of markers (Figure 6). Incorporating genomic information generally enhanced accuracy compared to pedigree-based estimation, except with 500 SNPs where PBLUP exhibited higher accuracy than GBLUP (Figure 6). With maximum training population and SNP subsets, genomic evaluation improved accuracy by 17%, 19%, 25%, and 19% for 1,000, 1,500, annotated (2,204), and all SNPs (3,406), respectively, compared to PBLUP. When comparing GBLUP and ssGBLUP models, the prediction accuracy was consistently favored the ssGBLUP model, except when using annotated SNPs in the GBLUP model (Figure 6).

To the best of our knowledge, our study is the first to use the multi species medium-density 60K SNP-array (Nascimento-Schulze et al., 2023) to estimate genetic and genomic parameters in blue mussel (M. edulis). After all quality controls, we identified high-quality SNPs from 23,252 SNPs present on the array across 768 individuals. The necessity for stringent filtering of genotyping data is highlighted by the prevalence of poor-quality markers. One of the distinguishing features of the Mytilus species complex is its high degree of heterozygosity and abundance of mobile elements, which pose challenges in accurately identifying polymorphic SNPs (Gerdol et al., 2020; Sun et al., 2017). Compounding this difficulty, the current SNP array was developed using a whole-genome low coverage approach, which can lead to loss of polymorphic sites, especially across genetically diverse populations. A prior study validating this SNP array found 89.8% of the 23,252 SNPs to be polymorphic due to the close genetic match between the populations used for development and validation (Nascimento-Schulze et al., 2023). In contrast, after quality control in our study, only 14.6% of the SNPs were retained as polymorphic, likely due to genetic differences between our sampled population and the one used for array creation. This discrepancy highlights the influence of population structure on SNP retention and substantial genetic complexity and high polymorphism exhibited in mussel genome (Gerdol et al., 2020; Smietanka et al., 2014). This complexity stems from its evolutionary history of extensive gene flow among related species and a genome characterized by significant dispensable content. Indeed, about 30% of genes in Mytilus are dispensable, likely impacting SNP distribution and usability across populations (Gerdol et al., 2020). Out of 23,252 SNPs identified in M. edulis, without any other quality control, only 16,213 (70%) were annotated on the recently published reference genome of M. edulis. Assembly errors in the reference genome may rise from several factors, such as exceptionally high genetic polymorphism levels, non-Mendelian segregation of marker loci in paired crosses, and a significant occurrence of null alleles in genetic markers (Hedgecock et al., 2015). If we consider the other quality controls done, only 2,204 SNPs (64.8% of polymorphic, good quality SNPs, 9.4% of the total available SNPs) mapping successfully. This limited SNP coverage and sparse distribution across the linkage map could potentially lead to the omission of QTLs in specific regions, particularly in regions where markers are underrepresented. To improve SNP array quality and increase the number of informative SNPs, future work could include strategies such as pooled sequencing with a greater number of individuals, incorporating genetically diverse populations, or performing high-coverage whole-genome resequencing numbers of individuals from each population (Nascimento-Schulze et al., 2023). These adjustments would increase the SNP array’s robustness, thereby improving genome coverage, QTL detection power, and marker-trait association accuracy in M. edulis.

The bi-species Axiom Affymetrix 57K SNP array has been used in Pacific oysters, where applying the AxAS software’s best practice workflow led to a notable reduction in the number of informative SNPs. Specifically, Gutierrez et al. (2018) reported 23,000 informative SNPs from 820 individuals, Vendrami et al. (2019) identified 21,499 SNPs from 232 individuals, and Jourdan et al. (2023) obtained 14,500 SNPs from 2,420 individuals. This reduction is largely attributed to the complex genetic structure of molluscs, stemming from the highly polymorphic nature of their genomes (Gerdol et al., 2020; Jiao et al., 2021; Song et al., 2021), and is further influenced by the genetic relationship between the training population used for array design and the breeding candidates in selective breeding program (Houston et al., 2020). However, recent studies on bivalves have demonstrated that a moderate number of high-quality markers (1,000–3,000) could suffice for accurate predictions (Gutierrez et al., 2018; Kriaridou et al., 2020; Penaloza et al., 2022).

Linkage disequilibrium (LD) at the genome level plays a crucial role in the efficacy of breeding programs, influencing genetic variance and the accuracy of association analyses (Goddard and Hayes, 2009; Siol et al., 2017). In our study, values for r ² ranged between 0.07 and 0.09 for SNPs within a distance of 10 kb and from 0.03 to 0.08 within 50 kb in the populations studied, suggesting limited short-range LD at current SNP marker densities. LD levels decreased to less than 0.05 at 100 kb in two populations. The high LD in OLE-PON is probably an overestimate due to a smaller effective population size (N_e) than in the other populations. Overall, LD between adjacent markers within each population was predominantly less than 0.1 within 2 kb, indicating a rapid decline in LD within the blue mussel genome. Therefore, these LD estimates should be treated with caution. This swift decay suggests a historically large effective population size and high recombination rate, reflecting substantial genetic diversity within the population (Ellegren and Galtier, 2016). Moreover, LD values are population-specific, and influenced by evolutionary factors such as natural selection, mutation, genetic drift, origin and migration, as well as molecular forces including historical recombination events, and breeding history such as historical effective population sizes, intensity and direction of artificial selection, population admixture, and mating patterns (Du et al., 2007). Our findings confirm the low LD in M. edulis populations, consistent with previous studies on bivalves (Jones et al., 2013b; Vera et al., 2022; Jourdan et al., 2023). However, the sparse SNP distribution and generally low LD values (r ² < 0.1) suggest that increasing marker density is essential to improve prediction accuracy and QTL detection.

F_ST is widely applied to evaluate genetic differentiation between/among populations (Hu et al., 2022). The low F_ST values (F_ST < 0.03) observed in our study suggest minimal genetic differentiation among mussel populations, indicating a lack of significant genetic structure. This phenomenon may be attributed to similar selection pressure and limited gene flow among the mussel populations, irrespective of geographic location. Similar to our findings, wild populations of M. edulis sampled from the Atlantic coast of France and from the Northern France to Germany had low Fst value (<0.003) using ancestry-informative SNPs (Simon et al., 2020). In addition, similar findings have been reported in other studies, such as pairwise F_ST (<0.02) among wild edible cockle using SNPs information (Vera et al., 2022) and among wild populations of Pacific oyster using allozymes and microsatellites markers (Appleyard and Ward, 2006). The PCA result showed that the WIM breeders formed a substructure with allele frequencies that differed from those of the other populations. This divergence may reduce the accuracy of EBV prediction in WIM populations and may be a cause of masking of the genetic architecture, as observed in a study of resistance to nodavirosis in sea bass (Griot et al., 2021). However, this clustering is inconsistent with the low mean F_ST values (<0.03), suggesting limited genetic differentiation and a degree of genetic relatedness among the populations. The two families whose offspring showed greater isolation from others in the WIM population may be due to the peculiar characteristics of the parents, which drive the first axis of the PCA.

Our study presents the first report of heritability estimates for resistance to V. splendidus experimental infection in M. edulis based on genome-wide SNPs. We observed moderate heritability for V. splendidus resistance (0.22–0.36), which are higher compared to our previous study using the same population. This increase may be attributed to the inclusion of a third experimental infection in this study, despite the overall lower mortality rate (Ajithkumar et al., 2024b). Disease resistance to pathogens in bivalves seems to be a heritable trait, with moderate to high heritability in oysters, clams, and abalone, ranging from 0.21 to 0.63 (Degremont et al., 2015; Brokordt et al., 2017; Smits et al., 2020). Studies on oysters have shown varying levels of heritability (h²: 0.09–0.54) against different Vibrio spp. pathogens at different life stages (Azema et al., 2017; Nordio et al., 2021; Zhai et al., 2021; Dietrich et al., 2022). When comparing heritability estimates across different methods, both GBLUP and ssGBLUP consistently showed higher heritability estimates than pedigree-based methods. This difference is likely due to the genomic relationship matrix constructed based on genome-wide SNPs information can capture both within and between-family genetic variance, whereas traditional pedigree-based selection only captures genetic variance between families (Boudry et al., 2021). The higher heritability was observed with GBLUP, which may be due to the selective genotyping approach (348 dead and 348 alive individuals genotyped), potentially biasing heritability estimates and limiting the capture of family genetic variation. However, the ssGBLUP approach, which includes both genotyped and non-genotyped individuals, provides a more precise heritability estimate by capturing family variation from both genotyped and phenotyped individuals. To date, numerous studies across aquaculture species have similarly demonstrated that GBLUP methods provide higher estimated heritability and greater accuracy compared to PBLUP (Tsai et al., 2015; Gutierrez et al., 2018). These results underscore the presence of genetic variation for resistance to V. splendidus in our mussel populations, and highlight significant opportunities for enhancing disease resistance through selective breeding programs, whether using pedigree-based or genomic selection strategies.

QTL detection in our populations posed challenges due to limited number of markers and individuals. Given the data in the current study do suggest a polygenic nature of resistance to V. splendidus, exploiting all markers to calculate genomic breeding values for resistance may be the most effective approach. Our association analyses suggest that resistance against V. splendidus exhibits a polygenic architecture without major QTLs. Similar findings have been reported for bacterial disease resistance in various aquaculture species including, Atlantic salmon Salmo salar (Correa et al., 2015), Coho salmon Oncorhynchus kisutch (Barría et al., 2018), Gilthead Sea Bream (Palaiokostas et al., 2016), European seabass (Oikonomou et al., 2022), and Pacific Oyster (Yang et al., 2022). For instance, a study on catfish identified four QTLs associated with columnaris resistance using a high-density SNP array (Geng et al., 2015), highlighting the importance of high-density SNP array for GWAS studies. Furthermore, multiple studies on rainbow trout concluded that a panel of approximately 5,000 evenly distributed SNPs across the genome was effective in identifying multiple QTLs for resistance to bacterial cold water disease due to the strong long-range LD in breeding populations (Vallejo et al., 2017; Vallejo et al., 2018; Vallejo et al., 2022). The SNP array used in this study had a limited number and quality of SNPs, and the mussel genome showed a very low level of LD. Our study used 2,204 SNPs, which may not provide sufficient coverage given the rapid LD decay and unevenly distributed SNPs across the genome, potentially leading to the omission of important QTLs. Future studies would benefit from incorporating high-quality SNP markers (>30K SNPs) and a larger number of phenotyped individuals (>1,000), which could significantly enhance the power and precision of QTL detection for V. splendidus resistance in blue mussels (Jones et al., 2013b; Barría et al., 2018). Additionally, developing an F2 generation by crossing parental lines with clearly contrasting phenotypic traits (resistant and susceptible) would provide an ideal basis for QTL mapping (Gu et al., 2011). In F2 populations, genotype segregation combined with divergent phenotypic differences provides strong conditions for QTL identification (Geng et al., 2015). Family-based mapping in these populations further strengthens QTL detection power by minimizing recombination between QTLs and linked markers, thereby improving the accuracy of marker-trait associations (Geng et al., 2015; Gonen et al., 2015). Moreover, family-based association mapping is particularly effective for detecting rare alleles associated with QTLs, which may be challenging to identify in other population structures (Geng et al., 2015).

The accuracy of genomic selection is affected by several factors, including the relationship between training and validation animals, sample size in the reference population, marker density, effective population size, LD structure, underlying trait architecture and heritability of the trait (Yáñez et al., 2023). Therefore, the lower prediction accuracies observed here may reflect the marker density or underlying genetic architecture of the trait. The choice of genomic selection model for breeding programs requires a prior understanding of the genetic architecture of the selected trait(s). In the current study on M. edulis populations, the genetic contribution to the observed variation in resistance to V. splendidus was largely polygenic in nature. For the improvement of polygenic traits, GBLUP is the most reliable model and typically provides the highest prediction accuracy for highly polygenic traits, while the Bayesian models are preferable for traits controlled by few large effect loci in genomic selection (Legarra et al., 2015; Yáñez et al., 2023).

In our study, genomic selection improves accuracy of up to 19% compared to pedigree selection. A key consideration for the commercial implementation of genomic selection in shellfish aquaculture is the high cost of genotyping. Reference population size and marker density are two key factors for effectively reducing the cost of genomic selection (Song et al., 2022). One strategy to enhance the economic viability of genomic selection is to use a low-density SNP panel, which can reduce costs while achieving prediction accuracy similar to that of medium-density panels (Kriaridou et al., 2020). The prediction accuracies for genomic models in our study ranged from 0.32 to 0.48 for resistance to V. splendidus (with SNP densities ranging from 500 to ∼3,400), whereas the accuracy of PBLUP was 0.36. This result is slightly lower than the ranges reported for disease-related traits in other bivalve species. For instance, genomic selection prediction accuracies from GBLUP models for resistance to Ostreid herpesvirus (OsHV-1-lvar) ranged from 0.68 to 0.76 in the Pacific oyster (Gutierrez et al., 2020). Prediction accuracies for growth-related traits using the GBLUP model in other bivalves are relatively similar, e.g., 0.52–0.73 in the Pacific oyster (Gutierrez et al., 2018; Jourdan et al., 2023), 0.67–0.79 in the Portuguese oyster (Vu et al., 2021), 0.25–0.48 in the Zhikong scallop (Wang et al., 2016) and >0.83 in European flat oyster (Penaloza et al., 2022). Other reports on genomic prediction accuracies for disease-related traits in finfish aquaculture species show the prediction accuracies as low as 0.21, reviewed in Houston et al. (2020). However, this result highlight that genomic selection is a useful approach to increase resistance to V. splendidus in our blue mussel populations.

Overall, our results showed that genomic methods predict better accuracy (25%–33%) for resistance to V. splendidus using ∼2000 SNPs in a family-based design compared to pedigree-based estimation. This indicates that substantial improvements in the rate of genetic gain can be achieved through genomics-based selection techniques. It also increases the feasibility of using a low-density SNP panel for genomic selection in mussel breeding for Vibrio resistance, as low-density genotyping is significantly more cost-effective than high-density SNP arrays. Furthermore, studies on disease resistance in the Pacific oyster, growth traits in the European flat oyster, and heat tolerance in the Pacific abalone have shown that low-density SNP panels of around 1,000–2000 SNPs can achieve EBV accuracies similar to those obtained with medium-density arrays (Gutierrez et al., 2020; Kriaridou et al., 2020; Liu et al., 2022; Penaloza et al., 2022). Similar findings in multiple aquatic species have shown that low-density panels can achieve higher accuracies than the pedigree-based approach, making them a viable alternative for identifying candidates with the highest genetic merit for complex traits such as growth and disease resistance (Kriaridou et al., 2020).

Although the mussel genome is 1.4 Gb in size, our study suggests that a relatively low number of genetic markers can still achieve high prediction accuracy, with a rapid LD decay observed across all populations. In this study, we highlight the possibility of reducing the genotyping costs associated with genomic prediction approaches, caution should be exercised regarding the smallest marker density. Our study found that using only 500 SNPs in the GBLUP model resulted in an estimated decrease in the accuracy of genomic breeding values (GEBVs) for resistance to V. splendidus by 11% compared to PBLUP. It is important to note that when both pedigree information and genotypes are available, using ssGBLUP is preferable, as it demonstrates superior accuracy compared to PBLUP. Furthermore, our findings emphasize that annotated SNPs on the M. edulis genome provided more information about the studied population and led to higher prediction accuracy than using all SNPs in either the GBLUP or ssGBLUP model. This difference could be due to even distribution of phenotypes among genotyped individuals (50% mortality), or unannotated markers may introduce noise, thereby affecting the accuracy of GEBV estimation. Further investigations using more SNPs (>10K SNPs), and larger reference population (>1,000 individuals) hold potential for genomic selection to further increase the prediction accuracy for host resistance to V. splendidus in farmed mussel populations.