Genetic Variation, GWAS and Accuracy of Prediction for Host Resistance to Sparicotyle chrysophrii in Farmed Gilthead Sea Bream (Sparus aurata)

Gilthead sea bream (Sparus aurata) belongs to a group of teleost which has high importance in Mediterranean aquaculture industry. However, industrial production is increasingly compromised by an elevated outbreak of diseases in sea cages, especially a disease caused by monogeneans parasite Sparicotyle chrysophrii. This parasite mainly colonizes gill tissues of host and causes considerable economical losses with mortality and reduction in growth. The aim of current study was to explore the genetics of host resistance against S. chrysophrii and investigate the potential for genomic selection to possibly accelerate genetic progress. To achieve the desired goals, a test population derived from the breeding nucleus of Andromeda Group was produced. This experimental population was established by crossing of parents mated in partial factorial crosses of ∼8 × 8 using 58 sires and 62 dams. The progeny obtained from this mating design was challenged with S. chrysophrii using a controllable cohabitation infection model. At the end of the challenge, fish were recorded for parasite count, and all the recorded fish were tissue sampled for genotyping by sequencing using 2b-RAD methodology. The initial (before challenge test) and the final body weight (after challenge test) of the fish were also recorded. The results obtained through the analysis of phenotypic records (n = 615) and the genotypic data (n = 841, 724 offspring and 117 parents) revealed that the resistance against this parasite is lowly heritable (h² = 0.147 with pedigree and 0.137 with genomic information). We observed moderately favorable genetic correlation (R_g = −0.549 to −0.807) between production traits (i.e., body weight and specific growth rate) and parasite count, which signals a possibility of indirect selection. A locus at linkage group 17 was identified that surpassed chromosome-wide Bonferroni threshold which explained 22.68% of the total genetic variance, and might be playing role in producing genetic variation. The accuracy of prediction was improved by 8% with genomic information compared to pedigree.

Introduction

Gilthead sea bream (Sparus aurata) is an economically important farmed fish species, specifically to the Mediterranean aquaculture industry with an annual production of ∼ 160,563 metric tons (FEAP Secretariat, 2017). The enhanced production and the demand of gilthead sea bream is also accompanied by production challenges with an increase of disease outbreaks (viral, bacterial, and parasitic) in sea cages (Food and Agriculture Organization of the United Nations [FAO], 2014), especially those which are caused by monogeneans parasites (Sitjà-Bobadilla, 2004; Nowak, 2007) with direct life cycles on the host. Sparicotyle chrysophrii (S. chrysophrii), is a common monogeneans parasite which infests both wild and cultured gilthead seabream and instigates lethal epizootics in sea cages and is considered as one of the main diseases which threatens Mediterranean aquaculture. This parasite mainly colonizes gill tissues of the host and causes considerable economical losses with mortality ranging from 10 to 30% (Athanassopoulou et al., 2005; Muniesa et al., 2020), and, more importantly by reduction in growth of the farmed stocks triggered by emaciation and anemic state of the survivors (Sitjà-Bobadilla and Alvarez-Pellitero, 2009). Treatment of farmed fish against S. chrysophrii is very expensive, labor intensive, and requires resource demanding bath applications using anthelmintics which are also associated with environmental impacts (Basurco and Rogers, 2009; Rigos et al., 2013). Moreover, bath treatments may cause handling stress to the fish which may further increase losses.

Selection for genetic resistance against parasitic diseases is a highly valuable tool to help prevent or diminish disease outbreaks, especially where effective therapeutic agents or vaccines are limited or lacking (Bishop and Woolliams, 2014). The selection response is dependent on several factors including heritability, weight given to a trait, accuracy of selection and the generation interval. The pace of genetic progress can be increased through selective breeding (Ødegård et al., 2011) given that the trait is of moderate to high heritability (Odegård et al., 2014). Breeding programs in gilthead sea bream primarily included growth rate, and parasite/pathogen resistance with emphasis on minimizing productions costs (Chavanne et al., 2016; Janssen et al., 2017). Increasing the number of traits in the breeding goal could be complex and may impact selection response due to a reduced weight per trait and possible unfavorable genetic correlations among traits (Haffray et al., 2012; Yáñez et al., 2016).

The discovery of molecular markers [e.g., single nucleotide polymorphisms (SNPs)] and their associations with the traits of economic importance could prove very helpful to overcome unfavorable correlations among traits. Moreover, selection response in animal breeding programs especially for the traits which are difficult to improve (e.g., lowly heritable traits, meat quality traits, sex limited traits, etc.) by traditional selection can benefit from genome-wide distributed genetic markers and the knowledge about genetic markers linked to genes affecting quantitative traits (Meuwissen and Goddard, 1996; Pyasatian et al., 2007). The genetic markers have revealed significant associations with the quantitative traits of economic importance in Atlantic salmon (Everett and Seeb, 2014; Gonen et al., 2015; Moen et al., 2015).

Single nucleotide polymorphisms are popular molecular genetic markers widely used in aquaculture and livestock research (Gray et al., 2000; Aslam et al., 2012). Restriction-site associated DNA based sequencing (RAD-Seq) is a reduced representation high-throughput sequencing technique for the simultaneous detection and genotyping of individuals for the detected SNPs (Baird et al., 2008). The technique has been elaborated and used in several studies, briefly it employs (i) fragmentation of genomic DNA using restriction enzyme(s), (ii) selection of fragment size(s) from the digested DNA, unique individual specific nucleotide based barcoding and library preparation (iii) pooling the libraries from multiple samples, and high throughput sequencing. RAD-Seq technique has contributed in the development of genomic resources and studies on aquaculture species, i.e., discovery of SNPs, construction of high-density linkage maps (Kakioka et al., 2013; Palaiokostas et al., 2013a, b, 2015; Gonen et al., 2014) and genome wide association studies (GWAS) in a cost-efficient manner (Campbell et al., 2014; Palti et al., 2015; Tsai et al., 2016). A simplified version of RAD-Seq is known as 2b-RAD which requires type IIB restriction enzymes to digest genomic DNA from both ends (up and downstream) to the enzyme recognition site (Wang et al., 2012). The fragment size obtained with type IIB restriction enzymes are of uniform length (∼36 bp tags), and thus avoids the sampling error because size selection is not needed, which is required in RAD-Seq during the size selection step (Puritz et al., 2014). 2b-RAD approach has also been used and assessed by researchers for exploring genetic basis of traits in gilthead sea bream and other fish species (Fu et al., 2016; Palaiokostas et al., 2016; Pecoraro et al., 2016).

Moreover, there have been recent developments in genomic resources and tools for gilthead sea bream including the availability of a reference genome assembly (Pauletto et al., 2018; Pérez-Sánchez et al., 2019) and SNPs based genotyping array (Lozano, 2019). The development of these genomic resources/tools can also provide opportunities for advanced selection methods such as genomic selection. Genomic selection (GS) is an advanced method for predicted breeding values of candidates which utilizes the effects of genome-wide distributed markers potentially influencing quantitative traits with an objective to achieve reliable selection (Meuwissen et al., 2001). The use of GS becomes highly important for the traits which are difficult to measure on live selection candidates, e.g., carcass and disease resistance traits. The GS contribute through significant improvement in accuracy for predicting breeding values compared to accuracies obtained via traditional pedigree-based approaches (Palaiokostas et al., 2016; Tsai et al., 2016; Bangera et al., 2017; Vallejo et al., 2017).

The aim of current study was to explore genetic basis of host resistance to S. chrysophrii by estimating the role of genetic variation explaining resistance to this gill parasite, generate a linkage map, perform a GWAS to detect QTL, and finally compare accuracies of predictions based on genomic and pedigree data.

Materials and Methods

Resource Population

All the fish used originated from the breeding nucleus of the Andromeda Group which was developed with a cross of parents mated in ∼8 × 8 partial factorial crossbreeding design using 58 sires and 62 dams. The brooders of this species are known to be fluent spawners under natural photoperiod without any need of hormonal and/or light regimes. Parents were allowed to perform natural mass spawning events in eight different broodstock tanks. All the parents were tissue sampled for later parentage assignment and genomic analysis. The collection of fertilized eggs was performed in four consecutive days with equal number of eggs collected in each day from each of the eight tanks. The collected eggs were pooled each day to avoid any environmental variances. For this specific experiment, an average of around 250 fingerlings per tank (total N = 2,000) were randomly sampled from each tank and kept until they reached the average body weight of 16 g. The descendent fish taken from this mating design had approximately the same size.

Challenge Test

