© H.J. Alves Monteiro et al., Published by EDP Sciences 2022

1 Introduction

An important goal in endangered species conservation programs and the reintroduction of locally depleted or extirpated populations is preserving and maximising genetic diversity. Stock enhancement of marine invertebrate species has been reported to be successful in terms of boosting short- to medium-term census population size, e.g. in Japanese scallop,  Mizuhopecten yessoensis (Kitada 2020). However, in spite of potential effects on long-term success (Reed and Frankham, 2003; Weeks et al., 2011), there is generally little information about levels of genetic diversity represented in stocking material of marine invertebrates. The prediction is that using a (small) broodstock with limited genetic variation to replenish local populations will accelerate rates of inbreeding and loss of genetic variation, reducing the long-term fitness of a population (Ryman and Laikre, 1991). Modelling analyses demonstrate potentially severe effects of failing to maximise genetic diversity in stocking material used to re-establish or replenish local populations (Waples et al., 2016), and genetic properties should hence always be taken into consideration in programs for marine stock enhancement, not least in bivalve keystone species with strong small-to-large scale impact on nutrient cycling, habitat stabilisation, food-webs and providing ecosystems services. Although census population size, N, and level of genetic variation commonly are very large under wild conditions, the genetically effective size, N_e, of marine bivalve populations may be significantly reduced by strong reproductive skew among individuals in combination with recurring ‘sweepstakes’ effects (Hedgecock et al., 2007). Under a ‘sweepstakes’ scenario, the genetic variation propagated across generations is strongly affected by relatively few reproductive “winners” and many “losers”. This leads to a very small and varying N_e/N ratio (Hedgecock and Pudovkin, 2011), making it difficult to predict population levels of genetic variation from census size. Supplementation of a small, genetically depleted natural population by genetically diverse hatchery-reared spats (i.e., juvenile settled oysters) as a form of “genetic rescue” is generally expected to have beneficial effects (Frankham, 2015). However, to be successful, all types of supplementations require comprehensive guidelines for genetic management, which is lacking for most of the depleted bivalve species, including the European flat oyster (Ostrea edulis L.).

