The Effect of the Deformity Genetic Background of the Breeders on the Spawning Quality of Gilthead Seabream (Sparus aurata L.)

Fish egg quality is strongly related with the ability of the egg to be fertilized and develop a normal embryo with good survival and a lack of abnormalities. Large variations in the spawning quantity or quality impact directly in the competitiveness and sustainability of hatcheries, which create an overly large broodstock in order to satisfy the on-growing companies’ demand for undeformed fry. The present study reports, for the first time in relation to gilthead seabream, the effect of the genetic background of breeders for presence or absence of deformity on their spawning quality and the importance of considering this when creating broodstock. The spawning quality of crosses of breeders with genetic background for presence or absence of deformity (EBVdef), were evaluated during a whole spawning season, through study of the following traits: oocyte yield, fertilization rate, viability rate, hatching rate, larval survival rate, fertilized eggs, viable eggs, hatched eggs, and number of alive larvae. Breeders with a genetic background for deformity and a normal phenotype had shorter spawning periods, lower oocyte yield and, consequently, produced a lower number of alive larvae. In these two traits, the genetic background of breeders was of greater importance during intermediate spawning periods, when spawning is generally considered optimal for the industry, while environmental factors were more important at the beginning and end of the spawning season. In conclusion, these results demonstrate the importance of controlling the breeders’ genetics when creating broodstock.