A challenge test was performed following the controllable cohabitation infection model developed by Rigos et al. (2015) using donor fish as a prime source of disease transfer. The infection model is detailed and discussed in Rigos et al. (2015). In brief, 2,000 healthy gilthead sea bream weighing 16 ± 4 g were transferred to the facilities of the Hellenic Centre for Marine Research (Athens, Greece). Prior to the transfer, the fish were maintained in land tanks to avoid S. chrysophrii infections. Upon arrival, 20 fish were randomly examined by a stereoscopical analysis of their two external gill arches to ascertain the absence of the parasite. One thousand fish, considered as the recipient group, were transferred from the hosting tanks into a net cage (6 m³) located within an inland cement tank (50 m³). The remaining fish were kept in identical conditions in an adjacent cement tank and constituted the control group. One week later, 250 gilthead sea bream (50 ± 13 g) naturally infected with S. chrysophrii (infection confirmed by stereoscopy but level not counted) were introduced into the cement tank harboring the recipient fish and were left to swim freely outside of (around and beneath) the net for the whole experiment’s duration. With this design, donor fish were not interfering with the feeding environment of the recipient fish. Fish tanks were continuously supplied with aerated flow-through sea water. All fish were hand-fed a commercial diet at a daily rate of 2–2.5% of body weight. The whole experiment lasted for 10 weeks from the time when donor fish were introduced in the system until the final recording of gill parasites. The water temperature during the experimental period ranged from 27 to 20°C (summer-fall). At the beginning of the experiment, fish were weekly monitored (random sample of 10 fish) and subjected to stereoscopical examination until it was confirmed that the recipient fish had received the parasite. This sampling ended after 2 weeks at the second sampling event when the observed prevalence was found to be 100%.

Traits

At the end of the experiment, all surviving recipients (n = 807) along with the control fish (n = 1,000) were individually weighed. The weighing of control fish was performed to evaluate the impact of infection vs. infection free environment on the growth. All recipient fish were slaughtered and subjected to stereoscopical examination (described above) for total counting of adult parasites of their 2 external gill arches. At the end of challenge test we had recordings on parasite count (PC), initial and final body weights (BW1 and BW2, weights at the beginning and at the end of the challenge test, respectively), and the specific growth rate (SGR) within the window (∼ 68 days) of the challenge test which was calculated using BW1 and BW2 information. The formula used for the calculations of SGR is given below.

$\text{SGR}=(\ln \text{BW}2-\ln \text{BW}1)/68$

Tissue samples of all the phenotyped fish were collected and stored in ethanol (100%) for further genomic work.

DNA Isolation, Library Preparation and Sequencing

The genomic DNA was extracted from the collected tissue samples of ∼ 20 mg. The genomic DNA extraction method and library preparation protocol are detailed in Aslam et al. (2018a). In-short, 935 (120 parents and 815 juveniles) 2b-RAD libraries were constructed following the protocol reported by Wang et al. (2012) with minor modifications. The digestion of genomic DNA (300 ng) was performed using AlfI (Thermo Fisher Scientific, United States) enzyme. The digestion step was followed by ligation of library-specific adaptors, sample specific barcoding, and the enrichment of 2b-RAD fragments using PCR amplifications. The individual specific libraries were then pooled using two different multiplexing strategies for parents (64 libraries per pool) and offspring (128 libraries per pool) with an objective to have higher depth of sequencing in parents. The pooled libraries were quality controlled using Agilent 2100 Bioanalyzer, and finally sequencing of pooled libraries was performed on an Illumina NextSeq500 platform (Illumina, San Diego, CA, United States) with 50 base single-end sequencing (v2 chemistry, high output kit – 50 cycles).

Genotyping 2b-RAD Loci

The genotyping from the individuals specific raw sequence reads was performed using following four steps,

(i) Adapter trimming and quality filtering was performed using custom developed scripts from the 2bRAD pipeline v2.0 (Wang et al., 2012). Poor quality reads were filtered out if mean Phred quality score within a sliding window of 4 bp was less than 15, and after the trimming and filtering, sequence reads of 36-bp length were kept.

(ii) The development of de novo reference sequence was performed by clustering of quality sequence reads of parents, which was accomplished with the use of Perl based scripts together with CD-HIT program (Fu et al., 2012; Wang et al., 2012).

(iii) Individual specific quality reads were then aligned to the developed reference sequence using bwa samse (V = 0.7.13-r1126; Li and Durbin, 2009) and the ‘mpileup’ function of SamTools version 1.2 (Li et al., 2009) was used to call variants and the call option of bcftools (Li et al., 2009) was used to call the genotype at each variant site for each fish.

(iv) The detection of reliable putative SNPs in the population was carried out using genotype quality (≥20), read depth (≥5), minor allele count (≥50) in the population at a particular site, and finally, minimum of 40 individuals were required to have genotype calls at the particular locus to fulfill mentioned criteria.

The SNP that passed the above-mentioned criteria were considered as a putative SNP for further analyses.

Parentage Assignments

The genotype data was further filtered to obtain a panel of highly informative markers for parentage assignments. Hence, the SNP genotype data was pruned for both parents and offspring using minor allele frequency criteria (MAF ≥ 0.35). The filtering criteria retained 916 highly informative SNPs which were then used to construct the pedigree with the likelihood ratio method implemented in CERVUS version 3.0 (Kalinowski et al., 2007). The assignments obtained from likelihood ratio method were verified using the opposite homozygote count method (Hayes, 2011; Ferdosi et al., 2014) with the full set of SNPs.

Genetic Linkage Mapping

The data was further pruned for the genetic linkage analysis, both at the marker and the individual levels. The filtering of markers was performed on, locus specific genotyping rate (≥70%), deviation from expected Mendelian segregation (P < 0.001), and the Hardy–Weinberg equilibrium exact test (P < 1.0 × 10^–7). Filtering of individuals was performed using individual specific genotyping rate (≥70%). An additional criterion was used for building linkage only where individuals or families were removed if they had less than five full-sibs per family. This criterion was used to include informative families with sufficient information on recombination, and to avoid computational problems in building the map.

Linkage mapping was performed using Lep-Map v2 (Rastas et al., 2015) with the analysis steps described in Aslam et al. (2018a). Briefly, the module “SeparateChromosomes” was used with minimum LOD threshold value of 19 to obtain linkage groups, “JoinSingles” module was used at LOD score limit of 5 in combination with LOD score difference of 2, and finally the module “OrderMarkers” was used to estimate the order and distance between the markers in centiMorgans (cM). The linkage maps are reported as sex averaged maps unless otherwise indicated and map figures were plotted using R package LinkageMapView (Ouellette et al., 2018).

Statistical Analyses

Descriptive Statistics

A generalized linear model was used in R to test the effects of recorded traits and to obtain initial evaluation of models. The growth traits (BW1, BW2, and SGR) did not show any significant effect on PC and therefore were not used as covariates. The initial body weight (BW1) showed a significant (P ≤ 0.001) effect on both SGR and BW2. The BW1 was not used as a covariate for BW2 to avoid genetic parameter deviations for BW2 due to potential collinearity between traits (BW1 and BW2) but it was used as a covariate for SGR, which should adjust for the same starting weight for all the fish and then obtain genetic variation for SGR. The raw data for PC presented a positively skewed distribution. Thus, to normalize this distribution a log transformation was performed on count values, i.e.,

$\text{LPC}={\log }_e(\text{PC}+1)$

where LPC is a log transformed parasite count, and PC is the raw parasite count. The addition of 1 to PC was performed to avoid taking logarithm of 0 (since some fish had PC = 0).

Variance Components Estimation

The heritabilities and genetic correlations for LPC, SGR, and BW2 were estimated using ASReml 4.0 (Gilmour et al., 2015) with genomic and pedigree-based relationship matrices (G-matrix and A-matrix, respectively) using bivariate (LPC vs. SGR and LPC vs. BW2) mixed models.

$y_{\text{i}}=\mathrm{\mu }+Xb+Zu+e(1)$

where y_i is a vector of ‘n’ records on trait i (LPC, SGR, or BW2), μ, b, and u are overall mean, vector of fixed effects of BW1 applied on SGR only, and the vector of additive genetic effects, respectively. The vector of breeding values is assumed to follow either

u∼MVN(0,G⊗G₀),oru∼MVN(0,A⊗G₀), where u = $\left[\begin{array}{c}{u}_1\\ u_2\end{array}\right]$, $G_0=\left[\begin{array}{cc}\hfill {\mathrm{\sigma }}_{u_{11}}^2\hfill & \hfill {\mathrm{\sigma }}_{u_{12}}^2\hfill \\ \hfill {\mathrm{\sigma }}_{u_{21}}^2\hfill & \hfill {\mathrm{\sigma }}_{u_{22}}^2\hfill \end{array}\right]$, G = genomic relationship matrix, A = pedigree relationship matrix, X = incidence matrix for fixed effects, Z = incidence matrix for additive effects, and e is the vector of random residuals with e∼MVN(0,I⊗R₀), where $R_0=\left[\begin{array}{cc}\hfill {\mathrm{\sigma }}_{{\text{e}}_1}^2\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\mathrm{\sigma }}_{{\text{e}}_2}^2\hfill \end{array}\right]$.

The genomic relationship matrix was constructed using VanRaden (VanRaden, 2008) method 1 as $\frac{P{P}^{\prime }}{2^{*}{\sum }_{i=1}^{Nsnp}{p}_i(1-p_i)}$; where P, Nsnp, and p_i are the matrix of centralized genotypes, total number of SNP markers, and the allele frequency of the reference allele, respectively. The narrow sense heritability was computed as the ratio of additive genetic variance to total phenotypic variance and genetic correlations were obtained from the above described bivariate analyses implemented in ASReml 4.0 (Gilmour et al., 2015).

Genome Wide Association Analysis (GWAS)