O. edulis, commonly called flat oyster or native European oyster, is a sessile bivalve mollusc in the Ostreidae family. The species is endemic to the European coasts and can live on various substrate types at depths down to ∼50 m. O. edulis beds provide numerous important ecosystem services and have been included in the OSPAR list of endangered and/or declining species since 2003 (https://www.ospar.org/documents?v=7183). This has resulted in several recent initiatives of conservation and restoration projects in Europe (Pogoda et al., 2019). O. edulis abundance decreased dramatically in the 20th century, so much so, that the concept of a shifting ecosystem baseline, also known for several important fisheries (Pauly, 1995), could be applied for the species (Pogoda et al., 2019). Many historical flat oyster populations are now extinct or highly depleted (Berghahn and Ruth, 2005; Lotze, 2005; Thurstan et al., 2013). The decline in abundance is the direct result of long-lasting overfishing coupled with, more recently, disease outbreaks (Grizel, 1985; Harrang, 2013; Harrang et al., 2015). In several places, the small remnants of O. edulis populations are fragmented and at risk of extinction due to inbreeding depression (Beck et al., 2011). Numerous, partially unregulated, translocations of flat oysters have occurred between European countries to support production (Bromley et al., 2016), which have increased the spread of protozoan parasites Bonamia ostreae and Marteilia refringens in natural populations (Comps et al., 1980; Sas et al., 2020). Hence O. edulis is at risk from genetic, biotic, as well as abiotic stressors (elevated sea temperature and oceans acidification, see Kamermans and Saurel, 2022), that could further reduce the numbers of last wild populations.

In response to dwindling oyster populations, the last decade saw revised restoration strategies for reintroducing the species in severely depleted sea beds using various spat production systems: hatchery with limited numbers of breeders per tank (n < 50), ponds with large numbers of breeders per tank (n > 100) and spat collectors in the wild (see Fariñas-Franco et al., 2018). Several knowledge gaps challenging restoration efforts have been identified with respect to reproductive biology and the transmission of genetic variation through hatchery production (Colsoul et al., 2021, 2020). For example, control of genetic diversity and inbreeding are difficult to achieve in hatchery operations due to potential uneven parental contribution. As a result, there is scarce information about how various hatchery procedures affect genetic variation among produced spats.

Parentage assignment with genetic markers can provide information on the reproductive skew and the magnitude of cross-fertilization in aquaculture broodstocks (Lallias et al., 2010a; Varney and Wilbur, 2020) and has also been used to demonstrate inbreeding depression (depressed growth) in response to increasing relatedness levels in hatchery broodstock selected for Bonamia resistance (Naciri-Graven et al., 1999). O. edulis parentage contribution for wild and hatchery-reared progeny has traditionally been assessed using a limited number (n < 6) of genetic markers that tend to show high frequencies of non-amplifying alleles (‘null alleles’; Hedgecock et al., 2007; Lallias et al., 2010a, 2010b). The presence of ‘null alleles' is a commonly reported challenge in bivalves (Gonzalez-Wanguemert et al., 2014; Hedgecock et al., 2007; Reece, 2004), and could introduce errors into empirical assessments of biological parentage and mating systems also in O. edulis by means of incorrect parentage exclusions (e.g. incur substantial biases on average exclusion probabilities; see Dakin and Avise, 2004). However, even with recently improved access to genomic resources (Gutierrez et al., 2017), microsatellite marker-based parental reconstruction is still the most commonly applied method, being cheap and quick, with numerous microsatellites markers available (Lallias et al., 2009; Launey, 2002; Naciri et al., 1995; Vera et al., 2015).

This study aims to address identified knowledge gaps related to transfer of genetic variation in flat oyster hatchery production for restocking and aquaculture (Colsoul et al., 2021, 2020; Nascimento‐Schulze et al., 2021) by experimentally assessing parentage under controlled hatchery conditions. In hatchery production, reproductive success is expected to be equal among broodstock individuals through common hatchery conditioning (Reynaga-Franco et al., 2020). However, the validity of this assumption is likely to vary in connection with several factors. As part of the new strategy implemented in the Danish O. edulis hatchery to produce parasite (Bonamia) free spat and preserve high genetic diversity, we developed and tested an effective and easily applicable method to assess reproductive skew and genetic diversity in O. edulis hatchery produced spat. Four broodstocks of varying number were constructed using wild adult oysters collected from two localities in the Limfjorden (Denmark). We applied microsatellite marker-based genetic parentage analysis in combination with genetic diversity analysis based on a suite of estimators (allelic richness, heterozygosity, relatedness, N_e, and effective number of breeders N_b) to determine genetic parameters in offspring produced in five separate swarming events. Secondly, we compared genetic diversity in hatchery-produced spat with estimates from wild population samples, including from natural settlement on a spat collector. Finally, we use these results to provide recommendations for hatchery practice.

2 Materials and methods

2.1 Biological material and experiments.

The Danish Shellfish Centre at DTU Aqua dredged flat oysters for broodstock in 2020 at two locations, Nissum and Løgstør, considered to belong to a single natural population (Nielsen and Petersen, 2019) in the Limfjorden, Denmark (Fig. 1a). Fifty-two adult flat oysters were selected for broodstock by phenotypic features indicative of maturity and general good quality (based on length, width, total fresh weight (Supplementary notes: Hatchery's procedures, broodstock conditioning) and of being pathogen-free (Bonamia ostreae and Marteilia refringens). In the hatchery, they were distributed among four separate experimental reproduction tanks as follows:

• a large broodstock tank: 29 parents from Løgstør (‘Large’), 400L April 2020.

• a small broodstock Tank 1: 6 parents from Løgstør (‘Small 1’), 24L April 2020.

• a small broodstock Tank 2: 7 parents from Løgstør (‘Small 2’), 24L April 2020.

• a small broodstock Tank 3: 10 parents from Nissum (‘Small 3’), 100L March 2020.

Small tanks were named taking into consideration the broodstock size: from ‘Small 1’ with the least individuals to ‘Small 3’ with the most individuals. This design allowed a comparison between parentage contribution patterns within and among relatively large (‘Large’) and relatively small (‘Small 1–3’) broodstocks. For O. edulis, the release of male and female gametes is defined as the spawning event. The subsequent fertilization is internal: oocytes positioned in the female pallial cavity are fertilized by inhaled spermatozoa (His et al., 1999). After fertilisation, embryos are incubated internally, before being released from the female oyster as planktonic larvae 6–10 days later. The term “spawning” is used for gamete release, and larval releases are commonly referred to as “swarming” events (Colsoul et al., 2021; His et al., 1999). Five swarming events occurred and were sampled for analysis: respectively, two in ‘Large’ (i.e., “Large sp1” in June and “Large sp2” in July 2020), and one in each of ‘Small 1–3’ (respectively in July 2021, June 2020 and May 2020, Fig. 1b and c, Supp. Fig. 1). All released larvae were collected using a mesh net of 80–180 μm (see details Supp. notes: Hatchery's procedures, broodstock conditioning). Subsamples of each swarming event were taken to produce spat (juveniles that have reached the substrate settlement stage). All broodstock and random subsets of ∼48 offspring per swarming event were sampled for genotyping (Tab. 1). We inferred an instance of sex reversion between two separate swarming events in the Large tank (see Results). To explore temporal development in parentage patterns further we sampled twice the number of spats (n = 99) from the second swarming event (Tab. 1). Tissue samples consisted of minute gill biopsies for adults, and total tissue for spat (directly stored in absolute ethanol after sampling until processing). To examine the developed genotyping approach for analyses of larvae (individuals of approximately 120 μm, i.e., small input volume tissues), we also analysed material from two additional hatchery setups from 2019, where 24 broodstock from Nissum and small larvae collected on their first day of swarming (24 larvae for both setups, Supp. Tab. 1) were sampled and genotyped for parentage assignment.

To compare genetic parameters between hatchery produced spats and the wild broodstock donor population, and hereunder from a natural settling event, we analysed three additional samples: 66 and 42 adult oysters of different size classes were collected from, respectively, Løgstør and Nissum, and from Nissum, 30 newly settled spats were collected using a traditional “Chinese hats” spat collector (van den Brink et al., 2020) deployed by Oyster Boat ApS company. The devices were deployed in June 2018, and spats settled on it were collected for genetic analysis. We hypothesized that these spats were the result of a natural swarming season at this given location, spreading from mid-June to mid-July 2018.

Fig. 1 Hatchery setting and breeding pairs. (a) Flat oyster sampling locations in the Limfjorden connected to the North Sea in west and to the Kattegat in east. Coloured symbols indicate locations for collections of broodstock at Løgstør (green), Nissum (blue), and for the natural settling event (‘Natural settlement’, purple). (b and c) Experimental set-up in four parental broodstock tanks (‘Small 1-3’ in panel (b) and ‘Large’ with two separates swarming events analysed in panel (c)), and offspring from the five analysed swarming events. For individual tanks, hexagons show broodstock sizes and origins (green = Løgstør, blue = Nissum). For individual swarming events, each coloured segment of the chord diagram represents a paternal full-sib family, and the thickness of the chord represents the number of offspring identified, as shown by peripheral tick-marks. Individual swarming collections all had one single dam, the label of which is in bold. Two instances of sex reversion (one dam to sire and one sire to dam) in parental breeders are indicated by the pink label symbol “R”. Sires contributing to two consecutive swarming events are indicated by the blue label “C”. Tags inside chord diagrams give the total number of observed breeders (dam and sires) as well as the number of effective breeders Nb.

2.2 Microsatellite genotyping

Genomic DNA was extracted from individuals using 15 mg of gill tissue from adults, entire larvae, and entire soft tissue (without the shell) for the spat. Except for larvae, we used the E.Z.N.A.^® DNA Kit (Omega Bio-Tek, Norcross, GA, USA) for DNA extractions, diluting DNA products at 50 µl. Given their small volume, individual larval samples were separated under a dissection microscope and placed in individual tubes. Several DNA-extraction protocols (e.g., using short-term exposure to liquid nitrogen to break larval membranes, centrifugation with ceramic beads, decreased elution volume) were compared to produce adequate DNA yield (details available upon request). Maximal DNA concentrations (5.73–452 ng/μl) were obtained with the NucleoSpin Tissue XS kit (Macherey-Nagel, Düren, Germany) following the manufacturer's protocol. Nineteen species-specific microsatellites markers were selected from the literature (Lallias et al., 2009; Launey, 2002; Naciri et al., 1995; Vera et al., 2015). One additional marker was developed for the current study, designed on the first scaffold of the species' unpublished genome (Pers. Comm. A. Tanguy) using the Primer-Blast (https://www.ncbi.nlm.nih.gov/tools/primer-blast/) and Primer3 tools (https://primer3.ut.ee/). The twenty markers were first tested on subsets of adults, larvae, and spats using a Qiagen® Multiplex PCR protocol for amplification of microsatellite loci (Supp. Tab. 2). Polymerase chain reactions were conducted in volumes of 12.5 µl, with an annealing temperature of 57 °C for all loci. Three markers exhibited evidence of large heterozygote deficiencies during screening and were discarded. All individuals (collections from the wild, hatchery broodstock, larvae, and spats) were then genotyped for the remaining seventeen microsatellites contained in four multiplexes. PCR products were analysed utilizing GeneScan® 500 Liz Size Standard in a SeqStudio Genetic Analyser (Applied Biosystems, Foster City, CA, USA) and genotypes were called in GeneMapper v2.4. Each microsatellite was calibrated and tested on ∼50 of the wild oyster collections (Nissum and Løgstør). To assess the reliability of genotyping, five per cent of the samples, covering 11 broodstock, and 13 spats, including spats from parental sex reversions (Results), were PCR amplified and genotyped twice. Thirteen hatchery-reared spats with more than three missing genotypes across the 17 loci were discarded from the analysis (2.6% of 492 genotyped individuals).

2.3 Genetic diversity and differentiation

Number of alleles (Na), percentage missing genotypes and test for departures from Hardy-Weinberg equilibrium in each sample (1000 replicates for the Monte Carlo procedure) were assessed using GenAlEx for each sample (Smouse and Peakall, 2012). The R package hierfstat (version 0.5–7, Goudet, 2005) was used to estimate per sample the observed heterozygosity (Ho) and the inbreeding coefficient (Fis) (Nei's method done with the function basic.stats). The hierfstat function allelic.richness was used to calculate the mean absolute allelic richness (Ar) per sample.

Null allele frequency per locus and sample was estimated with both ML-NULLFREQ (Kalinowski and Taper 2006) and MICRO-CHECKER (Van Oosterhout et al., 2004). The two methods gave similar results and only the null alleles frequency values from MICRO-CHECKER are reported for each locus (Supp. Tab. 3). In MICRO-CHECKER, tests for scoring error and presence of null allele were performed by randomizing genotype data a thousand times and corrected for repeated sampling for each locus (following Bonferroni). Estimates of genetic differentiation (F_st, Weir and Cockerham, 1984) and their confidence intervals (0.025 and 0.975 quantiles) were generated for each pair of samples using a procedure with 10,000 bootstraps over loci with hierfstat function boot.ppfst. We considered pairwise F_st values different from zero if the lower bound of the confidence interval was larger than 0.

2.4 Parental reconstruction

Following Jones et al., (2010), parental reconstruction augmented by a full probability model (maximum likelihood) was chosen for parentage assignment. The genotype data were analysed both with the user-friendly (GUI) version of COLONY 2.0 (version 2.0.6.626, Wang, 2004) and checked for consistency using Cervus (version 3.0.728, (Kalinowski et al., 2007) to obtain parent pair assignments for each hatchery-reared spat and larvae. Scoring error rate per locus was set to 0.5% but was increased to 0.7% for loci presenting non-negligible null allele frequency estimates (see Results). Identified parent pairs were compared between methods, as well as the degree of confidence obtained. COLONY settings included a polygamous mating system, strong sib-ship prior, and a 95% probability the true sire and dam to be included among parental candidates. For Cervus, a simulation was performed using genotypes from 108 wild oysters (Nissum and Løgstør, combined) to estimate allele frequencies. Then, to assign parents at a designated level of confidence in runs of parentage analysis for each tank, we used the “parent pair analysis where the sexes of candidate parents are unknown” module. We followed the author's recommendation of only accepting assignments made at the 95% or greater confidence level. The assignments made with COLONY were deemed consistent with Cervus when the same parent-pair was identified with a probability of >.80, following Kalinowski et al., 2007. O. edulis is an alternating hermaphrodite (mainly protandrous), with highly flexible reproductive tactics. Both male and female gametes can be observed simultaneously within individuals (Da Silva et al., 2009; Maneiro et al., 2020). We were not able to non-invasively identify the functional sex of broodstock for individual spawning events. Instead, we inferred dam, versus sire, identity using parentage analysis results in combination with the parsimonious principle that the single parent producing all offspring within an individual swarm (see Results) was the dam, and the other contributing individual of the parent-pair, was the sire. We calculated the numbers of progeny produced by each dam-sire pair per swarming event and used chord diagrams (R package “circlize”: Gu et al., 2014) to illustrate results. To do so, we created contingency matrices (one for each swarming event) of the number of progeny produced by each parent-pair, tracing numbers of spat by each inferred sire to their inferred dam (Figs. 1b and 1c).

We assessed the number of offspring needed to produce reliable estimates of the total number of sib-ship families produced per swarming event by plotting the cumulative number of identified full-sib families as a function of numbers of analysed spats (N), using parentage data from the Large broodstock, which contained the largest number of potential parents (n = 29). The value of N, where the curve flattened, gave a visual indication of whether enough offspring were genotyped to encounter most sib-ship families produced in one swarming event (Supp. Fig. 2).

2.5 Effective number of breeders

The effective number of breeding males per swarming event was calculated with a method modified from Ardren and Kapuscinski, (2002):

$$N_{bm}=(N_m{{\displaystyle \overline{k}}}_m-1)/[{{\displaystyle \overline{k}}}_m+({\sigma }_{km}^2/{{\displaystyle \overline{k}}}_m)]\text{,}$$

where N_bm is the effective number of males in a swarming event, k represents the number of progeny produced by an individual and N_m represent the number of contributing males in a swarming event (Figs. 1b and 1c). As all swarming events were inferred to be produced by a single dam (see Results), the effective number of breeders, N_b, was estimated as N_bm + 1.

2.6 Relatedness

Sample specific relatedness estimates were determined using COANCESTRY (version 1.0.1.1021 (Wang, 2011) and the Lynch and Ritland estimator (Lynch and Ritland, 1999). First, the allele frequencies in wild collections (Nissum and Løgstør, n = 108) were used to bound the estimation of relatedness parameters for each hatchery swarming event and the “Natural settlement” samples. If an allele present in a swarming event was not detected in the wild, it was given an arbitrary frequency of 0.001. Mean relatedness coefficients per sample and their 95% confidence intervals were estimated with the 1000 bootstraps procedure. Pairwise-Wilcoxon signed rank tests were performed for each sample pair to test whether the mean relatedness differed. For each swarming event, marker-based relatedness estimates were compared with parentage based average relatedness estimates, assuming that full-sibs and half-sibs were related by 0.50 and 0.25, respectively.

2.7 Effective population size

Contemporary effective population size (N_e) was estimated using the single sample linkage disequilibrium method (Hill, 1981; Waples, 2006; Waples and Do, 2010) as implemented in “NeEstimator (version 2.1, Do et al., 2014) under a random mating model. We chose the linkage disequilibrium model to infer effective population size and excluded singleton alleles from the dataset (lowest allele frequency used: 0.01).

3 Results

3.1 Genotyping results and null alleles

The attempt to genotype individuals before they reach the spat stage, i.e. using samples of hatchery produced larvae from 2019 (Supp. Tab. 1), was not successful as all 48 larval samples consistently displayed high degrees of missing data (average 8.6% missing loci). Further, as preliminary parental assignment analyses indicated low classification success, indicative of scoring error, all larvae were discarded from the parentage analysis.

With larvae removed, we genotyped a total of 492 individuals, 332 spats (30 from the natural settlement sample and 302 hatchery-reared spats) and 160 adults (108 from the wild samples and 52 from the four broodstocks). All adults passed the genotyping success criteria, however 13 hatchery-reared spats had more than 17.6% of missing data and were therefore excluded, resulting in 479 genotypes in total (Tab. 1). The proportion of missing data per locus in the retained individuals ranged from 0 to 33%, with an average missing data per locus of 1.76% for all samples cumulated. Evidence of null alleles was detected in four loci (Oed 165, OE-27, OEe-45 and OEe-46) out of the seventeen loci (Supp. Tab. 3), as indicated by significant heterozygote deficiencies and/or detection of null alleles with MICROCHECKER in the wild populations. The loci were retained but parentage analyses were done by increasing the expected genotype scoring error value from 0.5% to 0.7% for these four loci, which is the recommended approach under, respectively, low and moderate to high frequencies of null alleles (J. Wang Pers. Comm.).

3.2 Genetic diversity

Mean allelic richness, Ar, estimates were similar for wild samples (5.265 and 5.253 for Løgstør and Nissum respectively), as was mean observed heterozygosity (Ho, 0.796 and 0.793, respectively, Tab. 1). A statistically significant reduction in Ar was found between broodstock and spats in all hatchery setups (pairwise Wilcoxon signed-rank test p < 0.05). Inbreeding coefficient, F_is, values were negative for all hatchery-reared collections, ranging from −0.245 to −0.094, while broodstocks showed estimates ranging from 0.026 to 0.067.

3.3 Genetic differentiation

All pairwise F_st estimates between wild collections (Nissum, Løgstør, and natural settlement) were low, ranging from −0.004 to 0.002, corroborating the assumption that they were representative of a single local population (Supp. Fig. 3, Supp. Tab. 4). Pairwise F_st's between hatchery-reared spat swarms ranged from 0.004 to 0.184. No significant differences between the wild population (Løgstør, Nissum) and the broodstocks was detected (unbounded lower confidence intervals, Supp. Tab. 4).

3.4 Parentage

Evaluation of the total number of sib-ships identified as a function of spat sample size showed that detection of all, or at least the majority, of full-sib families in each swarming event was attainable with genotyping of 30 offspring (Supp. Fig. 2). It indicates that our applied sample size of spat genotyped is sufficiently large to identify all the breeding pairs involved in each swarming event. All 289 spats were consistently assigned to two parents at p > 0.999 in both COLONY and Cervus (Supp. Tab. 5). Each swarming event was inferred to be produced by a single dam and multiple sires (Figs. 1b and 1c). The numbers of progeny varied between 1 and 44 in sires, and 44 and 99 in dams. Individual males exhibited strong reproductive skew, siring between 1 and 92% of the analysed offspring per swarming event. Across tanks, 65% of broodstock was inferred to produce at least one progeny. Within tanks, 30–66% of broodstock did not contribute to sampled spat. The two consecutive swarming events in the Large broodstock were produced by two different dams, and sired by, respectively, sixteen and nine males: of which four sired offspring in both events (Fig. 1c). Two instances of sex reversal were indicated in the Large broodstock tank (Fig. 1c). Thus, the dam of the first swarming event was inferred to sire one offspring in the second swarming event (i.e., female reverting to male), and a sire of one offspring in the first swarming event became the dam of the second swarming event (i.e., male reverting to female), after a period of 1 month in between the two swarming events. In these two offspring inferred to be the results of sex reversed parents, we took care to reassess hatchery sampling and laboratory labelling setups, as well as re-genotyped offspring (see above) to control for potential error. Parentage assignment probabilities were >0.99 for each of the two identified parents and <0.00 for all other broodstock, for both spat. When inspecting genotype data for these two spats and all potential parental candidates we found that mis-assignment of the correct parent would require an unlikely high occurrence of 9‒15 point mutations or genotyping errors (average 11.7) per individual. In comparison, there were no genotype mismatches between the inferred sex reversed parents and their offspring. In the two small tanks with less than ten broodstock, respectively, two and one males sired the majority of spats, and the tank with ten broodstock showed a higher number of sires (six) (Fig. 1b).

3.5 Effective number of breeders

Estimates of the effective number of sires, N_bm, ranged almost an order of magnitude (0.91–9.03) across swarming events, and were on average approximately half (55%) of the census number of contributing sires (Figs. 1b and 1c). When comparing total N_b with total broodstock size, the average ratio came out at 31% (ranging 14–46%) across swarming events. Interestingly the second swarming event of the Large broodstock (the “Large sp2”) had a lower number of effective breeders than the ‘Small 3’ broodstock's only swarming event, respectively, 4.08 and 6.48 effective breeders (It is worth noting that ‘Small 3’ broodstock's individuals were older with total fresh weight and shell length bigger than the other broodstocks, see Supp notes: Hatchery's procedures, broodstock conditioning).

3.6 Relatedness

Marker-based relatedness estimates were on average high (>0.25) in all hatchery spat collections, and low and statistically non-significant from zero in all wild collections (Fig. 2). Relatedness estimated in the “Natural settlement” was slightly higher (0.02) than in the two other wild samples (statistically significant different, pairwise Wilcoxon signed-rank test p < 0.05). A post hoc sibship analysis with COLONY, and supported by COANCESTRY, identified five individuals as sibs, supporting that the natural settlement sample also included some related individuals (Fig. 2). All hatchery-produced spat samples were significantly different, except ‘Large sp1’ and ‘Small 3’ spats (Wilcoxon signed-rank test p > 0.05). Sib-ship based average relatedness estimates overall corresponded with marker-based estimates (Fig. 2).

Fig. 2 Mean Relatedness for each sample. Lynch and Ritland’s (1999) relatedness estimator calculated in COANCESTRY (Wang, 2011). Each point represents the relatedness of a given dyad. Average relatedness estimates per sample are connected for easy overview with blue dots and lines. Parentage-based relatedness for each hatchery-reared spat sample is depicted with a red triangle. Negative relatedness values are consequences of unbiased estimation of relatedness and represent a relatedness value of 0. Point sizes for outliers are increased. The lines dividing each box into two parts are medians, ends of the box are the upper and lower quartiles.

3.7 Effective population size

N_e point estimates in hatchery spats ranged from 3.1 to 26.1 per swarming event and were all bounded at the 95% confidence limits (95% confidence interval range 2.8–28.6, Tab. 1). Estimates for the two wild collections of adults came out as unbounded for both Løgstør and Nissum; indicative of a very large N_e. The N_e estimate for the natural settlement was lower, estimated at 189.3, and with a relatively broad, but bounded, 95% confidence interval (121.0–412.2).

4 Discussion

Our study represents one of still few examples of genetic monitoring of hatchery-reared spats in O. edulis, providing insights into broodstock parentage contributions with unprecedented statistical accuracy. Prior parentage analyses in the species reported incomplete parentage assignment, likely due to insufficient marker resolution, coupled with technical error due to inconsistent amplification (null alleles) in microsatellite markers (Hedgecock et al., 2007; Lallias et al., 2010b). Our microsatellite marker assay can likely be successfully applied across O. edulis management and aquaculture scenarios, although it is noted that null-allele frequencies could vary among populations and strains (Van Oosterhout et al., 1991). Our results strengthen knowledge about transfer of genetic variation from broodstock to spat sources produced under standard hatchery conditions. Across swarming events we found that the realised genetic variation in produced spat was consistently lower than expected from the census number of broodstock. Thus, average N_b/N ratio was 0.31, and as low as 0.14 in one swarming event with 29 broodstock individuals in the tank. Except for average heterozygosity, all estimators of genetic diversity supported the notion that hatchery-reared spats showed reduced variation in comparison with their parental broodstock source. Although specific hatchery conditions, hereunder a suite of biotic and abiotic factors, likely have effect on reproductive outcomes, our results are in line with results from other genetic parentage studies testing natural spawning behaviour under hatchery conditions. Thus, Lallias et al. (2010b) used a single large (n = 62) broodstock to test parentage in six swarming events produced over 15 days. There, the majority of broodstock (∼70%) produced at least one offspring. However, individual swarming events exhibited N_b/N ratios similar to ours (average ∼ 0.28, range 0.19–0.34). In our study, one broodstock (‘Large’) with 29 individuals produced two swarming events separated by one month, showing highly heterogeneous outcomes in terms of reproductive skew (respective N_b estimated at 10 and 4). Moreover, a quarter of the broodstock sired offspring in both swarming events. This shows that although spat production based on multiple swarming events from the same broodstock is likely to increase genetic variation, there is still a degree of ‘genetic redundancy’ that would not have been experienced with separate broodstocks. Yet, it is general hatchery practice to produce spats from several swarming events from the same brood tank, and our results indicate that this practice may increase total genetic diversity through the facilitation of broader parental contribution. As expected, brood tanks with lower broodstock size (n ≤ 10) generally produced spats with less total genetic variation and higher levels of relatedness. Nonetheless, for the ‘Small 3’ swarming event with ten broodstock, both N_e and N_b estimates exceeded those of ‘Large sp2’. The ‘Small 3’ broodstock was composed of larger (on average 11% longer) and heavier (on average 45% heavier) oysters compared to the other broodstocks, which may have increased their individual fecundity and lowered reproductive variance (see details Supp. notes: Hatchery's procedures, broodstock conditioning). Whatever the cause, results underline the strong parentage variation and show that without genetic monitoring it is difficult to predict the amount of genetic variation represented in individual swarming events.

Our analyses supported that the two O. edulis samples collected across ∼100 km in the Limfjorden belong to the same wild population. Although census population sizes show strong fluctuations in this area (Nielsen and Petersen, 2019), our study showed evidence of very large effective population sizes, indicative of lack of severe local sweepstakes effects, and no imminent risk of inbreeding depression. Changing habitat conditions were suggested to explain decreasing N_e's in local populations of Pacific Oyster Crassostrea gigas (Hedgecock and Pudovkin, 2011). From an applied point of view, this indicates that collecting O. edulis in the Limfjorden for broodstock composition in the hatchery can contribute a high genetic diversity that can support the new strategy of keeping a Bonamia-free broodstock. However, our results also show that Ryman-Laikre effects are likely to affect natural populations through the skewed parental contributions of hatchery reared offspring if released back into the wild, and that these effects will be amplified if the same broodstock is used through consecutive breeding seasons. Ideally, new broodstock should be collected from the wild each breeding season to increase diversity and limit adaptation to hatchery environments (Lynch and O'Hely, 2001). However, for practical reasons, this may be difficult in a situation with Bonamia infected natural populations where extensive screening and quarantine procedures are needed to secure disease-free broodstock. Consequently, it is recommended to maximize diversity in the Bonamia-free broodstock through introducing new genetic diversity from the wild on a regular basis, in combination with monitoring of genetic diversity and relatedness in the hatchery to limit subsequent negative effects in natural populations.

An N_e of 50 (or 100) is considered the minimum threshold for a cultured population, in order to minimise immediate risk of inbreeding depression (Frankham et al., 2014). In our case, all individual swarming events showed N_e << 50. Thus, in a restocking context, the estimated number of swarming events required for directed conservation releases would amount to either three tanks with ∼30 broodstock, or ten tanks with ∼ ten broodstock, to meet a minimum criterion of N_e = 50. To minimise reproductive skew and maximize N_e it is always preferable to employ the largest possible number of broodstock tanks. However, even if practical considerations limit tank number, broodstock management should incorporate strategies to maximize N_e in stocking material to at least meet these minimum criteria. There is increased focus on the importance of expanding marine conservation efforts to include considerations of genetic resources, rather than simply focusing on maximising production (Grant et al., 2017). Throughout Europe, flat oyster brood tanks in hatcheries may vary in terms of individuals per broodstock. Furthermore, broodstock individuals may originate from a genetically diverse population or an inbred one. Therefore, the use of a high-powered tool to estimate parental contributions and population effective size in the hatcheries is needed to inform practices on a case-by-case basis.

Interestingly, the N_e estimate for the ‘Natural settlement’ sample which was assumed to represent a discrete local swarming season, returned a bounded result that was more than one order of magnitude higher than estimates from hatchery samples. Spat collectors were placed in an open bay and although exposed to natural recruitment over a whole season, 15% of the sampled spats were identified to be siblings. This result attests both to the hugely different levels of genetic variation represented in a natural versus a hatchery swarming event, but also highlights that genetic variation associated with local natural settling events only describes a fraction of the genetic variation within the total population. These results further demonstrate that our method provides a tool to monitor potential local recruitment processes and decreases in effective population sizes in the wild, as well as in a hatchery context.

Two hatchery swarming events separated by one month, showed two instances of sex reversion. Sex change under hatchery settings has been observed at relatively low frequency under pair-mating breeding designs (Toro and Newkirk, 1990). The sex ratio of a sex-changing species is expected to be biased towards the ‘first sex’ (i.e. male in the protandrous O. edulis) leading to a reduction in N_e in relation to N (Coscia et al., 2016). Two broodstock out of 29 (∼7%) in a tank were inferred to reverse sex and we found no reason to suspect parentage assignment error to have affected this estimate. Sex ratio dynamics could thus have strong impact on parentage patterns under hatchery conditions, but we emphasize that further studies are needed to determine the impact of sex change on genetic diversity in oyster spat production.

There is a consensus about the necessity for active and dynamic management actions in order to restore this once abundant species in European coastal waters (Matthiessen, 2008; Olsen, 1883; Thurstan et al., 2013), with the aim to simultaneously safe-guard ecosystem function and to improve access to low-carbon emission food resources. Regional flat oyster conservation projects are thus emerging throughout Europe (Colsoul et al., 2020; Fariñas-Franco et al., 2018; Pogoda et al., 2019; van den Brink et al., 2020; zu Ermgassen et al., 2020) strengthening the need for comprehensive information on the local genetic diversity of remnant populations, as well as for genetic monitoring of spat production methods. Loss of genetic diversity in O. edulis hatchery strains is an inherent consequence of a finite broodstock and increasing genetic diversity of spats is a central focus for restoration. Our results are aligned with this expectation and demonstrate that Ryman-Laikre effects (Ryman and Laikre, 1991; Waples et al., 2016) are likely to represent a real challenge in flat oyster restoration and conservation programs. While the lowest number of potential breeders did produce the lowest diversity in offspring, it is also clear that it may be difficult to predict the outcome of even controlled breeding set-ups in a hatchery, and that monitoring of genetic diversity may thus provide useful information for broodstock management. Here, we developed an easily applicable method to accurately estimate parentage, relatedness, and genetic variation in both wild and hatchery scenarios as a useful tool to control and implement new protocols in the hatchery.

Data Availability

Raw genotypes can be downloaded from: https://doi.org/10.5281/zenodo.6826349. R scripts are available upon request.

Funding

This project was funded by the European Regional Development Fund (Interreg, MarGen_II Project, #175806)

Competing Interests

The authors declare that they have no conflict of interest.

Author contributions

J.H.H., D.B. and C.S. conceived the project and obtained financial support. J.H.H., D.B., H.J.A.M. and C.S. designed the study. J.H.H., D.B., H.J.A.M. and C.S. led the project. C.S. collected broodstock and spats as well as conducted the experiment in the hatchery. H.J.A.M. performed the vast majority of DNA extraction and M.J. performed the majority of the genotyping. H.J.A.M. and D.B. conducted the bioinformatic analyses and all authors interpreted the results. H.J.A.M. and D.B. wrote the first draft of the manuscript with input from C.S. and J.H.H. All authors critically reviewed and approved the final manuscript.

Supplementary Material

Supplementary Table 1. The inferred parent pairs for parental assignements using larval material in 2019.

Supplementary Table 2. Description and diversity statistics summary, by locus, for wild adults flat oysters from Limfjorden.

Supplementary Table 3. Description and diversity statistics summary, by locus, for wild adults flat oysters from Limfjorden.

Supplementary Table 4. Pairwise (Fst) confidence intervals.

Supplementary Table 5. The inferred parent pairs for each offspring with their probabilities are obtained from the full likelihood method in COLONY.

Supplementary Figure 1. Summary of the swarming events in the hatchery.

Supplementary Figure 2. Cumulative curve. Offspring's genotypes needed to identify all the sib-ship families in a swarming event.

Supplementary Figure 3. (Fst) Heatmap. Samples are represented by each column, whereas pairwise (Fst)values for each pair of sample (Weir and Cockerham) are represented by the cells.

Acknowledgements

We would like to thank our local lab manager Dorte Meldrup for the prompt support during the execution of the project. We thank Carl St. John for his help with the figures. We are grateful for the technical support in the hatchery by Pascal Barreau and Lene Laursen Christensen. Thanks to Dr. Lone Madsen and Karí Mouritsen help regarding the DNA extraction of some broodstock. Helpful comments have been provided by Dr. Matthew P. Hare and Dr. Pierre Moënne-Loccoz. We would like to thank Dr. Arnaud Tanguy for allowing us an early access to the Ostrea edulis genome assembly. Furthermore, we would like to thank the members of the Native Oyster Restoration Alliance (NORA) and the MarGen_II project (European Regional Development Fund, Interreg) for their valuable discussions during the project.

References

All Tables

All Figures

Fig. 1 Hatchery setting and breeding pairs. (a) Flat oyster sampling locations in the Limfjorden connected to the North Sea in west and to the Kattegat in east. Coloured symbols indicate locations for collections of broodstock at Løgstør (green), Nissum (blue), and for the natural settling event (‘Natural settlement’, purple). (b and c) Experimental set-up in four parental broodstock tanks (‘Small 1-3’ in panel (b) and ‘Large’ with two separates swarming events analysed in panel (c)), and offspring from the five analysed swarming events. For individual tanks, hexagons show broodstock sizes and origins (green = Løgstør, blue = Nissum). For individual swarming events, each coloured segment of the chord diagram represents a paternal full-sib family, and the thickness of the chord represents the number of offspring identified, as shown by peripheral tick-marks. Individual swarming collections all had one single dam, the label of which is in bold. Two instances of sex reversion (one dam to sire and one sire to dam) in parental breeders are indicated by the pink label symbol “R”. Sires contributing to two consecutive swarming events are indicated by the blue label “C”. Tags inside chord diagrams give the total number of observed breeders (dam and sires) as well as the number of effective breeders Nb.

In the text

Fig. 2 Mean Relatedness for each sample. Lynch and Ritland’s (1999) relatedness estimator calculated in COANCESTRY (Wang, 2011). Each point represents the relatedness of a given dyad. Average relatedness estimates per sample are connected for easy overview with blue dots and lines. Parentage-based relatedness for each hatchery-reared spat sample is depicted with a red triangle. Negative relatedness values are consequences of unbiased estimation of relatedness and represent a relatedness value of 0. Point sizes for outliers are increased. The lines dividing each box into two parts are medians, ends of the box are the upper and lower quartiles.

In the text