Heritability Estimates and Genetic Correlation for Growth Traits and LCDV Susceptibility in Gilthead Sea Bream ( Sparus aurata )

: The lymphocystis disease (LCD) is a viral infection with a high economic impact in gilthead sea bream aquaculture. In this study, genetic estimates associated with lymphocystis disease virus (LCDV) disease susceptibility and growth were determined in sea bream juveniles. Two ﬁsh batches (named batch 1 and batch 2) were built from mass spawning and reared under industrial conditions until disease outbreak. At the moment of the sampling ( n = 500 specimens for each batch), all animals had the typical LCDV lesions in the skin. For phenotyping, animals were weighted and photographed for image analysis (surface covered and lesion intensity). LCDV DNA copies were quantiﬁed in the liver by qPCR. Batch 1 had a higher surface covered and lesion intensity than batch 2, and the body caudal region was the lowest a ﬀ ected region in both batches. The average LCDV DNA copies in liver were higher in the batch 1 than batch 2, and they were positively correlated with severity index (SI) categories ( r 2 = 0.90–0.94). The total number of families evaluated were 150 and 128 for batch 1 and batch 2, respectively, with a high bias in o ﬀ spring contribution by family and broodstock. Heritabilities for weight and length were 0.18 and 0.14 in batch 1 and 0.06 and 0.05 in batch 2, respectively. Heritability for the number of viral DNA copies was low ( < 0.08) in both batches. Heritabilities for SI in binary scale were 0.32 / 0.33 and 0.21 / 0.24 (underlying liability / Bayesian approach) for batch 1 and batch 2, respectively. Genetic correlations were very high and positive when growth traits (weight and length) or disease traits (LCDV DNA copies and SI) were compared. batch results indicate the genetic selection for LCDV susceptibility and growth is feasible in sea bream juveniles, although estimates are highly dependent on the age. The information provided is relevant to designing selective breeding programs in sea bream.


Introduction
The lymphocystis disease (LCD) is a highly contagious viral infection responsible for high economic losses in the aquaculture industry worldwide [1]. The causative agent is the lymphocystis disease virus (LCDV) that belongs to the genus Lymphocystivirus, family Iridoviridae. This virus is hosted in the infected fish for long periods causing no clinical Iridoviridae signs (known as asymptomatic carriers) until host immunocompetence decreases due to different factors such as temperature shifts, grading operations, transport or high tank densities. Then, the virus actively multiplies in host cells giving rise to the typical lesions (named lymphocysts) consisting of irregularly shaped, randomly distributed white masses in the integumentary surfaces that in severe infections aggregate in nodular structures [2]. Experimental challenges have demonstrated that at least three weeks after LCDV injection are required for the observation of clinical signs [1,3,4]. Once the lesions appear, the disease spans for approximately 20 days, the time required for the hypertrophic dermal cells to mature and break resulting in high morbidity but low mortality rates unless secondary infections occur [2,3,5,6]. In spite of being a self-limiting benign process [3], LCDV infection has a high economic impact in hatcheries due to the important delay in fish growth, reducing feed conversion rates and making fish unavailable for commercialization due to external lesions and the lack of acceptance by the consumers.
Several approaches have been established to control LCDV outbreaks in the hatcheries in the absence of efficient commercial vaccines. The most widespread approach is the monitoring and removal of asymptomatic carriers by using highly sensitive diagnostic methods [7,8]. This strategy has been successfully applied to prevent vertical transmission, but it requires a continuous evaluation of broodstocks and does not prevent horizontal transmission. A second approach is the use immunostimulants supplied through the diet to enhance immune response and block virus replication [9,10]. This strategy is especially useful when stressful operations in the hatcheries are planned and preventive measures can be implemented. A third approach is genetic selective breeding for resistance lineages. In Paralichthys olivaceous, a major locus associated with resistance against LCDV genotype II was identified [11] and used for the design of marker-assisted breeding programs and LCDV-resistant broodstocks [12]. This latter strategy has been revealed as one the most sustainable strategies to reduce and minimize the impact of LCDV outbreaks in aquaculture, although it requires evaluation of genetic variation in species-specific populations.
The gilthead sea bream (Sparus aurata) is one the most cultured fish species in the Mediterranean basin. This species is infected by genotype VII, which is frequently observed in other Mediterranean species such as Senegalese sole [13,14]. LCDV outbreaks are recurrently observed every year and, occasionally, associated with high mortalities [3,15]. Currently, most ongoing genetic programs for this species are mostly focused on the improvement of growth, morphology and carcass quality [16][17][18][19]. Genetic estimates for disease resistance were only reported for the bacterial pathogen Photobacterium damselae subsp. piscicida [20,21]. However, the genetic variation controlling susceptibility to LCDV in sea bream remains still unexplored. This study aimed to estimate genetic parameters (heritabilities and genetic correlations) for LCDV susceptibility (as evaluated by the number of DNA viral copies and the severity of lesions as function of surface covered and intensity) and growth (weight and length) in gilthead sea bream. The results are the first step to create new genetically-improved broodstocks with a lower susceptibility to this disease that can enhance the sustainable production of this important aquacultured species.