A univariate linear mixed model was used to perform genome wide association analysis for the LPC trait only. Model components were all the same as in model (1) except that the marker effects were also computed. The linear mixed model included random polygenic effects which accounted for family relationships (covariance between relatives) in the GWAS analysis. The marker – trait associations were performed using GCTA program (Yang et al., 2011) with the ‘–mlma-loco’ function which allows the estimation of a SNP effect by accounting the additive genetic variance expressed by all the markers distributed over all the linkage groups other than linkage group that contains the SNP.

The SNPs were considered as genome wide significant when they exceed the Bonferroni threshold for multiple testing P-value of P ≤ 2.29 × 10^–06 with –log₁₀(P) = 5.63, and if they surpassed P-value of P ≤ 5.74 × 10^–05 with –log₁₀(P) = 4.24 then they were graded as chromosome-wide significant. The significant values (chromosome and/or genome-wide) were computed as described in Aslam et al. (2020).

Quantile–Quantile (q–q) Plot

The q–q plot with the distribution of the observed –log₁₀(P-values) for each SNP, and the expected –log₁₀(P-values) from a theoretical distribution was plotted. The inflation factor (lambda, λ) was calculated using $\mathrm{\lambda }=\frac{\text{meadian}({\mathrm{\chi }}^2)}{0.456}$, chi-square statisticsχ² are obtained from significance-based P-values, and 0.456 is chi-square expected under the null hypothesis.

Estimation of the SNP(s) Variance

The proportional contribution of the chromosome and/or genome-wide significant SNP(s) to the total genetic variance was computed in ASReml 4.0 (Gilmour et al., 2015) as the reduction in the total genetic variance due to the addition of fixed effect of the SNP(s) in model (1). However, the G matrix used in this model was constructed with all other SNPs except the SNP(s) used as fixed effect in model (1).

Breeding Value Estimation

Breeding values were computed for LPC trait with univariate model (as described under the genome-wide association study) to quantify and compare accuracy of predictions obtained through pedigree (PBLUP) vs. genomic (GBLUP) information based BLUP models. The genomic relationship matrix was computed using all the quality markers (n = 21,773) left after filtering in genotyping process.

Accuracy of Prediction

For the estimation and the comparison of accuracies of prediction, we used records of individuals which had at-least one parent assigned and had genotype information along with the LPC record available. This limitation left us with 549 individuals available, and the accuracy of prediction for LPC was calculated as following.

A cross validation scheme was designed where ∼18% of the total 549 individuals were randomly masked. The random masking retained 449 individuals with information on LPC which were used as training set, while, randomly masked 100 individuals with missing phenotype records were used as validation animals. The breeding values were computed using genomic vs. pedigree relation matrices (Gvs.P). The mean accuracy of 20 replicates was computed as correlation (r_corr) of the estimated breeding value (pedigree/genomic) with the true phenotype, which were scaled by the square root of the heritability as $r_{\text{corr}}=\frac{\mathrm{\rho }(G\left[P\right]EBV,y)}{\sqrt{h^2}}$; where ρ, G[P]EBV, y, and h² are correlation coefficient, breeding values estimated using genomic or pedigree information, partially masked phenotypes, and genomic or pedigree based heritability estimates, respectively.

Results

Descriptive Statistics

The data statistics of recorded traits are presented in Table 1, together with the total number of observations for different traits (PC, LPC, BW1, BW2, SGR). There are some missing records due to some uncertain observation/counts. The distribution of parasite count and the log transformed parasite count is depicted in Supplementary File 1 (Supplementary Figures S1.1, S1.2).

TABLE 1

Table 1. Descriptive statistics for the recorded growth and parasite count traits.

Genotyping RAD Alleles

Approximately 12 runs of sequencing for RAD libraries using on NextSeq 500 platform yielded a total of 4.271 billion reads. The generated total amount of sequence reads were distributed as 1.217 and 3.054 billion across parents and offspring, respectively. The mean number of raw reads for parents were 10.06 (±2.306) millions, while the offspring had low average mean raw reads with a value of 4.22 (±2.464) millions. The trimming and quality filtering step of raw reads slightly reduced the number with a loss of 74.5 million (0.17%) reads which resulted in an average number of 9.864 (±2.458) and 4.215 (±2.301) million quality reads available for parents and offspring, respectively.

The catalog which was developed from the quality reads of parents comprised of 269,660 unique loci, that was used as a reference sequence. The process of SNP calling discovered 33,684 2b-RAD tags which revealed at-least one SNP identified. The SNP data was further pruned based on MAF (≥0.02) and the locus specific genotyping rate (≥30%) with the aim of increasing overall informativeness and decreasing the amount of missing or erroneous information. This pruning step retained 21,773 quality SNPs for further analyses. Individuals were also filtered out if the individual specific rate of genotyping was < 30%, which resulted to a drop of 94 individuals (91 offspring + 3 parents), and ultimately SNPs based genotype data comprised of 21,773 loci typed on 841 individuals (724 offspring 117 parents).

Parentage Assignments

The analysis of parentage assignments revealed 154 full-sib families with 1–22 sibs per family (Supplementary File 1 and Supplementary Figure S1.2). Out of the 724 offspring, 659 got both parents assigned, 62 got assignment of a single parent (sire or dam), while 3 offspring did not get any assignment. The validation of parentage assignment results using likelihood vs. opposite homozygotes count methods revealed 95% concordance. Out of the total 154 full-sib families, 45 families were used for the construction of linkage mapping which had a minimum of five sibs per family.

Linkage Map

The process of linkage mapping using 21,773 SNP markers in total, generated a linkage map comprised of 15,184 SNPs which were grouped into 24 linkage groups (SA01-SA24). The remaining 6,589 SNP markers did not get assignments to any group. The obtained total sex average map length was 1406.02 cM (Table 2, Supplementary Figure S1.3, and Supplementary File 2). The highest number of markers was found in linkage group SA02 with 754 markers, while the lowest number of markers were found in linkage group SA24 with 496 markers. The number of SNPs and the corresponding chromosome map length showed an average (for 24 linkage groups) correlation of 0.412. The sex specific maps (male vs. female) showed that the female genetic map length is larger than the males with total map length of 1,598.19 and 1,236.35 cM, respectively.

TABLE 2

Table 2. The genetic map of gilthead sea bream, Sparus aurata (SA) with sex specific and average map information.

The loci of assigned linkage groups were examined for their corresponding positions and order on physical map using alignment of tags to the genome assembly, Saurata_v1 (Pauletto et al., 2018). The results revealed more than 90% concordance in order of marker with coverage of each chromosome shown in Table 2. Moreover, linkage disequilibrium (LD) across markers from each linkage group were plotted using genotype data for parents only to present the distribution of markers along with the pattern of LD decay with increasing marker distances. The pattern showed decreasing trend of LD with increasing marker distances and the distribution of markers display the coverage of each chromosome (Table 2, Supplementary File 1, and Supplementary Figure S1.4).

Estimates of Variance Components and GWAS

The heritability estimates of LPC, SGR, and BW2 were very similar across the models (univariate vs. bivariate) and information sources (genomic vs. pedigree). The analysis using pedigree information resulted in heritability estimates of 0.147 ± 0.069, 0.155 ± 0.068, and 0.169 ± 0.072, for LPC, SGR, and BW2, respectively. Heritability estimates obtained with genomic information were 0.137 ± 0.061, 0.229 ± 0.077, 0.226 ± 0.074 for LPC, SGR, and BW2, respectively. The LPC trait showed a strong favorable negative genetic correlation of −0.807 ± 0.242 and −0.769 ± 0.224 (pedigree and genomic, respectively) with the SGR while −0.549 ± 0.269 and −0.701 ± 0.208 (pedigree and genomic, respectively) with BW2 (Table 3).

TABLE 3

Table 3. Estimates of variance components for PC and SGR using model 1 with pedigree vs. genomic information.

The genome-wide association analysis revealed a single SNP at linkage group SA17 which crossed the Bonferroni-corrected chromosome-wide significance threshold with a P-value of 4.02 × 10^–5 and allele substitution effect α of −0.374 for LPC (Figure 1 and Table 4). The allele substitution effects, and the minor allele frequencies for the top 5 SNPs are presented in Table 4. The proportion of the total genetic variance explained by the highest significant SNP was 22.68% estimated using the top SNP as fixed effect in the model.

FIGURE 1

Figure 1. Manhattan plot for distribution of P-values across different linkage groups. The unmapped markers are annotated as “UN” group in the plot. The red solid line represents the genome-wide while the blue dashed line displays the chromosome-wide significant Bonferroni thresholds. The SNP that crossed chromosome-wide Bonferroni threshold is highlighted green, and genes within ± 100Kb region of the top significant SNP are highlighted with arrows.

TABLE 4

Table 4. The statistics of the top 5 significantly associated SNPs with effect sizes.

Quantile–Quantile Plot

Quantile–quantile plot with the distribution of observed vs. expected P-values is depicted in Figure 2. The genomic inflation factor (lambda, λ) is a genomic control to get an estimate on false positives. The obtained λ-value from the fitted GWAS model with all markers was 1.023 (Figure 2). The observed P-values did not seem elevated in the q-q plot, which suggest that the GWAS analysis was not suffering from some unknown structure in the data, and thus also not from spurious associations. This is also confirmed by the value of λ= 1.023, which is close to 1.

FIGURE 2

Figure 2. QQ-Plot of P-values – Log Parasite Count.

Accuracy of Prediction