INTRODUCTION
The gilthead seabream (Sparus aurata L.) is a demersal and eurythermal species with a large reproductive capacity. It is found in subtropical coastal areas, from the Eastern Atlantic to the Mediterranean (Froese and Pauly, 2019). In the wild spawning starts in October in the Gulf of Cadiz or December-January in eastern Mediterranean and lasts for 3-4 months (Arias, 1980;Kissil et al., 2001). Every female releases over two million eggs per kilogram they weigh (APROMAR, 2017). In aquaculture, larval rearing of the gilthead seabream is carried out at 16-24 • C (Tandler et al., 1989;Koumoundouros et al., 1997a,b;Shields, 2001), an optimal temperature for the hatching rate and larval survival rate at the yolk sac resorption stage (Polo et al., 1991). The wide natural distribution and environmental conditions experienced by the gilthead seabream have made it one of the most important species in the European marine fish farming and production is still expanding worldwide. In a European context, gilthead seabream was fifth in terms of annual production, at 91,964 tons, and fourth in terms of economic value (APROMAR, 2020).
Breeding programs are essential tools that contribute to the long-term development of aquaculture in order to improve the quality of the farmed strains. Aquaculture production based on genetic selection schemes has increased from 1% in 1997 to 80-83% in 2017 (Gjedrem, 1997(Gjedrem, , 2004(Gjedrem, , 2012Janssen et al., 2017). In gilthead seabream, genetically selected seed production only accounts 31-44% of total fry in the European market (Chavanne et al., 2016), with an estimated increase in growth of 5-29% (Knibb, 2000;Brown, 2003;Thorland et al., 2015). However, in this species quality of morphology and growth traits are deemed to have the same level of importance in the industry (Janssen et al., 2017) as the fish are sold whole and due to the loss of market value throughout the value chain (Bardon et al., 2009;Boglione and Costa, 2011;Boglione et al., 2013).
Skeletal deformity is the most important trait when defining the morphological quality of the fish, reducing the physiological ability of fish for correct development i.e., reducing their growth rate, increasing their mortality rate and significantly affecting their welfare (Andrades et al., 1996;Karahan et al., 2013;García-Celdrán et al., 2016). The presence of deformities can affect up to 30% of production, causing an annual loss of more than €50M in the European aquaculture industry (Castro et al., 2008;Fernández et al., 2008;Haga et al., 2011). In gilthead seabream, deformity appears early in development and its prevalence increases with age (Lee- Montero et al., 2015).
Abiotic, biotic, xenobiotic, nutritional and environmental factors have been studied as the causes of deformities (Afonso and Roo, 2007), as well as genetic factors or interactions (Andrades et al., 1996;Afonso et al., 2000;Fernández et al., 2008). Although environment has a large effect on deformity, different families or cohorts have shown a variable response to specific environmental stressors, suggesting additive genetic variation in these traits (Kause et al., 2007;Lee-Montero et al., 2015;García-Celdrán et al., 2016), the presence of QTLs (Negrín-Báez et al., 2015a,b) or inbreeding depression as a consequence of the consanguinity level (Aulstad and Kittelsen, 1971;Kincaid, 1983). In gilthead seabream, segregation of a major gene (Astorga et al., 2004) and polygenic inheritance (Afonso et al., 2000;Lee-Montero et al., 2015) have been suggested as explanations for the incidence of skeletal deformities. Negrín-Báez et al. (2015b) reported a significant excess of descendants with severe deformities from directed crosses involving parents with the same deformities, suggesting the elimination of deformed fish from a breeding nucleus or the inclusion of this trait in breeding programs.
In sparid species, evidence of mineralization can be discerned from 30 DPH (Socorro, 2006), where the temperature-dependent differential expression of genes associated with bone formation suggests the long-term regulation of osteogenesis (Riera-Heredia et al., 2018). In gilthead seabream, jaw and mouth bones are formed after the onset of exogenous feeding (Faustino and Power, 2005). In fact, additive components already appear from the larvae (Fragkoulis et al., 2018(Fragkoulis et al., , 2020 and fingerling stages (García-Celdrán et al., 2016).
Therefore, a major issue in the aquaculture industry is the production of a large number of eggs through controlling gamete quality in order to satisfy market demand for seeds through breeding programs (Bobe and Labbé, 2010;Migaud et al., 2013). Fish egg quality, measured through traits such as fertilization rate, oocyte yield and alive offspring number, is considered very important in accordance with European scientific criteria (AQUAEXCEL, 2013), including good survival and lack of abnormalities (Kjørsvik et al., 1990;Brooks et al., 1997). Gametes, both in terms of their quantity and quality, are closely associated with the competitiveness and sustainability of hatcheries and farms (Bromage et al., 1992;Migaud et al., 2013;Theodorou et al., 2016).
The aims of this study were to research, for the first time in gilthead seabream, the additive genetic component effect of skeletal deformities on spawning quality and to establish the main traits for evaluating reproductive ability at an industrial level. For this purpose, estimated breeding values for the presence or absence of deformity when the breeders were a commercial size (EBV def ) was calculated and their reproductive ability traced throughout an entire spawning season.

MATERIALS AND METHODS
All animal experiments were conducted in accordance with the European Union Directive (2010/63/EU) on the protection of animals for scientific purposes, at the Fundación Canaria Parque Científico Tecnológico (FCPCT) of the Universidad de Las Palmas de Gran Canaria (Canary Islands, Spain).

Biological Material and Experimental Design
A total of 4,108 gilthead seabream adults from the Canary Islands region and the 3 rd generation of the Spanish National Breeding Program (PROGENSA R ) (Lee-Montero et al., 2015) at AQUANARIA S.L. were visually asessed for the presence or absence of deformity at commercial size (ATOL:0000087), in accordance with AquaExcel-ATOL (AQUAEXCEL, 2013). The genotyping of all animals evaluated for morphology was also established through the molecular characterization of microsatellite markers (SMsa1), with PCR methodology, following the instructions of its authors (Lee-Montero et al., 2013). The genetic relationship between descendants (study population) and their parents (parental assignment) was reconstructed, comparing breeders' and descendants' genotypes by using VITASSIGN software (Vandeputte et al., 2006). Genetic parameters were estimated for both data matrix (phenotypic and genotypic information) by using the restricted maximum likelihood method with VCE 6.0 software (Neumaier and Groeneveld, 1998;Groeneveld et al., 2010). A stock of breeders was then preselected on the basis of its EBV def for conditioning and sexing before the spawning season, and moved to the FCPCT-ULPGC facilities. Fish with a deformed phenotype were discarded.
The EBV def of the 4,108 gilthead seabream stock ranges between −0.10515 for normality and +0.07821 for deformity (Figure 1), with an average value of −0.01551 and a standard deviation value of 0.02681. After discarding the dead and misshapen, a total of 32 breeders were selected on the basis of their EBV def , gender, weight and relationship coefficient. In this way, it was possible to establish two groups of breeders with opposite values in their EBV def , which corresponds to a value of +0.05443 in the group of genetically deformed (gD) breeders, and −0.05964 in the group of genetically normal (gN) breeders. Between EBV def values of gD and gN selected breeders, was contained almost 96% of the evaluated population (Figure 1).
Two different types of crosses, with different EBV def , were established in 2 m 3 cylindroconical tanks, with duplicates to minimize the environmental factors. The first was a cross of gN fish (N1 and N2) and the second a cross of gD fish (D1 and D2), with the same intensity of selection but reverse EBV def (EBV def = −0.05 vs. EBV def = +0.05, respectively). The crosses were established by taking into account the EBV def , relationship coefficient, biomass, and the weight ratio of males and females, so that the duplicates were similar and the differences between crosses were exclusively due to the EBV def . The biomass sex-ratio (kg male/kg female) was similar between tanks; 1.43, 1.3, 1.18, and 1.1 for N1, N2, D1, and D2, respectively, with a mean of 1.28 kg male/kg female, as recommended by Fernández-Palacios et al. (1995, 1997. All breeders had normal phenotypes and the same age, between and within crosses, aiming to avoid the impact of females age on spawning quality (Jerez et al., 2012), by keeping the biomass sex-ratio constant. Thus, breeders' maturation and spawning occurred spontaneously under natural conditions. The coefficient of relationship, EBV def , total biomass, male and female biomass, number of males and females were measured in breeders in the four experimental tanks (Figure 2). Tanks were supplied with 16 L/min seawater at 7.98 ± 0.10 pH and kept under natural light photoperiod (11-13 h light). Breeders were fed a commercial diet (Skretting ARC, Stavanger, Norway).
During the spawning season, the monthly average temperature decreased between December and February from 20.9 to 19.9 • C, followed by an increase to 22.8 • C in June. Food intake, expressed as a percentage of biomass, remained constant throughout the experiment with a range from 0.1 to 1.1% and an average value of 0.4% in all tanks.
FIGURE 1 | The distribution of 4,108 breeders from the 3 rd generation of the Spanish National Breeding Program (PROGENSA R ) according to their estimated breeding value for the presence or absence of deformity at commercial size (EBV def ). In yellow is the mean EBV def of the genetically normal (gN) preselected breeders and the percentage of breeders with a lower value. In blue is the mean EBV def of the genetically deformed (gD) preselected breeders and the percentage of breeders with a higher value. In red is the percentage of breeders with intermediate EBV def . FIGURE 2 | Experimental design representation that shows the two genetic groups (genetically deformed, gD; genetically normal, gN) with their duplicates (D 1 , D 2 and N 1 , N 2 ), the number of breeders per tank and their data: kinship (a), estimated breeding value for the presence or absence of deformity at commercial size (EBV def ), biomass, female weight (♀weight), and male weight (♂weight).
Samples of eggs were collected and observed under a binocular stereo microscope. The oocyte yield, fertilization rate and viability rate were estimated volumetrically after samples were counted. In addition, two 96-well plates were used to hold 192 floating eggs (1 per well) with 200 µl filtered seawater and incubated at 19 • C. After 24 and 96 h, the number of alive larvae was counted in order to calculate the hatching rate (number of hatched larvae/number of fertilized eggs) and the survival larval rate (%).
The number of fertilized eggs, viable eggs, hatched eggs and the number of alive larvae were also calculated as associated traits, since they represent the ultimate aim of the industry, a collection of quantitative and qualitative data that determines the final number of larvae available for rearing.

Statistical Analysis
The data were tested for normality using the one sample Kolmogorov-Smirnoff test, as well as for homogeneity of variance using Levene's test. When a normal distribution and/or homogeneity of variance was not achieved, data were subjected to the Kruskal-Wallis non-parametric test (Zar, 1984). To study the variation of spawning quality traits throughout the whole spawning season, comparisons between types of crosses (gD vs. gN) and between experimental tanks (D1, D2, N1, N2) were carried out, using the day as a unit of measurement ( Table 1). The variation of spawning quality traits was also studied in short time periods, structuring the whole spawning season into twelve fortnights. Each spawning quality trait was compared between types of crosses (gD vs. gN) within each fortnight. Additionally, each trait was compared within types of crosses (gD or gN) between fortnights ( Table 2). The following general linear model was used: Where, µ is the mean of the population, α i is the effect of the genetic factor (types of crosses), β ij is the effect of the tank factor within genetic factor, and ε ijk is the residual error.
For the oocyte yield and number of alive larvae traits, within fortnights, the influence of genetic and tank factors was measured by the following statistical algorithm (SA): execution of a standard regression analysis, and estimation of variance components of principal factors. Within fortnights, the same statistical procedure SA was also applied to explain the number of alive larvae trait using as predictors the oocyte yield, fertilization rate, viability rate, hatching rate and larvae survival rate traits. Statistical analysis was performed using the SPSS package, version 22.0 (SPSS Inc, Chicago, United States). All statistically significant comparisons of this work, Tables 1, 2, include the correction factor of Bonferroni's to minimize the Type-I error (Gordon et al., 2007; see Supplementary Material for differences between Bonferroni and Benjamini-Hochberg tests). Thus, significant differences were considered from 0.05, as confident level.

RESULTS
In all experimental tanks, breeders spawned daily between December 2018 and June 2019. The numbers of spawning days registered were 86, 116, 176, and 177 days for tanks D1, D2, N1, and N2, respectively. In Table 1, the mean values of all spawning quality traits are reported per tank and type of cross.
In Table 2, the comparisons for all spawning quality traits are presented: within fortnight per types of crosses, and within type of cross per fortnights.

Oocyte Yield
Oocyte yield contribution among tanks, for both types of crosses, was generally higher and more reproducible in the middle of the spawning season and lower at either ends of it (Figure 3). Across the whole spawning season the values between tanks were only significantly different within the gN group. Between genetic groups, the gN group reported a value 5.27 times higher than the gD group (Table 1). This significant difference was maintained across all fortnights (Table 2 and Figure 4A). Over the fortnights, within each genetic group both registered their best values in the first phase of the spawning season (fortnights two-five for gD, and fortnights one-seven for gN). Within fortnights oocyte yield was explained by genetic factor, from 29.1% in fortnight twelve to 90.9% in fortnight six ( Figure 5A) and a mean value across fortnights of 59.6%. The tank factor explained less in the middle (0% in fortnight six) and more at the beginning (44% in fortnight two) and at the end (65% in fortnight eleven) of the spawning season, with a mean value of 23.5%.

Fertilization, Viability, and Hatching Rates
In the whole spawning season these three traits did not report statistical differences between genetic groups, while statistical differences were estimated within genetic groups between tanks ( Table 1).
Within fortnights, statistical differences were detected between genetic groups in fortnight two for fertilization rate, in fortnights three, five, six, eight, and nine for viability rate, and in fortnights one, three, and nine for hatching rate ( Table 2 and Figures 4B-D). The fertilization rate of gD was lower than for gN across the fortnights. The viability rate of gD was higher than gN across fortnights at the beginning of the spawning season (fortnights one, two, and three) and lower than gN in the rest of the fortnights. The trend for hatching rate between gD and gN groups did not follow any regular pattern ( Table 2).
Between fortnights, within genetic groups, fertilization rate for both groups showed a mild increase toward the end of the spawning season. In the case of gN group, this rate decreased in the last two fortnights. For both genetic groups, the viability rate reported increases from the beginning of the spawning season, peaking in fortnight three and then decreasing until the end. The viability rate was more stable for the gN group. The hatching rate did not report any differences (Figures 4B-D).

Larval Survival Rate
For the whole spawning season, the gN group reported a value 7.2% higher than the gD group (Table 1). Within the same genetic group between tanks, only significant differences were found for the gD group.
Within fortnights, at the beginning of the spawning season (fortnights one-three), the gD group showed higher values than the gN group, a trend which was reversed for the rest of the spawning season (Table 2 and Figure 4E). Regarding the evolution of this trait across the fortnights, it was observed that the values of the gN crosses were generally constant throughout the spawning period, while in the gD crosses the values declined slightly at the end of the spawning period ( Figure 4E).

Fertilized Eggs, Viable Eggs, Hatched Eggs, and Number of Alive Larvae
In the whole spawning season these traits all differed significantly between genetic groups, the gN group being almost 6 times higher than gD group (Table 1). Within the same genetic group between tanks there were significant differences only in the gN group. Within fortnights the gN group was superior throughout the whole spawning season (Table 2 and Figure 4F). Between  fortnights, the gD and gN groups showed the same pattern in the distribution of values across all traits: low values at the beginning, high values in the middle and low values at the end ( Table 2).
For the number of alive larvae, the explanation of its variability within fortnights due to genetic and tank factors is shown in Figure 5B. Thus, genetic factors were highest in the middle of the spawning season (from fortnights five to eight), ranging from 63 to 87%. Conversely, the impact of the tank was high after the beginning of the spawning season, low in the middle and at its highest at the end (80-84%). However, environmental factors (not explained) were especially influential at the beginning (56-46% in fortnights one to two, respectively) and at the end (66% in fortnight twelve) of the spawning season. Figure 5C shows the importance of the oocyte yield, fertilization rate, viability rate, hatching rate, and larvae survival rate on the number of alive larvae. It can be seen that oocyte yield and viability rate explain the majority of variance. The influence of oocyte yield was at its maximum between fortnights three and seven (70.8-89.1%). Conversely, viability rate was at its maximum between fortnights nine and eleven (58-89.7%). Taking all traits into account, the best linear model explained between 84.5 and 98.7% of total variation of the number of alive larvae.

DISCUSSION
Body malformation is one of the most important traits defining the quality of batches. Deformities have a negative effect on the profit of aquaculture companies, especially hatcheries, because on-growing companies do not accept deformed fry (Afonso and Roo, 2007). Most skeletal anomalies in marine aquaculture species occur during their embryonic, larval and metamorphosis stages (Koumoundouros, 2010), and poor spawning quality is often associated with the appearance of such deformities. In gilthead seabream, variation in the inbreeding of breeders produces differences in spawning quality (Astorga, 2005).
In this study the effect of the EBV def on spawning quality traits (oocyte yield, fertilization rate, viability rate, hatching rate, larval survival rate, fertilized eggs, viable eggs, hatched eggs and number of alive larvae) is reported for first time in gilthead seabream. Breeders started spawning naturally at the beginning of the species' spawning season in December-January and continued until April-June (Arias, 1980;Fernández-Palacios et al., 1995;Kissil et al., 2001;Ibarra-Zatarain and Duncan, 2015). The number of spawning days from all tanks ranged between 86 and 116 days for the gD group, and 176-177 for the gN group. These results agree with values reported by Jerez et al. (2012) (121-130 days of spawning). The difference in spawning days observed between both genetic groups was robust because it was reproduced in the same way within the genetic groups between tanks, suggesting at least a genetic effect based on the EBV def .
The oocyte yield trait reported the greatest variability, showing differences between genetic groups (fed with commercial diets). In any event, their values ranged according to the locality and environmental conditions (Fernández-Palacios et al., 1995). While the oocyte yield of gN group was in line with   Xu et al. (2019), oocyle yield of gD group was 5.3 times lower in comparison to the gN group for the same time period. These differences were kept within fortnights. On the other hand, the gD group reported values lower than those reported by Jerez et al. (2012), and much lower than those reported by Xu et al. (2019), using diets with nutritional deficiencies. Oocyte yield trait was widely influenced by genetic factor throughout the whole spawning season, while environmental factors was only relevant at the beginning and end of the spawning season. In line with this, breeders with high genetic content in morphological deformity stimulate poor spawn quality through a lower level of oocyte yield and would be characterized for their EBV def content.
Traits as fertilization, viability, and hatching rates are widely used for controlling spawning quality but in this study, they did not show significant differences between genetic groups, although viability rate and hatching rate estimates were, respectively, lower and similar to other studies (Fernández-Palacios et al., 1995, 1997Scabini et al., 2011;Xu et al., 2019;Ferosekhan et al., 2020). The number of breeders and sex-ratio used in both types of crosses per tank (gN and gD groups) were in concordance with value used by Ibarra-Zatarain and Duncan (2015), in reproductive season of gilthead seabream. In fact, the fertilization rate registered was in line with values reported by other authors, also using mass spawning to study this species (Fernández-Palacios et al., 1995, 1997Scabini et al., 2011;Jerez et al., 2012;Ibarra-Zatarain and Duncan, 2015), and higher than other research using mating with only one male or female (Gorshkov et al., 1997;Xu et al., 2019;Ferosekhan et al., 2020).
Larval survival rate reported in this study was within range for this species (Carrillo et al., 1989;Fernández-Palacios et al., 1995, 1997Scabini et al., 2011;Xu et al., 2019;Ferosekhan et al., 2020) and the higher survival of the gN group compared to the gD group suggest a major susceptibility of the latter to environmental factors.
The large differences observed for fertilized eggs, viable eggs, hatched eggs and the number of alive larvae between genetic groups were similar along the whole spawning season, having significant differences between groups within fortnights (6 times higher in the gN group vs. the gD group). Moreover, between fortnights and independently of the genetic group, the same pattern in the distribution of values across these traits was reported: significantly low values at the beginning, high values in the middle, and low values at the end. Oocyte yield values described a distribution pattern practically equal but with differences of 5.3 times between groups, showing the high responsibility of it on all those traits, as expected.
The number of alive larvae is determinant for defining the effective production of hatcheries, in order to satisfy the on-growing companies' demand, which requires planning. It is therefore essential to know the spawning window within which the spawning quality is optimal. Genetic improvement is accumulative, permanent, and extendable to the whole production chain (Falconer and Mackay, 2001). In this study, the genetic factor showed a significant impact on the number of alive larvae (63-87%), mainly in the middle of the spawning season, in line with the period when spawning is considered optimal at an industrial level (Navarro et al., 2009).
A linear regression modeling of spawning for number of alive larvae, revealed that oocyte yield and viability rate traits explained the majority of its variance, being maximum in the middle of the spawning season by the oocyte yield. These results are in concordance with the empiric suggestions reported by Fernández-Palacios et al. (1995).
In conclusion, this study demonstrates how the deformity genetic background of phenotypically normal breeders affects and explains their spawning quality, for the first time in gilthead seabream. Shorter spawning periods, lower oocyte yield and, consequently, a lower number of alive larvae were found among gD breeders. The number of alive larvae, the most determinant trait of spawning quality, peaked in the middle of the spawning season on the basis of oocyte yield and genetic factor, reflecting the importance of controlling the breeders' genetics and the utility of oocyte yield trait for tracking spawning quality. Conversely, the viability rate and tank factor supported the number of alive larvae at the beginning and end of the spawning season. Given the results obtained, for oocyte yield trait and number of alive larvae, the gN breeders showed values 5.3 and 6 times higher than gD breeders, respectively. For future studies, it would be useful to estimate the differential gene expression patterns between breeders with different EBV def values, as well as between their offspring, in order to discover the genetic structure involved in skeletal deformity determination, in gilthead seabream.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

ETHICS STATEMENT
Ethical review and approval was not required for the animal study because all the animal experiments were performed according to the European Union Directive (2010/63/EU) on the protection of animals for scientific purposes, at Fundación Canaria Parque Científico Tecnológico (FCPCT) of the Universidad de Las Palmas de Gran Canaria (Canary Islands, Spain).

ACKNOWLEDGMENTS
We would like to thank the technical staff of the IU-ECOAQUA, where research was conducted. We would also like to thank AQUANARIA S.L. for the maintenance of breeders before preselection and Hipólito Fernández-Palacios for his recommendations concerning breeder management.