Biometric Data
Two fish experimental batches (batch 1 and batch 2) were evaluated for LCDV susceptibility after a disease outbreak in a commercial hatchery. These fish batches were built by collecting the eggs during four consecutive days (4DL model) for family production from three sea bream broodstocks, and they were always managed as a single unit under the same culture conditions from larval stages Fishes 2020, 5, 2 3 of 12 until disease outbreak and sampling. At the moment of disease outbreak, average weight was 1.3 ± 0.5 g for batch 1 (~80 days post-hatch-dph) and 3.8 ± 1.1 g for batch 2 (~140 dph) (p < 0.05; Figure 1a). Average length was 3.8 ± 0.4 cm for batch 1 and 5.3 ± 0.4 cm for batch 2 (p < 0.05; Figure 1b).
(a) (b) Figure 1. Weight (a) and length (b) for batch 1 and batch 2. The mean ± standard deviation is indicated (total n = 498 for each batch). Asterisk denotes significant differences between batches.

Disease Phenotyping
All sampled animals had external skin lesions compatible with the LCDV infection (100% morbidity). To quantify the percentage of surface covered by lesions and their intensity, each animal was photographed and later analyzed by image analysis for three regions (R1-R3) (see Section 4.2 and Figure S1). Phenotype data for surface and intensity traits are depicted in Figure 2. Wilcoxon test showed that surface covered in batch 2 was significantly lower than batch 1 (p < 0.05; average 2.55 and 2.49 for batch 1 and batch 2, respectively). Analysis of the body region revealed that the surface of regions 1 and 2 was more affected (2.73 ± 0.47 and 2.66 ± 0.55) than region 3 (2.27 ± 0.77) (p < 0.05) in batch 1 ( Figure 2). In contrast, a gradient from region 1 (2.86 ± 0.35) to region 3 (1.92 ± 0.83) (p < 0.05) was detected in batch 2. Lesion intensities were also significantly (p < 0.05) higher in batch 1 than batch 2 (average 1.9 ± 0.7 vs. 1.5 ± 0.6) ( Figure 2). No differences in intensity among regions were observed (1.79-1.93 for batch 1 and 1.4-1.6 for batch 2).  The mean ± standard deviation is indicated (total n = 498 for each batch). Asterisk denotes significant differences between batches.

Disease Phenotyping
All sampled animals had external skin lesions compatible with the LCDV infection (100% morbidity). To quantify the percentage of surface covered by lesions and their intensity, each animal was photographed and later analyzed by image analysis for three regions (R1-R3) (see Section 4.2 and Figure S1). Phenotype data for surface and intensity traits are depicted in Figure 2. Wilcoxon test showed that surface covered in batch 2 was significantly lower than batch 1 (p < 0.05; average 2.55 and 2.49 for batch 1 and batch 2, respectively). Analysis of the body region revealed that the surface of regions 1 and 2 was more affected (2.73 ± 0.47 and 2.66 ± 0.55) than region 3 (2.27 ± 0.77) (p < 0.05) in batch 1 ( Figure 2). In contrast, a gradient from region 1 (2.86 ± 0.35) to region 3 (1.92 ± 0.83) (p < 0.05) was detected in batch 2. Lesion intensities were also significantly (p < 0.05) higher in batch 1 than batch 2 (average 1.9 ± 0.7 vs. 1.5 ± 0.6) ( Figure 2). No differences in intensity among regions were observed (1.79-1.93 for batch 1 and 1.4-1.6 for batch 2).
(a) (b) Figure 1. Weight (a) and length (b) for batch 1 and batch 2. The mean ± standard deviation is indicated (total n = 498 for each batch). Asterisk denotes significant differences between batches.