Prediction accuracies obtained using genomic information was 0.55 ± 0.165 which was higher than the accuracies obtained using pedigree information with estimates of 0.51 ± 0.147.

Discussion

Sparicotyle chrysophrii, a monogeneans gill parasite is one of the most severe threats for sea-cage production of gilthead sea bream in the Mediterranean. The infestation by S. chrysophrii affects all age groups in the sea which not only cause poor growth of fish but also mortalities that can go as high as 30% (Muniesa et al., 2020). The mortality at later stages when the fish is at the market size is likely to cause much more economic loss to the farmers. The treatment against this parasite includes chemicals such as formalin, hydrogen peroxide, and praziquantel which are very effective but are considered environment unfriendly and unsafe. Moreover, these applied available treatments affect consumer’s perception on product quality, ethics of production especially when there are increasing demands toward ecological products. The treatment cost against S. chrysophrii can go substantially high (€∼0.3 per kg fish produced) and the stocks may require multiple treatments in a production cycle. The cost of treatment for amoebic gill disease, which is caused by amoebic parasite, Paramoeba perurans accounts for 10–20% of the production cost in Australian farmed A. salmon (Nowak, 2012). Hence, treatments are not long-term remedy toward sustainable production. However, genetic improvement of commercial bream populations for resistance against S. chrysophrii (defined as low parasite count) using advanced selection methods may prove to be very useful, effective and preferred strategy especially when the trait is lowly heritable. Upgrading the overall robustness or resistance level against pathogens and/or parasites would avoid treatments using hazardous chemicals which should ultimately reduce environmental risks and promote sustainable production.

In the current study, we found significant difference in average weight (BW2) of challenged vs. control fish which was 43.93 g (SD = 11) and 58.01g (SD = 12), respectively, which suggests that the S. chrysophrii leads to a poor welfare of fish with effects on growth. Limited information is available for genetic variation of host resistance against S. chrysophrii parasite. However, there are several studies on genetic variation for resistance against different ectoparasites in Atlantic salmon and livestock species. In Atlantic salmon, resistance against amoebic gill disease where amoebic parasite colonizes the gill tissues and causes distress for the host is also reported to have a moderate heritability (Taylor et al., 2007; Kube et al., 2012; Robledo et al., 2018; Lillehammer et al., 2019). Similar estimates were obtained for host resistance to sea lice; ∼0.2 to 0.3 for the North Atlantic sea louse (Gjerde et al., 2011; Tsai et al., 2016), and 0.1–0.3 for the Pacific sea louse (Yáñez et al., 2014; Correa et al., 2017). Likewise, heritability estimates reported for resistance to common ticks in cattle were moderate (Shyma et al., 2015). Heritability estimates for resistance against S. chrysophrii gill parasite, i.e., LPC obtained in our study were low (0.147 and 0.137 with pedigree vs. genomic information, respectively) but significant, suggesting that the selection for improved resistance could be affective.

The analysis of genetic correlations showed favorable associations of LPC with both growth traits, i.e., SGR and BW2 (−0.549 to −0.807, respectively), depicting the decreasing trend of parasite counts toward heavier and fast growing families. Favorable relations among the traits should allow selection for growth and resistance to be performed simultaneously. Moreover, selection for growth should help to get correlated response for increasing resistance against S. chrysophrii. A similar favorable genetic correlation (–0.88) has been described for the severity gill disease due to Paramoeba perurans (gill parasite) and the body weight in Atlantic salmon (Gjerde et al., 2019). The reported genetic correlations for lice counts and growth traits have also been either close to zero or slightly favorable in Atlantic salmon (Tsai et al., 2016).

The current study produced a medium-density linkage map comprising 15,184 SNP markers distributed over 24 linkage groups. The generated numbers of linkage groups are consistent with the number of chromosomes/karyotype of gilthead seabream (Cataudella et al., 1980; Pauletto et al., 2018). The linkage map length in our study was 1,406.02 cM, which is shorter than what was obtained in previous studies by Tsigenopoulos et al. (2014) and Aslam et al. (2018a), where the sex-average map length was reported to be 1,970.29 and 3,899 cM, respectively. The difference in map lengths could also be explained by the difference in genome coverage, map resolution and population differences. However, the linkage map length reported by Palaiokostas et al. (2016) was ∼2.77 times longer than the map length obtained in the current study, which could be due to differences in adopted parameters and methods. Conversely, it was also observed that the high recombination rates involving false recombination detected due to data issues (e.g., genotyping errors, null alleles and nonrandom missing genotypes) may also cause the inflated genetic lengths (Jorgenson et al., 2005; Ronin et al., 2012).

Sex biased recombination has been observed in many fish species including Atlantic salmon, rainbow trout, zebrafish, etc., with significant reduction in recombination rate in heterogametic sex (Sakamoto et al., 2000; Singer et al., 2002; Moen et al., 2008). Mammal species, e.g., human, dog, pigs, etc. have also shown similar trend (Broman et al., 1998; Wong et al., 2010; Tortereau et al., 2012) where heterogametic male expressed lower recombination rates than the female. However, gilthead sea bream is a sequential protandrous hermaphrodite with majority of fish being able to produce either kind of gamete (sperm or ova) at different stages of their life. The S. aurata individuals in their first reproductive cycle mature as males (at the age of 2 years), though ovaries start differentiating in the larval stage but are replaced by the testes (Mylonas et al., 2011) as fish matures at 2 years of age, and then in the subsequent cycles, the testes regress and functional ovary develops in some males. Regardless of the hermaphrodite nature, we observed sex biased recombination in S. aurata with female-to-male recombination rates of 1.29:1.0 (Table 2) and the total female-specific map was 361.84 cM longer than the male specific map. The results obtained from our study on heterochiasmy were concordant with other studies (Franch et al., 2006; Tsigenopoulos et al., 2014; Aslam et al., 2018a). However, Palaiokostas et al. (2016) did not observe much difference in male vs. female recombination (1:1.05, female-to-male recombination rate) with a male specific map (4,010 cM) being slightly longer than the female map (3,822 cM). This discordance could be due to population specific differences, and/or different methods used for constructing the linkage maps. There are several other reports on the existence of heterochiasmy in hermaphroditic animals where females have shown a higher recombination rate than males (Theodosiou et al., 2016). The phenomenon of heterochiasmy in hermaphroditic species has been argued to be an effect of the female meiotic drive, which is considered as a mechanism in the evolution of neo-sex chromosomes (Yoshida and Kitano, 2012; Theodosiou et al., 2016; Abbott et al., 2017).

Our genome wide association analysis revealed no major QTL region(s) that reached the genome-wide significance threshold for LPC (Figure 1). However, there was a signal for a suggestive QTL with one SNP detected on linkage group 17 (SA17) that crossed the chromosome-wide Bonferroni threshold (Figure 1). This chromosome-wide significant SNP had a MAF of 0.136 with an allele substitution effect (α) of -0.374 in favorable direction (Table 4), which means that the substitution of allele should have a parasite reducing effect. Though, the second SNP of SA17 (Table 4) had low MAF (0.054) but showed an α value of similar magnitude to the top significant SNP of SA17 (Table 4). Three other SNPs, one from SA19 and the two from unknown (Un) linkage group also showed strong significance P-value but did not cross the threshold line (Figure 1 and Table 4). These SNPs and possibly multiple other loci may have crossed the threshold significance line with increased power (using more individuals with phenotype information) given that the association of these SNPs are not merely by random chance. As, these SNPs are distributed across different linkage groups and have shown some level of effect on trait, though, only one SNP is crossing suggestive chromosome-wide level of significance. From these results, one could argue about polygenic nature of trait with possibility of a few loci with relatively large and many loci with potentially small effects. In Atlantic salmon, resistance against parasitic diseases such as amoebic gill disease and sea lice resistance have also been reported to be very polygenic in nature (Tsai et al., 2016; Robledo et al., 2018). We detected a slight LD of 0.15 between the two markers of linkage group SA17 (Supplementary Figure S1.5). However, markers of unknown groups in Table 4 expressed relatively higher LD of 0.33 with each other but showed close to zero LD with markers of SA17 and SA19. The LD information among markers of different groups of Table 4 further strengthens the hypothesis of a polygenic nature of LPC, as close to zero LD for the markers of different linkage groups exempts the chance of misplacement of these markers during linkage mapping.

The λ value obtained from the P-values of SNPs from GWAS analysis/model was 1.023 (Figure 1), which was close to 1.0 indicating that the P-values are not inflated by any population stratification. Hence, no adjustments to the model was required when lambda is less than or equal to 1.0, as explained by Hinrichs et al. (2009).

The observed reduction in genetic variance was 0.022 (0.097–0.075) when an additional fixed effect of chromosome-wide significant SNP (Locus ID: 116287_33) was used in model 1. Hence, the proportion of genetic variance explained by the SNP was 22.68% which seems quite high compared to the QTL signal in the Manhattan plot (Figure 1). The tagged SNP (116287_33) might be highlighting an important QTL region explaining variation either directly or through LD with the true QTL. However, a seemingly inflated high proportion of genetic variation explained by the SNP could be a result of the Winner’s curse (Nakaoka and Inoue, 2010). We had only 549 recorded individuals with genotypes available for this complex parasite resistance trait and the top significant SNP (116287_33) did not pass the genome-wide threshold and therefore a validation of this QTL in an independent population is needed. The concordance of the results from validation study would signify the importance of the discovered QTL region and improve confidence for possibly performing efficient and economical QTL based selection for the industry. Moreover, the top 5 SNP markers of the current study that presented large effect sizes had relatively low MAF, which signifies the importance of rare mutations as explained in studies from Dickson et al. (2010) and Eichler et al. (2010) on synthetic associations and missing heritability due to rare mutations.