Disease Phenotyping
All sampled animals had external skin lesions compatible with the LCDV infection (100% morbidity). To quantify the percentage of surface covered by lesions and their intensity, each animal was photographed and later analyzed by image analysis for three regions (R1-R3) (see Section 4.2 and Figure S1). Phenotype data for surface and intensity traits are depicted in Figure 2. Wilcoxon test showed that surface covered in batch 2 was significantly lower than batch 1 (p < 0.05; average 2.55 and 2.49 for batch 1 and batch 2, respectively). Analysis of the body region revealed that the surface of regions 1 and 2 was more affected (2.73 ± 0.47 and 2.66 ± 0.55) than region 3 (2.27 ± 0.77) (p < 0.05) in batch 1 ( Figure 2). In contrast, a gradient from region 1 (2.86 ± 0.35) to region 3 (1.92 ± 0.83) (p < 0.05) was detected in batch 2. Lesion intensities were also significantly (p < 0.05) higher in batch 1 than batch 2 (average 1.9 ± 0.7 vs. 1.5 ± 0.6) ( Figure 2). No differences in intensity among regions were observed (1.79-1.93 for batch 1 and 1.4-1.6 for batch 2).  To better characterize the genetic variation and the correlations between productive traits, the severity index (SI) was used to classify the animals into four categories: Weak, moderate, severe or very severe ( Figure 3). The average SI was 5.0 ± 2.6 and 3.9 ± 2.2 for batch 1 and batch 2, respectively. In batch 1, 60% of fish were considered as severe or very severe infected and this percentage dropped to 40% in batch 2 ( Figure 3).
To better characterize the genetic variation and the correlations between productive traits, the severity index (SI) was used to classify the animals into four categories: Weak, moderate, severe or very severe ( Figure 3). The average SI was 5.0 ± 2.6 and 3.9 ± 2.2 for batch 1 and batch 2, respectively. In batch 1, 60% of fish were considered as severe or very severe infected and this percentage dropped to 40% in batch 2 ( Figure 3).

Quantification of Viral DNA Copies
Viral DNA copies were quantified in the liver of fish from batch 1 and batch 2 ( Figure 4a). The average number of LCDV copies was 9.2 × 10 7 ± 7.2 × 10 8 µ g total DNA −1 in batch 1 and 6.3 × 10 5 ± 5.9 × 10 6 µ g total DNA −1 in batch 2. Interestingly, a high and positive correlation between the average number of viral DNA particles in liver and the SI categories (r 2 = 0.90-0.94) was determined ( Figure  4b).

Quantification of Viral DNA Copies
Viral DNA copies were quantified in the liver of fish from batch 1 and batch 2 ( Figure 4a). The average number of LCDV copies was 9.2 × 10 7 ± 7.2 × 10 8 µg total DNA −1 in batch 1 and 6.3 × 10 5 ± 5.9 × 10 6 µg total DNA −1 in batch 2. Interestingly, a high and positive correlation between the average number of viral DNA particles in liver and the SI categories (r 2 = 0.90-0.94) was determined ( Figure 4b).
To better characterize the genetic variation and the correlations between productive traits, the severity index (SI) was used to classify the animals into four categories: Weak, moderate, severe or very severe ( Figure 3). The average SI was 5.0 ± 2.6 and 3.9 ± 2.2 for batch 1 and batch 2, respectively. In batch 1, 60% of fish were considered as severe or very severe infected and this percentage dropped to 40% in batch 2 ( Figure 3).

Quantification of Viral DNA Copies
Viral DNA copies were quantified in the liver of fish from batch 1 and batch 2 ( Figure 4a). The average number of LCDV copies was 9.2 × 10 7 ± 7.2 × 10 8 µg total DNA −1 in batch 1 and 6.3 × 10 5 ± 5.9 × 10 6 µg total DNA −1 in batch 2. Interestingly, a high and positive correlation between the average number of viral DNA particles in liver and the SI categories (r 2 = 0.90-0.94) was determined ( Figure  4b).

Parental Contribution and Genetic Estimates
The 95.9% of the progeny in batch 1 were assigned to one breeder and 86% to a unique parent pair. These figures were 93.4% and 84.0% in batch 2. The total number of families were 150 and 128 for batch 1 and batch 2, respectively with the 87% offspring assigned F2 broodstock, 9.7% to broodstock F1 + F2 and 3.3% to the non-selected control broodstock. A total of 54 out of the 82 breeders contributed offspring although with a high bias in contribution percentages ( Figure S2). The reconstructed breeder sex ratio contributing offspring was 1M:2.8F and 1M:1.6F for batch 1 and batch Fishes 2020, 5, 2 5 of 12 2, respectively. Males showed a higher bias (p < 0.05) in offspring contribution per breeder than in females with maximal percentages ranging between 31% and 51% for batch 1 and batch 2, respectively (9% and 16%, respectively, in females). The analysis of SI and the number of LCDV DNA copies by family indicated a large variation that was even more evident in batch 2 ( Figure 5). The SI ranged from 1.7 to 9 in both batches and the log number of viral DNA copies between 4.6-7.9 (in logarithmic scale) in batch 1 and 1.8-6.2 in the batch 2.
The 95.9% of the progeny in batch 1 were assigned to one breeder and 86% to a unique parent pair. These figures were 93.4% and 84.0% in batch 2. The total number of families were 150 and 128 for batch 1 and batch 2, respectively with the 87% offspring assigned F2 broodstock, 9.7% to broodstock F1 + F2 and 3.3% to the non-selected control broodstock. A total of 54 out of the 82 breeders contributed offspring although with a high bias in contribution percentages ( Figure S2). The reconstructed breeder sex ratio contributing offspring was 1M:2.8F and 1M:1.6F for batch 1 and batch 2, respectively. Males showed a higher bias (p < 0.05) in offspring contribution per breeder than in females with maximal percentages ranging between 31% and 51% for batch 1 and batch 2, respectively (9% and 16%, respectively, in females). The analysis of SI and the number of LCDV DNA copies by family indicated a large variation that was even more evident in batch 2 ( Figure 5). The SI ranged from 1.7 to 9 in both batches and the log number of viral DNA copies between 4.6-7.9 (in logarithmic scale) in batch 1 and 1.8-6.2 in the batch 2. The heritability and genetic correlation estimates for weight, length, LCDV DNA copies and SI in binary scale (SI) are presented in Table 1. Heritabilities for weight and length were 0.18 and 0.14 in batch 1, respectively, and they decreased to 0.06 and 0.05 in batch 2. Heritability for the number of viral DNA copies was low in both batches (0.04 and 0.08 for batch 1 and batch 2, respectively).   The heritability and genetic correlation estimates for weight, length, LCDV DNA copies and SI in binary scale (SI) are presented in Table 1. Heritabilities for weight and length were 0.18 and 0.14 in batch 1, respectively, and they decreased to 0.06 and 0.05 in batch 2. Heritability for the number of viral DNA copies was low in both batches (0.04 and 0.08 for batch 1 and batch 2, respectively).
To estimate the genetic susceptibility to LCDV infection under the threshold model, animals were binary coded according to SI: Weak/moderate vs. severe/very severe. In the observed scale, the heritability for the binary SI was 0.20 ± 0.10 and 0.14 ± 0.08 for batch 1 and batch 2, respectively. When heritabilities were transformed to the liability scale, the values increased to 0.32 and 0.21, respectively. A Bayesian linear mixed model was also carried out to estimate heritability using the binary-coded LCDV susceptibility data. The Bayesian narrow-sense heritability estimate was 0.33 and 0.24 for batch 1 and batch 2, respectively, with a 95% Bayesian credibility interval for the additive genetic component of 0.12-0.67 and 0.01-0.53 ( Figure 6).
The genetic and phenotypic correlations between the weight, length, LCDV DNA copies and SI are given in Table 1. Genetic correlations were very high and positive between growth traits (weight and length) and between disease traits (number of LCDV DNA copies and SI) in both batches. However, when growth and disease traits were compared, the genetic correlations were positive and moderate-high in the batch 1 but negative in batch 2. It should be noted that for length and SI, the model converged but did not provide errors. Phenotypic correlations between weight and length were Fishes 2020, 5, 2 6 of 12 high and positive and moderate and positive between LCDV DNA particles and SI. In batch 1, the phenotypic correlations between biometric and susceptibility traits were close to zero while in batch 2 they were low and negative. Table 1. Heritabilities (in bold) and the phenotypic (below the diagonal) and genetic correlations (above the diagonal). Means ± standard deviations (SD) are indicated for body weight, fork length, LCDV DNA copies and SI in binary scale (SI).  To estimate the genetic susceptibility to LCDV infection under the threshold model, animals were binary coded according to SI: Weak/moderate vs. severe/very severe. In the observed scale, the heritability for the binary SI was 0.20 ± 0.10 and 0.14 ± 0.08 for batch 1 and batch 2, respectively. When heritabilities were transformed to the liability scale, the values increased to 0.32 and 0.21, respectively. A Bayesian linear mixed model was also carried out to estimate heritability using the binary-coded LCDV susceptibility data. The Bayesian narrow-sense heritability estimate was 0.33 and 0.24 for batch 1 and batch 2, respectively, with a 95% Bayesian credibility interval for the additive genetic component of 0.12-0.67 and 0.01-0.53 ( Figure 6). The genetic and phenotypic correlations between the weight, length, LCDV DNA copies and SI are given in Table 1. Genetic correlations were very high and positive between growth traits (weight and length) and between disease traits (number of LCDV DNA copies and SI) in both batches. However, when growth and disease traits were compared, the genetic correlations were positive and moderate-high in the batch 1 but negative in batch 2. It should be noted that for length and SI, the model converged but did not provide errors. Phenotypic correlations between weight and length