Cellular and molecular mechanisms underlying the host response to S. chrysophrii infection remain largely unknown (Sitjà-Bobadilla and Alvarez-Pellitero, 2009; Henry et al., 2015). Long-term infection apparently induces severe anemia, lower humoral response and impaired innate immune response (higher respiratory burst, but lower antimicrobial activities). The tag sequence (36 bp) of the highest significantly associated tag/SNP (116287_33) was aligned to the sea bream genome assembly, Saurata_v1 (Pauletto et al., 2018) to potentially obtain molecular insight which might help to understand possible causes in trait variation. The tag sequence aligned to chromosome 15 with a significant e-value = 4.00 × 10^–9. The region of approximately ± 100 Kb surrounding the highest significant SNP “116287_33” was searched for the underlying candidate genes which resulted in five genes (RFN122, XPO1, NFkB, SRF, and FAM160B1) encompassed in the candidate genomic region that have been reported to be involved in immune response in other species. RFN122 suppresses antiviral type I interferon (Wang et al., 2016) in virus-infected cells. Its potential role in the infection of an ectoparasite, however, needs to be determined. Exportin-1 (XPO1 or CRM1) is a nuclear export receptor involved in the transport of specific mRNAs. Inhibition of XPO1 was reported to down-regulated dendritic cell (DC) maturation (Chemnitz et al., 2010). As DCs are the most potent antigen-presenting cells, XPO1 differential expression/activity due to genetic variants might have broad effects on immune response. More generally, XPO1 has been shown to be involved in the regulation of NFkB, a key transcription factor in cell survival and immunity. Neurixn-1a is a single-pass type I membrane protein that is mainly expressed in the nervous system and there is no evidence so far of its involvement in immune function. Alternatively, serum response factor (SRF) is involved in many aspects of the immune system. Particularly relevant is SRF’s role in hematopoiesis and neutrophil migration during inflammation (Taylor and Halene, 2015), considering that in sea bream exposed to S. chrysophrii, increased activity of myeloperoxidase (MPO) was observed (Henry et al., 2015). Activity of MPO generally correlates well with levels of inflammatory responses sustained by neutrophilic granulocytes. Therefore, it might be possible that individual variation in SFR function has an effect on host-response to the parasite. Finally, FAM160B1 is a poorly characterized protein with no known involvement in immune response.

The accuracy of prediction obtained for LPC trait was fairly high with accuracy of 0.51 with pedigree and 0.55 using genomic information. Comparison of accuracies of predicting breeding values using PBLUP vs. GBLUP clearly shows that the genomic information based predictions are more accurate. Overall, we observed an ∼8% increase in accuracy with genomic information reflecting an expected advantage in genetic gain with genomic information. The complicated parasite resistance traits (e.g., lice and amoebic gill disease) in Atlantic salmon have also shown a similar trend with advantages of genomic over pedigree based predictions (Tsai et al., 2016; Robledo et al., 2018). The advantage of GBLUP over PBLUP is due to the fact that realized genomic-based relatedness between animals deviates from pedigree-based relationship coefficients. In addition, genomic breeding values are not affected by pedigree errors, though they can be affected by genotyping errors or introduction of mistakes during sample identifications. The current study had only rather shallow pedigree data available which might also have been a slight disadvantage for PBLUP.

The GS becomes even more beneficial compared to PBLUP when the trait is complex, difficult to record, lowly heritable and not recorded on the selection candidates. Additionally, GS can improve genetic gains through improving the accuracy of breeding values, and better use of variation within family (García-Ruiz et al., 2016; Zenger et al., 2019). Hence, in current era of genomics, the application of genomic selection is unavoidable in order to stay in industrial competition, to produce efficiently, and to follow market demand. The increase in accuracy of predictions (∼8%) for parasite count using genomics seem marginal with respect to the potential expenses on genotyping costs. However, the trait specific expected economic gain is very much dependent on the amount of production, production environment, given weight to the trait, and adopted strategy to develop production stock.

Conclusion

The current study shows that the 2b-RAD approach for the detection of SNPs and for the genotyping of a sea bream population is an effective method which resulted quality genotype data for ∼21K SNPs. The analysis reveals that the resistance against S. chrysophrii is significantly heritable (0.147 ± 0.069) and therefore genetic improvement is a valuable tool to reduce the prevalence of S. chrysophrii in farmed gilthead sea bream. The genetic architecture for host resistance to S. chrysophrii appears to be polygenic with one suggestive QTL detected at linkage group 17 which explains a large proportion (22.68%) of genetic variance. Since this SNP did not pass the genome-wide threshold, a validation of this QTL in an independent population is needed. The accuracy of predicting breeding values with genomic information was substantially higher (∼8%) than predictions using pedigree information. S. chrysophrii is a significant and increasing threat for Mediterranean sea bream production, and implementation of genomic selection in breeding programs may prove an efficient methodology to genetically improve host resistance to this parasite.

Data Availability Statement

Most of the data supporting the results are included within the manuscript as additional files. However, due to competitive nature and the privacy policy of commercial breeding companies, it is not permitted to share phenotypic, sequence, and genotypic data. The data can be made available for specific queries (on agreement) through the corresponding author.

Ethics Statement

The animal study was reviewed and approved by local ethical committee for the use of animals in the study. The reference code for the acquired approval is EL 25 BIO 37. Written informed consent was obtained from the owners for the participation of their animals in this study.

Author Contributions

MA analyzed the data and drafted the manuscript. RC performed the 2b-RAD library preparations and the other molecular laboratory-based analyses. AS, TM, CT, GR, LB, and KT were involved in the design of study, sampling, phenotype recording, and actively contributed in discussions. KT also assisted in arranging population material. KT, LB, CT, and GR played role in collection of data and coordination at the testing facility. AS and TM assisted in statistical data analysis. All authors read the manuscript, gave suggestions and comments for the improvement, and approved the final manuscript.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Funding

This work received funding from the European Union’s Seventh Framework Program (KBBE.2013.1.2-10) under grant agreement no. 613611 (FISHBOOST; www.fishboost.eu). The funding body played no role in the design of the study and collection, analysis and interpretation of data, and in writing the manuscript.

Acknowledgments