Discussion
Selective breeding for disease resistance is considered as a sustainable approach to support aquaculture growth. Previous studies in fish demonstrated that heritabilities for disease resistance traits oscillate from moderate to high, and hence they can be incorporated to breeding programs to minimize fish losses and enhance fish welfare [22]. In sea bream, genetic selection programs have focused mainly on growth, morphology and carcass quality traits [16][17][18]23]. To date, only heritabilities for fish survival after challenging with the pathogenic bacteria P. damselae subsp. piscicida were estimated, ranging from 0.12-0.28 [20,21]. In the case of LCDV, previous studies in P. olivaceous demonstrated that marker-assisted selection programs were successful to control the LCDV-C (genotype II) disease outbreaks [12]. In spite of the high incidence of LCDV in sea bream hatcheries, no report about the genetic variance for resistance is still available in part due to the difficulty to experimentally reproduce the disease, highly dependent on environmental conditions and fish immunological status. In the present study, two fish batches generated by mass spawning from three broodstocks were produced in a commercial hatchery in which LCDV outbreaks were previously reported. In both fish batches, the same environmental conditions under standard production procedures were carried out except for grading operations. A cross-sectional model was followed and fish were sampled for phenotyping and genotyping just when typical LCDV skin lesions were evident, assuming that LCDV infections occurred at least three weeks earlier [3].
As LCDV outbreaks are benign infections with low mortality rates [2,3,5,6], the measurement of traits associated to disease susceptibility should be carried out by quantifying the degree of clinical signs, such as body surface covered by lesions and their intensity. In our study, no mortality was observed at the sampling time, although all animals had the typical LCDV lesions ( Figure S1). To quantify the infection degree, an image analysis was carried out identifying a high variability in the body distribution and intensity of the lesions. Moreover, the numbers of LCDV DNA copies were calculated to determine the correlation with the severity of clinical signs. The qPCR analysis confirmed that all animals were positive to viral MCP protein [8], as expected by the lesions, although with a huge range of variation in viral loads. Moreover, the high correlation (>90%) between LCDV DNA copies and the four SI categories confirmed this index is a good indicator of viral replication. It should be noted that although the liver is not the target organ for LCDV viral replication, this virus is able to replicate on it before spreading in integumentary system [3,14,24]. The selection of this organ appears also as a key point to precisely quantify DNA copies to prevent cross-contaminations between skin lesions in highly infected fish populations.
The reconstructed familial structure revealed a high variance in the offspring contribution, a feature previously reported for mass-spawning in photoperiod-controlled sea bream broodstocks [20,25]. In the present study, only 59% of breeders contributed progeny with a higher contribution of female breeders that fit to the broodstock structure established at the beginning of the breeding period. The most striking was the strong bias for offspring contribution between F2 and the control broodstocks in spite of containing a similar number of breeders. Although the weight of breeders (average 1875 ± 240 vs. 731 ± 144 for F2 and control, respectively) could be a major factor behind the differences in egg production, the quadratic relationship reported for parental weights and progeny production [25] does not fully support this hypothesis. We cannot exclude that selective breeding for growth in F2 could also indirectly improved relative fecundity, a hypothesis that needs to be tested.
Heritability estimates for the hepatic number of LCDV DNA copies were almost negligible in both batches (0.04 and 0.08), but they showed a high genetic correlation with SI (>0.99). Nevertheless, the heritability in the liability scale for binary SI was estimated to be 0.32-0.33 and 0.22-0.24 for batch 1 and batch 2, respectively, indicating that this index could be a useful trait to evaluate the genetic susceptibility to LCDV. The heritabilities were moderate but similar to those reported for nodavirus resistance on the liability scale in the European sea bass (0.26 ± 0.11) [26] or against P. damselae subsp. piscicida in sea bream (0.12-0.28) [20,21], supporting the idea that selective breeding for LCDV susceptibility is feasible.
Previous reports indicated that heritabilities for furcal length in juvenile sea bream specimens are moderate (h 2 = 0.31 ± 0.07 for fish of 6.6 cm in length and age 130 dph [23] and h 2 = 0.38 ± 0.08 for fish of~2.3 cm in length and age 110 dph [20]. For weight, Navarro, Zamorano, Hildebrandt, Ginés, Aguilera and Afonso [23] showed that heritability at age 130 dph was 0.28 ± 0.07. Although the fish weight and length values in the present study are closer to those reported by Navarro et al. (2009), our heritability estimates are still a bit lower for both traits (0.18 and 0.14 for batch 1 for weight and length, respectively, and very close to zero in batch 2). These lower genetic estimates could be due to the energy reallocation mediated by LCDV from somatic growth to virus replication and fibroblast hyperplasia in the integumentary tissues. Moreover, the heritability differences between batches with different age, weight and SI suggest that bigger fish can mobilize energy towards the immune system response to fight against the LCDV infection masking the genetic effects on growth. Supporting this hypothesis are the genetic correlations between growth traits (weight and length) and the SI in both batches: Positive and high in batch 1 but negative in batch 2. One major constraint in this study is the lack of longitudinal information to span whole infection period, but the sampling of two aged batches with contrary genetic correlations indicates that susceptibility is highly dependent on fish age. Resistance to viral diseases is generally negatively correlated with weight [26,27]. In sea bream, the resistance to P. damselae subsp. piscicida was positively correlated to fork length, although the authors could not explain if the association was due to the survival time or the weight itself [20]. In our study, a high variation in weight was found within each batch, and the differences in genetic correlation seem to be more associated to development and a putative enhancement of immunocompetence than size itself, which allows fish to better deal with this benign infection. According to this hypothesis a mixed selection index that combines length and SI could be useful for selective breeding of superior animals, preserving growth potential but limiting the viral replication and phenotypic lesions.

Animals
Three broodstocks hosted at center IFAPA El Toruño (El Puerto de Santa María, Spain) belonging to the Spanish genetic breeding program PROGENSA ® were used for family production: (a) Non-selected animals or control broodstock (52 breeders, 4 years old); (b) F1 + F2 broodstock that included animals of two 1st (females) and 2nd (males) generations selected for growth (30 breeders; 7 years old in the case of females and 4 years old for males); (c) F2 broodstock that included only animals selected for growth of the 2nd generation (57 breeders; 4 years old). The female/male ratio was approximately 2:1 in the tanks. The three broodstocks were under a controlled photoperiod (8L:16D) to synchronize maturation and egg release that was initiated at the beginning of December 2016. During this period, animals were fed ad libitum Vitalis Cal (Skretting), and egg production was monitored daily. When total egg production became stable, two egg batches were established: One at the end of February and another in early April 2017. In both cases, eggs from the three broodstocks were collected and pooled for 4 consecutive days to maximize family production according to the 4DL model (Elalfy et al., unpublished). Incubation was carried out in cylinder conical tanks (1000 L) at a density of 500-1000 larvae·L −1 . Water conditions were as follows: Temperature 19.0 • C, salinity 34% and dissolved oxygen 6.4 mg·L −1 . All spawns were checked by PCR to confirm the absence of LCDV infections as previously reported [8]. As the LCD is very difficult to experimentally reproduce in the lab due to the highly dependence on environmental factors, viable non-hatched eggs and hatched larvae were transported to an industrial hatchery in Huelva (Spain) in which LCDV outbreaks were previously observed. Both larval batches were reared at an initial density of 100 larvae·L −1 in 5 m 3 tanks. Rearing protocols were those used in the hatchery at commercial scale, and these two batches were always managed as experimental units not graded until sampling.
The LCDV outbreak occurred in June 2017 with temperature 26 • C and salinity 36 ppt. At the moment of the sampling, clear lesions compatible with LCDV infections were evident in the animals. As our animals were cultivated under industrial conditions and all tanks were managed in the same way, no negative tanks were available. LCDV infection was confirmed by qPCR, and the degree of lesion intensities were determined image analysis. A total of 500 animals were randomly sampled from each batch and in situ sacrificed in slurring ice to preserve external lesions and thereafter transported to center IFAPA "El Toruño" for phenotyping and samplings. The two batches were named batch 1 (eggs from April; age~80 days post-hatch-dph) with a high number of typical lymphocysts and batch 2 (eggs from February;~140 dph) with less prominent lesions but all characteristic of LCDV infection.

Data Collection and Image Analysis
All animals were handled under strict biosecurity measures to avoid virus propagation in the facilities. All animals were photographed individually using a Canon EOs1300D camera following Fishes 2020, 5, 2 9 of 12 the methodology previously established in PROGENSA ® [19]. Then, animals were weighted and individually frozen at −20 • C for later tissue sampling. All procedures were authorized by the Bioethics and Animal Welfare Committee of IFAPA.
Image analysis was carried out using the Fiji 2.0.0-rc-69/1.52i [28]. All fish specimens were measured for furcal length and later evaluated for the extent of surface covered by lesions and their intensity by dividing the specimens in three corporal regions ( Figure S1): Region 1 (R1) that spanned from the caudal part of the operculum to the end of the pectoral fin; region 2 (R2) that spanned from the end of pectoral fin to the base of the caudal fin; region 3 (R3), caudal fin. Each region was ranked as 1 if lesions were covering less than 30% of the surface, 2 if between 30 and 70% and 3 if more than 70% ( Table 2). The head was not considered due to the difficulty of clearly visualizing the lesions. Moreover, lesion intensities were established for each region as 1 if they were not prominent, 2 if slightly prominent and 3 when lymphocysts were evident ( Table 2). All animals were independently evaluated by two researchers for phenotypic traits (surface and intensity), and both evaluations were used to build a consensus rating for each animal. Moreover, the severity index from surface (S) covered by lesions and their intensity (I) was calculated as follows: SI = S × I) using the mean of the three regions for S and I. Such an index oscillates between 1 and 9 and was used to define 4 categories for genetic analysis: Weak, moderate, severe or very severe ( Table 2).

Quantification of Viral DNA
Although the skin is the organ with the highest number of viral DNA copies, LCDV also replicates in the liver [3,24]. This organ was selected as the target organ to quantify the viral DNA copies in each specimen, since it provides enough biomass for the nucleic acid isolation minimizing the risk of cross-contamination between individuals. Livers were aseptically dissected, individually separated in eppendorf tubes and kept frozen at −20 • C until processing. All materials were cleaned with ethanol and distilled water between animals. In the case of batch 1, total DNA was isolated using Isolate II Genomic DNA Kit (Bioline) and in the case of batch 2, DNA samples were purified using phenol/chloroform/isoamyl alcohol (Sigma) following the manufacturer's instructions. DNA samples were treated with RNase A (Bioline) and quantified spectrophotometrically using the Nanodrop ND-8000. Absolute quantification of viral DNA copies was carried out according to the protocol specified by Valverde, Cano, Labella, Borrego and Castro [8] based on the quantification of the major capsid protein (mcp). Assays were run in a CFX96 TM Real-Time System (Bio-Rad) in a 10 µL final volume containing 200 ng of DNA, 300 nM each of specific forward and reverse primers and 10 µL of iQ™ SYBR ® Green Supermix (Bio-Rad). The amplification protocol used was as follows: Initial 7 min denaturation and enzyme activation at 95 • C, 40 cycles of 30 s at 95 • C and 1 min at 59 • C.