The results generated under this publication received funds from the European Union’s Seventh Framework Program (KBBE.2013.1.2-10) under the grant agreement no. 613611 (FISHBOOST). The genomic and the genetic analyses were performed on the Abel Cluster (national cluster of Norway) which was owned by the University of Oslo and the Norwegian metacenter for High Performance Computing (NOTUR), and ran by the Department for Research Computing at USIT, the University of Oslo IT-department (http://www.hpc.uio.no/). The part of contents of this manuscript was presented at the 11th World Congress on Genetics Applied to Livestock Production (Aslam et al., 2018b). Moreover, the lab based methodological procedures (2b-RAD library preparation) and the genomic analyses (sequence data processing and linkage mapping) performed in this manuscript were identical with our previous publication (Aslam et al., 2018a) which was also produced under the same project (FISHBOOST). We also thank reviewers of this manuscript for their comments and suggestions which have certainly improved the quality of the manuscript.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2020.594770/full#supplementary-material

Supplementary File 1 | Supporting information on results including information on parasite count trait, density of markers across different linkage groups and LD among markers.

Supplementary File 2 | The order of 2b-RAD tags on linkage map along with SNP positions and tag sequences.

Abbreviations

LPC, log transformed parasite count; h², heritability; SGR, specific growth rate; BW, body weight; GS, genomic selection; QTL, quantitative trait locus; SNPs, single nucleotide polymorphisms; RAD-Seq, restriction-site associated DNA sequencing; GWAS, genome wide association studies; QC, quality controls; LG, linkage groups; cM, centiMorgans.

Footnotes

1. ^ http://www.andromedagroup.eu/

References

Abbott, J. K., Nordén, A. K., and Hansson, B. (2017). Sex chromosome evolution: historical insights and future perspectives. Proc. R. Soc. B Biol. Sci. 284, 20162806.

Google Scholar

Aslam, M. L., Bastiaansen, J. W. M., Elferink, M. G., Megens, H.-J., Crooijmans, R. P. M. A., Blomberg, L. A., et al. (2012). Whole genome SNP discovery and analysis of genetic diversity in Turkey (Meleagris gallopavo). BMC Genomics 13:319. doi: 10.1186/1471-2164-13-391

PubMed Abstract | CrossRef Full Text | Google Scholar

Aslam, M. L., Boison, S. A., Lillehammer, M., Norris, A., and Gjerde, B. (2020). Genome-wide association mapping and accuracy of predictions for amoebic gill disease in Atlantic salmon (Salmo salar). Sci. Rep. 10:6435.

Google Scholar

Aslam, M. L., Carraro, R., Bestin, A., Cariou, S., Sonesson, A. K., Bruant, J.-S., et al. (2018a). Genetics of resistance to photobacteriosis in gilthead sea bream (Sparus aurata) using 2b-RAD sequencing. BMC Genet. 19:43. doi: 10.1186/s12863-018-0631-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Aslam, M. L., Carraro, R., Sonesson, A., Tzokas, K., Tsigenopoulos, C., Rigos, G., et al. (2018b). “Genetic basis of host resistance to S. chrysophrii in farmed gilthead sea bream (Sparus aurata) population”, Proceedings of the World Congress on Genetics Applied to Livestock Production, 2:39.

Google Scholar

Athanassopoulou, F., Ragias, V., Vagianou, S., Di Cave, D., Rigos, G., Papathanasiou, G., et al. (2005). Report of Sparicotyle (Microcotyle) chrysophrii Van Beneden and Hesse 1863, Atrispinum seminalis Euzet and Maillard 1973 and Polylabris tubicirrus Paperna and Kohn 1964 (Monogenea) on captive sea bream (Sparus aurata) and sharp snout sea bream (Diplodus puntazzo) in coastal Greece and Italy. Bull. Eur. Assoc. Fish Pathol. 25, 256–261.

Google Scholar

Baird, N. A., Etter, P. D., Atwood, T. S., Currey, M. C., Shiver, A. L., Lewis, Z. A., et al. (2008). Rapid SNP discovery and genetic mapping using sequenced RAD markers. PLoS One 3:e3376. doi: 10.1371/journal.pone.0003376

PubMed Abstract | CrossRef Full Text | Google Scholar

Bangera, R., Correa, K., Lhorente, J. P., Figueroa, R., and Yáñez, J. M. (2017). Genomic predictions can accelerate selection for resistance against Piscirickettsia salmonis in Atlantic salmon (Salmo salar). BMC Genomics 18:121. doi: 10.1186/s12864-017-3487-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Basurco, B., and Rogers, C. (eds) (2009). The Use of Veterinary Drugs and Vaccines in Mediterranean Aquaculture. Zaragoza: CIHEAM.

Google Scholar

Bishop, S. C., and Woolliams, J. A. (2014). Genomics and disease resistance studies in livestock. Livestock Sci. 166, 190–198. doi: 10.1016/j.livsci.2014.04.034

PubMed Abstract | CrossRef Full Text | Google Scholar

Broman, K. W., Murray, J. C., Sheffield, V. C., White, R. L., and Weber, J. L. (1998). Comprehensive human genetic maps: individual and sex-specific variation in recombination. Am. J. Hum. Genet. 63, 861–869. doi: 10.1086/302011

PubMed Abstract | CrossRef Full Text | Google Scholar

Campbell, N. R., Lapatra, S. E., Overturf, K., Towner, R., and Narum, S. R. (2014). Association mapping of disease resistance traits in rainbow trout using RAD sequencing. G3 4, 2473–2481. doi: 10.1534/g3.114.014621

PubMed Abstract | CrossRef Full Text | Google Scholar

Cataudella, S. P., Perin, R., and Sola, L. (1980). A chromosome study of eight Mediterranean species of Sparidae (Pisces, Perciformes). Genetica 54, 155–159. doi: 10.1007/bf00055985

CrossRef Full Text | Google Scholar

Chavanne, H., Janssen, K., Hofherr, J., Contini, F., Haffray, P., Komen, H., et al. (2016). A comprehensive survey on selective breeding programs and seed market in the European aquaculture fish industry. Aquac. Int. 24, 1287–1307. doi: 10.1007/s10499-016-9985-0

CrossRef Full Text | Google Scholar

Chemnitz, J., Turza, N., Hauber, I., Steinkasserer, A., and Hauber, J. (2010). The karyopherin CRM1 is required for dendritic cell maturation. Immunobiology 215, 370–379. doi: 10.1016/j.imbio.2009.05.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Correa, K., Lhorente, J. P., Bassini, L., López, M. E., Di Genova, A., Maass, A., et al. (2017). Genome wide association study for resistance to Caligus rogercresseyi in Atlantic salmon (Salmo salar L.) using a 50K SNP genotyping array. Aquaculture 472, 61–65. doi: 10.1016/j.aquaculture.2016.04.008

CrossRef Full Text | Google Scholar

Dickson, S. P., Wang, K., Krantz, I., Hakonarson, H., and Goldstein, D. B. (2010). Rare variants create synthetic genome-wide associations. PLoS Biol. 8:e1000294. doi: 10.1371/journal.pbio.1000294

PubMed Abstract | CrossRef Full Text | Google Scholar

Eichler, E. E., Flint, J., Gibson, G., Kong, A., Leal, S. M., Moore, J. H., et al. (2010). Missing heritability and strategies for finding the underlying causes of complex disease. Nat. Rev. Genet. 11, 446–450. doi: 10.1038/nrg2809

PubMed Abstract | CrossRef Full Text | Google Scholar

Everett, M. V., and Seeb, J. E. (2014). Detection and mapping of QTL for temperature tolerance and body size in Chinook salmon (Oncorhynchus tshawytscha) using genotyping by sequencing. Evol. Appl. 7, 480–492. doi: 10.1111/eva.12147

PubMed Abstract | CrossRef Full Text | Google Scholar

FEAP Secretariat (2017). FEAP, European Aquaculture Production Report 2008-2016. Brussels: FEAP Secretariat.

Google Scholar

Ferdosi, M. H., Kinghorn, B. P., Van Der Werf, J. H., Lee, S. H., and Gondro, C. (2014). hsphase: an R package for pedigree reconstruction, detection of recombination events, phasing and imputation of half-sib family groups. BMC Bioinform. 15:172. doi: 10.1186/1471-2105-15-172

PubMed Abstract | CrossRef Full Text | Google Scholar

Food and Agriculture Organization of the United Nations [FAO] (2014). National Aquaculture Sector Overview: Greece. Rome: FAO.

Google Scholar

Franch, R., Louro, B., Tsalavouta, M., Chatziplis, D., Tsigenopoulos, C. S., Sarropoulou, E., et al. (2006). A genetic linkage map of the hermaphrodite Teleost Fish Sparus aurata L. Genetics 174, 851–861.

Google Scholar

Fu, B., Liu, H., Yu, X., and Tong, J. (2016). A high-density genetic map and growth related QTL mapping in bighead carp (Hypophthalmichthys nobilis). Sci. Rep. 6:28679.

Google Scholar

Fu, L., Niu, B., Zhu, Z., Wu, S., and Li, W. (2012). CD-HIT: accelerated for clustering the next-generation sequencing data. Bioinformatics 28, 3150–3152. doi: 10.1093/bioinformatics/bts565

PubMed Abstract | CrossRef Full Text | Google Scholar

García-Ruiz, A., Cole, J. B., Vanraden, P. M., Wiggans, G. R., Ruiz-López, F. J., and Van Tassell, C. P. (2016). Changes in genetic selection differentials and generation intervals in US Holstein dairy cattle as a result of genomic selection. Proc. Natl. Acad. Sci. U.S.A. 113, E3995–E4004.

Google Scholar

Gilmour, A. R., Gogel, B. J., Cullis, B. R., Welham, S. J., and Thompson, R. (2015). ASReml User Guide Release 4.1 Structural Specification. Hemel Hempstead.

Google Scholar

Gjerde, B., Boison, S. A., Aslam, M. L., Løvoll, M., Bakke, H., Rey, S., et al. (2019). Estimates of genetic correlations between susceptibility of Atlantic salmon to amoebic gill disease in a bath challenge test and a field test. Aquaculture 511:734265. doi: 10.1016/j.aquaculture.2019.734265

CrossRef Full Text | Google Scholar

Gjerde, B., Ødegård, J., and Thorland, I. (2011). Estimates of genetic variation in the susceptibility of Atlantic salmon (Salmo salar) to the salmon louse Lepeophtheirus salmonis. Aquaculture 314, 66–72. doi: 10.1016/j.aquaculture.2011.01.026

CrossRef Full Text | Google Scholar

Gonen, S., Baranski, M., Thorland, I., Norris, A., Grove, H., Arnesen, P., et al. (2015). Mapping and validation of a major QTL affecting resistance to pancreas disease (salmonid alphavirus) in Atlantic salmon (Salmo salar). Heredity 115, 405–414. doi: 10.1038/hdy.2015.37

PubMed Abstract | CrossRef Full Text | Google Scholar

Gonen, S., Lowe, N. R., Cezard, T., Gharbi, K., Bishop, S. C., and Houston, R. D. (2014). Linkage maps of the Atlantic salmon (Salmo salar) genome derived from RAD sequencing. BMC Genomics 15:166. doi: 10.1186/1471-2164-15-166

PubMed Abstract | CrossRef Full Text | Google Scholar

Gray, I. C., Campbell, D. A., and Spurr, N. K. (2000). Single nucleotide polymorphisms as tools in human genetics. Hum. Mol. Genet. 9, 2403–2408. doi: 10.1093/hmg/9.16.2403

PubMed Abstract | CrossRef Full Text | Google Scholar

Haffray, P., Bugeon, J., Pincent, C., Chapuis, H., Mazeiraud, E., Rossignol, M.-N., et al. (2012). Negative genetic correlations between production traits and head or bony tissues in large all-female rainbow trout (Oncorhynchus mykiss). Aquaculture 368–369, 145–152. doi: 10.1016/j.aquaculture.2012.09.023

CrossRef Full Text | Google Scholar

Hayes, B. J. (2011). Technical note: efficient parentage assignment and pedigree reconstruction with dense single nucleotide polymorphism data. J. Dairy Sci. 94, 2114–2117. doi: 10.3168/jds.2010-3896

PubMed Abstract | CrossRef Full Text | Google Scholar

Henry, M. A., Nikoloudaki, C., Tsigenopoulos, C., and Rigos, G. (2015). Strong effect of long-term Sparicotyle chrysophrii infection on the cellular and innate immune responses of gilthead sea bream, Sparus aurata. Dev. Comp. Immunol. 51, 185–193. doi: 10.1016/j.dci.2015.03.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Hinrichs, A. L., Larkin, E. K., and Suarez, B. K. (2009). Population stratification and patterns of linkage disequilibrium. Genet. Epidemiol. 33, (Suppl. 1), S88–S92.

Google Scholar

Janssen, K., Chavanne, H., Berentsen, P., and Komen, H. (2017). Impact of selective breeding on European aquaculture. Aquaculture 472, 8–16. doi: 10.1016/j.aquaculture.2016.03.012

CrossRef Full Text | Google Scholar

Jorgenson, E., Tang, H., Gadde, M., Province, M., Leppert, M., Kardia, S., et al. (2005). Ethnicity and human genetic linkage maps. Am. J. Hum. Genet. 76, 276–290. doi: 10.1086/427926

PubMed Abstract | CrossRef Full Text | Google Scholar

Kakioka, R., Kokita, T., Kumada, H., Watanabe, K., and Okuda, N. (2013). A RAD-based linkage map and comparative genomics in the gudgeons (genus Gnathopogon, Cyprinidae). BMC Genomics 14:32. doi: 10.1186/1471-2164-14-32

PubMed Abstract | CrossRef Full Text | Google Scholar

Kalinowski, S. T., Taper, M. L., and Marshall, T. C. (2007). Revising how the computer program cervus accommodates genotyping error increases success in paternity assignment. Mol. Ecol. 16, 1099–1106. doi: 10.1111/j.1365-294x.2007.03089.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Kube, P. D., Taylor, R. S., and Elliott, N. G. (2012). Genetic variation in parasite resistance of Atlantic salmon to amoebic gill disease over multiple infections. Aquaculture 364, 165–172. doi: 10.1016/j.aquaculture.2012.08.026

CrossRef Full Text | Google Scholar

Li, H., and Durbin, R. (2009). Fast and accurate short read alignment with Burrows–Wheeler transform. Bioinformatics 25, 1754–1760. doi: 10.1093/bioinformatics/btp324

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N., et al. (2009). The sequence alignment/map format and SAMtools. Bioinformatics 25, 2078–2079. doi: 10.1093/bioinformatics/btp352

PubMed Abstract | CrossRef Full Text | Google Scholar

Lillehammer, M., Boison, S. A., Norris, A., Løvoll, M., Bakked, H., and Gjerde, B. (2019). Genetic parameters of resistance to amoebic gill disease in two Norwegian Atlantic salmon populations. Aquaculture 508, 83–89. doi: 10.1016/j.aquaculture.2019.04.072

CrossRef Full Text | Google Scholar

Lozano, D. (2019). Deliverable D3.1 – Development of a combined species SNP Array for the European Sea Bass and the Gilthead Sea Bream. Available online at: http://www.medaid-h2020.eu/index.php/2019/11/29/deliverable-d3-1/

Google Scholar

Meuwissen, T. H. E., and Goddard, M. E. (1996). The use of marker-haplotypes in animal breeding schemes. Genet. Select. Evol. 28, 161–176. doi: 10.1186/1297-9686-28-2-161

CrossRef Full Text | Google Scholar

Meuwissen, T. H. E., Hayes, B. J., and Goddard, M. E. (2001). Prediction of total genetic value using genome-wide dense marker maps. Genetics 157, 1819–1829.

Google Scholar

Moen, T., Hayes, B., Baranski, M., Berg, P. R., Kjøglum, S., Koop, B. F., et al. (2008). A linkage map of the Atlantic salmon (Salmo salar) based on EST-derived SNP markers. BMC Genomics 9:223. doi: 10.1186/1471-2164-9-223

PubMed Abstract | CrossRef Full Text | Google Scholar

Moen, T., Torgersen, J., Santi, N., Davidson, W. S., Baranski, M., Ødegård, J., et al. (2015). Epithelial cadherin determines resistance to infectious pancreatic necrosis virus in Atlantic Salmon. Genetics 200, 1313–1326. doi: 10.1534/genetics.115.175406

PubMed Abstract | CrossRef Full Text | Google Scholar

Muniesa, A., Basurco, B., Aguilera, C., Furones, D., Reverté, C., Sanjuan-Vilaplana, A., et al. (2020). Mapping the knowledge of the main diseases affecting sea bass and sea bream in Mediterranean. Transbound. Emerg. Dis. 67, 1089–1100. doi: 10.1111/tbed.13482

PubMed Abstract | CrossRef Full Text | Google Scholar

Mylonas, C. C., Zohar, Y., Pankhurst, N., and Kagawa, H. (2011). “Reproduction and broodstock management,” in Sparidae, Biology and Aquaculture of Gilthead Sea Bream and other Species, eds M. A. Pavlidis and C. C. Mylonas (New York, NY: John Wiley & Sons), 95–131. doi: 10.1002/9781444392210.ch4

CrossRef Full Text | Google Scholar

Nakaoka, H., and Inoue, I. (2010). The Winner’s Curse eLS. Chichester: John Wiley & Sons Ltd.

Google Scholar

Nowak, B. F. (2007). Parasitic diseases in marine cage culture – An example of experimental evolution of parasites? Int. J. Parasitol. 37, 581–588. doi: 10.1016/j.ijpara.2007.01.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Nowak, B. F. (2012). Neoparamoeba perurans, Fish Parasites: Pathobiology and Protection, eds P. T. K. Woo and K. Buchmann (Oxfordshier: CAB Interanational), 1–18.

Google Scholar

Ødegård, J., Baranski, M., Gjerde, B., and Gjedrem, T. (2011). Methodology for genetic evaluation of disease resistance in aquaculture species: challenges and future prospects. Aquacult. Res. 42, 103–114. doi: 10.1111/j.1365-2109.2010.02669.x

CrossRef Full Text | Google Scholar

Odegård, J., Moen, T., Santi, N., Korsvoll, S. A., Kjøglum, S., and Meuwissen, T. H. E. (2014). Genomic prediction in an admixed population of Atlantic salmon (Salmo salar). Front Genet. 5:402. doi: 10.3389/fgene.2014.00402

PubMed Abstract | CrossRef Full Text | Google Scholar

Ouellette, L. A., Reid, R. W., Blanchard, S. G., and Brouwer, C. R. (2018). LinkageMapView—rendering high-resolution linkage and QTL maps. Bioinformatics 34, 306–307. doi: 10.1093/bioinformatics/btx576

PubMed Abstract | CrossRef Full Text | Google Scholar

Palaiokostas, C., Bekaert, M., Davie, A., Cowan, M. E., Oral, M., Taggart, J. B., et al. (2013a). Mapping the sex determination locus in the Atlantic halibut (Hippoglossus hippoglossus) using RAD sequencing. BMC Genomics 14:566. doi: 10.1186/1471-2164-14-566

PubMed Abstract | CrossRef Full Text | Google Scholar

Palaiokostas, C., Bekaert, M., Khan, M. G. Q., Taggart, J. B., Gharbi, K., Mcandrew, B. J., et al. (2013b). Mapping and Validation of the major sex-determining region in Nile Tilapia (Oreochromis niloticus L.) Using RAD Sequencing. PLoS One 8:e68389. doi: 10.1371/journal.pone.0068389

PubMed Abstract | CrossRef Full Text | Google Scholar

Palaiokostas, C., Bekaert, M., Taggart, J. B., Gharbi, K., Mcandrew, B. J., Chatain, B., et al. (2015). A new SNP-based vision of the genetics of sex determination in European sea bass (Dicentrarchus labrax). Genet. Select. Evol. 47:68.

Google Scholar

Palaiokostas, C., Ferraresso, S., Franch, R., Houston, R. D., and Bargelloni, L. (2016). Genomic prediction of resistance to pasteurellosis in Gilthead Sea Bream (Sparus aurata) using 2b-RAD sequencing. G3 6, 3693–3700. doi: 10.1534/g3.116.035220

PubMed Abstract | CrossRef Full Text | Google Scholar

Palti, Y., Vallejo, R. L., Gao, G., Liu, S., Hernandez, A. G., Rexroad, C. E., et al. (2015). Detection and validation of QTL affecting bacterial cold water disease resistance in rainbow trout using restriction-site associated DNA sequencing. PLoS One 10:e0138435. doi: 10.1371/journal.pone.0138435

PubMed Abstract | CrossRef Full Text | Google Scholar

Pauletto, M., Manousaki, T., Ferraresso, S., Babbucci, M., Tsakogiannis, A., Louro, B., et al. (2018). Genomic analysis of Sparus aurata reveals the evolutionary dynamics of sex-biased genes in a sequential hermaphrodite fish. Commun. Biol. 1:119.

Google Scholar

Pecoraro, C., Babbucci, M., Villamor, A., Franch, R., Papetti, C., Leroy, B., et al. (2016). Methodological assessment of 2b-RAD genotyping technique for population structure inferences in yellowfin tuna (Thunnus albacares). Mar. Genomics 25, 43–48. doi: 10.1016/j.margen.2015.12.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Pérez-Sánchez, J., Naya-Català, F., Soriano, B., Piazzon, M. C., Hafez, A., Gabaldón, T., et al. (2019). Genome sequencing and transcriptome analysis reveal recent species-specific gene duplications in the plastic Gilthead Sea Bream (Sparus aurata). Front. Mar. Sci. 6:760. doi: 10.3389/fmars.2019.00760

CrossRef Full Text | Google Scholar

Puritz, J. B., Matz, M. V., Toonen, R. J., Weber, J. N., Bolnick, D. I., and Bird, C. E. (2014). Demystifying the RAD fad. Mol. Ecol. 23, 5937–5942. doi: 10.1111/mec.12965

PubMed Abstract | CrossRef Full Text | Google Scholar

Pyasatian, N., Fernando, R. L., and Dekkers, J. C. M. (2007). Genomic selection for marker-assisted improvement in line crosses. Theor. Appl. Genet. 115, 665–674. doi: 10.1007/s00122-007-0597-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Rastas, P., Calboli, F. C. F., Guo, B., Shikano, T., and Merilä, J. (2015). Construction of ultra-dense linkage maps with Lep-MAP2: stickleback F2 recombinant crosses as an example. Genome Biol. Evol. 8, 78–93. doi: 10.1093/gbe/evv250

PubMed Abstract | CrossRef Full Text | Google Scholar

Rigos, G., Fountoulaki, E., Cotou, E., Dotsika, E., Dourala, N., and Karacostas, I. (2013). Tissue distribution and field evaluation of caprylic acid against natural infections of Sparicotyle chrysophrii in cage-reared gilthead sea bream Sparus aurata. Aquaculture 408-409, 15–19. doi: 10.1016/j.aquaculture.2013.05.012

CrossRef Full Text | Google Scholar

Rigos, G., Henry, M., and Tsigenopoulos, C. (2015). Sparicotyle chrysophrii and gilthead sea bream – potential experimental infection model. Bull. Eur. Assoc. Fish Pathol. 35, 50–54.

Google Scholar

Robledo, D., Matika, O., Hamilton, A., and Houston, R. D. (2018). Genome-wide association and genomic selection for resistance to amoebic gill disease in Atlantic Salmon. G3 8, 1195–1203. doi: 10.1534/g3.118.200075

PubMed Abstract | CrossRef Full Text | Google Scholar

Ronin, Y., Mester, D., Minkov, D., Belotserkovski, R., Jackson, B. N., Schnable, P. S., et al. (2012). Two-phase analysis in consensus genetic mapping. G3 2, 537–549. doi: 10.1534/g3.112.002428

PubMed Abstract | CrossRef Full Text | Google Scholar

Sakamoto, T., Danzmann, R. G., Gharbi, K., Howard, P., Ozaki, A., Khoo, S. K., et al. (2000). A microsatellite linkage map of rainbow trout (Oncorhynchus mykiss) characterized by large sex-specific differences in recombination rates. Genetics 155, 1331–1345.

Google Scholar

Shyma, K. P., Gupta, J. P., and Singh, V. (2015). Breeding strategies for tick resistance in tropical cattle: a sustainable approach for tick control. J. Parasitic Dis. 39, 1–6. doi: 10.1007/s12639-013-0294-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Singer, A., Perlman, H., Yan, Y., Walker, C., Corley-Smith, G., Brandhorst, B., et al. (2002). Sex-specific recombination rates in Zebrafish (Danio rerio). Genetics 160, 649–657.

Google Scholar

Sitjà-Bobadilla, A. (2004). “Parasites in Mediterranean aquacultured fish: current impact and future research directions,” in Multidisciplinarity for Parasites, Vectors and Parasitic Diseases, ed. S. Mas-Coma (Bologna: Medimond), 301–311.

Google Scholar

Sitjà-Bobadilla, A., and Alvarez-Pellitero, P. (2009). Experimental transmission of Sparicotyle chrysophrii (Monogenea: Polyopisthocotylea) to gilthead seabream (Sparus aurata) and histopathology of the infection. Folia Parasitol. 56, 143–151. doi: 10.14411/fp.2009.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Taylor, A., and Halene, S. (2015). The regulatory role of serum response factor pathway in neutrophil inflammatory response. Curr. Opin. Hematol. 22, 67–73. doi: 10.1097/moh.0000000000000099

PubMed Abstract | CrossRef Full Text | Google Scholar

Taylor, R. S., Wynne, J. W., Kube, P. D., and Elliott, N. G. (2007). Genetic variation of resistance to amoebic gill disease in Atlantic salmon (Salmo salar) assessed in a challenge system. Aquaculture 272, S94–S99.

Google Scholar

Theodosiou, L., Mcmillan, W. O., and Puebla, O. (2016). Recombination in the eggs and sperm in a simultaneously hermaphroditic vertebrate. Proc. Biol. Sci. 283:20161821. doi: 10.1098/rspb.2016.1821

PubMed Abstract | CrossRef Full Text | Google Scholar

Tortereau, F., Servin, B., Frantz, L., Megens, H.-J., Milan, D., Rohrer, G., et al. (2012). A high density recombination map of the pig reveals a correlation between sex-specific recombination and GC content. BMC Genomics 13:586. doi: 10.1186/1471-2164-13-586

PubMed Abstract | CrossRef Full Text | Google Scholar

Tsai, H.-Y., Hamilton, A., Tinch, A. E., Guy, D. R., Bron, J. E., Taggart, J. B., et al. (2016). Genomic prediction of host resistance to sea lice in farmed Atlantic salmon populations. Genet. Select. Evol. 48, 47.

Google Scholar

Tsigenopoulos, C. S., Louro, B., Chatziplis, D., Lagnel, J., Vogiatzi, E., Loukovitis, D., et al. (2014). Second generation genetic linkage map for the gilthead sea bream Sparus aurata L. Mar. Genomics 18(Pt A), 77–82. doi: 10.1016/j.margen.2014.09.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Vallejo, R. L., Leeds, T. D., Gao, G., Parsons, J. E., Martin, K. E., Evenhuis, J. P., et al. (2017). Genomic selection models double the accuracy of predicted breeding values for bacterial cold water disease resistance compared to a traditional pedigree-based model in rainbow trout aquaculture. Genet. Select. Evol. 49:17.

Google Scholar

VanRaden, P. M. (2008). Efficient methods to compute genomic predictions. J. Dairy Sci. 91, 4414–4423. doi: 10.3168/jds.2007-0980

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, S., Meyer, E., Mckay, J. K., and Matz, M. V. (2012). 2b-RAD: a simple and flexible method for genome-wide genotyping. Nat. Meth. 9, 808–810. doi: 10.1038/nmeth.2023

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, W., Jiang, M., Liu, S., Zhang, S., Liu, W., Ma, Y., et al. (2016). RNF122 suppresses antiviral type I interferon production by targeting RIG-I CARDs to mediate RIG-I degradation. Proc. Natl. Acad. Sci. U.S. A. 113, 9581–9586. doi: 10.1073/pnas.1604277113

PubMed Abstract | CrossRef Full Text | Google Scholar

Wong, A. K., Ruhe, A. L., Dumont, B. L., Robertson, K. R., Guerrero, G., Shull, S. M., et al. (2010). A comprehensive linkage map of the dog genome. Genetics 184, 595–605. doi: 10.1534/genetics.109.106831

PubMed Abstract | CrossRef Full Text | Google Scholar

Yáñez, J. M., Bangera, R., Lhorente, J. P., Barría, A., Oyarzún, M., Neira, R., et al. (2016). Negative genetic correlation between resistance against Piscirickettsia salmonis and harvest weight in coho salmon (Oncorhynchus kisutch). Aquaculture 459, 8–13. doi: 10.1016/j.aquaculture.2016.03.020

CrossRef Full Text | Google Scholar

Yáñez, J. M., Lhorente, J. P., Bassini, L. N., Oyarzún, M., Neira, R., and Newman, S. (2014). Genetic co-variation between resistance against both Caligus rogercresseyi and Piscirickettsia salmonis, and body weight in Atlantic salmon (Salmo salar). Aquaculture 433, 295–298. doi: 10.1016/j.aquaculture.2014.06.026

CrossRef Full Text | Google Scholar

Yang, J., Lee, S. H., Goddard, M. E., and Visscher, P. M. (2011). GCTA: a tool for genome-wide complex trait analysis. Am. J. Hum. Genet. 88, 76–82. doi: 10.1016/j.ajhg.2010.11.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Yoshida, K., and Kitano, J. (2012). The contribution of female meiotic drive to the evolution of neo-sex chromosomes. Evol. Int. J. Organ. Evol. 66, 3198–3208.

Google Scholar

Zenger, K. R., Khatkar, M. S., Jones, D. B., Khalilisamani, N., Jerry, D. R., and Raadsma, H. W. (2019). Genomic selection in aquaculture: application, limitations and opportunities with special reference to marine shrimp and pearl oysters. Front. Genet. 9:693. doi: 10.3389/fgene.2018.00693

PubMed Abstract | CrossRef Full Text | Google Scholar