Genotyping and Parentage Assignment
DNA for offspring parentage assignment was isolated from liver as previously reported in Section 2.3. Breeders were sampled before the spawning season, and blood (~0.5 mL) was taken by puncturing in the caudal vein using a heparinized syringe, adding heparin (100 mU) and keeping at −20 • C until use. DNA from breeders was extracted from blood using the Isolate II Genomic DNA Kit and treated with RNase A (Bioline, London, UK) following the manufacturer's instructions.
Genotyping was carried out using multiplex PCR SMsa-1 and SMsa-2 [29]. PCRs were run on an ABI3130 Genetic Analyzer (Applied Biosystems, Foster City, CA, USA) and data analyzed using the Genemapper v3.8. Parentage assignment was determined using the exclusion method using Vitassign v8.2.1 [30]. The breeder gender was considered as unknown.

Statistical Analysis
Comparison of rates for surface covered by lesions and their intensity was carried out by nonparametric tests Wilcoxon test or Kruskal-Wallis. Data for weight and number of LCDV DNA copies were log-transformed to achieve adjustment to normality. Statistical analysis to test phenotype differences was carried out using Prism8 (GraphPad, San Diego, CA, USA). Variance components associated with the susceptibility to LCDV were estimated separately for batch 1 and batch 2 based on bivariate linear mixed models fitted by restricted maximum likelihood (REML) in Wombat (Meyer, 2007): y = Xβ + Zu + e, where y is the observed trait, β is the vector for fixed factor (broodstock origin), u is the vector for animal random factor and e is the error. Heritability and the phenotypic and genotypic correlations were calculated for four traits: Weight, length, number of LCDV DNA copies and LCDV susceptibility as defined by SI in binary scale: 0 = weak/moderate and 1 = severe/very severe. Some runs using the weight as covariate for SI were also carried out, but no effect on genetic estimates was determined and hence was not included in the analyses.
Heritability estimates for binary LCDV susceptibility in the observed scale were later transformed to the underlying liability scale according to Dempster and Lerner [31] using the formula h 2 = [p(1 − p)/z 2 ] × h 2 obs , where h 2 obs is the heritability on the observed scale (binary trait), h 2 is the heritability on the liability scale, p is the proportion of affected individuals and z is the density of a standard normal distribution at the pth quantile. In addition, heritability estimates for binary LCDV susceptibility were also calculated using a Bayesian framework under the threshold model as implemented in the MCMCglmm R package [32]. A χ 2 distribution with one degree of freedom was used as the prior distribution and we applied a probit scale and fixed the value of the residual variance to 1. MCMC ran for 1 million iterations with a thinning interval of 100 after a burn-in of 100,000 (de Villemereuil et al., 2013). Binary trait heritability was estimated as h 2 = Va/(Va + 1 + 1) as described in the course notes for the program MCMCglmm [32].

Conclusions
This study has optimized a methodology to phenotypically evaluate LCDV susceptibility and rank the animals according to SI, which is a mixed index combining surface covered by lesions and their intensity. Severity categories according to SI were highly correlated with hepatic viral DNA copies and demonstrated a moderate heritability, although they were highly dependent on the age, indicating that it could be a good trait to select animals with a lower LCDV susceptibility. The genetic correlations between growth and disease traits were moderate-high and positive in younger fish, indicating that a combined index should be established to successfully select stocks with low LCDV susceptibility and enhanced growth.