1 INTRODUCTION

Litopenaeus vannamei is a shrimp species that is cultured worldwide. It is believed to make up about 90% of shrimp aquaculture in the western hemisphere (Okpala et al., 2014) and greatly contributes to supplementing animal protein. Although several selective breeding projects for L. vannamei have been carried out (Andriantahina et al., 2013; Luan et al., 2015), broodstocks with high and stable growth traits and high resistance to disease are still urgently needed. Marker assistant selection (MAS) is an effective method for accelerating the selective progress in both plants and animals. With the development of highthroughput sequencing techniques, the high-density genetic linkage maps with single nucleotide polymorphisms (SNPs) covering the whole genome is now feasible in aquatic animals (Cui et al., 2015; Yu et al., 2015). For L. vannamei, advances in genome sequencing, high-throughput SNP discovery and high density linkage map construction (Yu et al., 2014; Yu et al., 2015) have established foundations for the genetic improvement of important traits based on MAS.

Genomic selection (GS) is one MAS strategy in which high-density markers covering the whole genome are used to predict the genomic estimated breeding value (GEBV) of each genotyped individual (Goddard and Hayes, 2007). GS is regarded as a powerful tool for complex traits, especially for low heritability traits, carcass quality traits, and disease resistance traits (Meuwissen, 2003). Compared with traditional breeding practices, GS can greatly accelerate the breeding process and reduce the costs of a breeding project because selection can be carried out at an early growth stage without the need for phenotypic measurements (Goddard and Hayes, 2007). As GS shows tremendous advantage over traditional breeding approach, it was promising to improve the accuracy of GS by using appropriate GS models.

Different GS models have been proposed to predict GEBVs. These models were established based on different assumptions about the distributions of marker effects. For example, the best linear unbiased prediction (BLUP) models, including RR-BLUP and GBLUP, assume that the effect of each marker follows a normal distribution with the same variance across all markers (Meuwissen et al., 2001; Habier et al., 2007). Bayesian models, such as BayesA and Bayesian LASSO, assume a different variance for each marker effect (Park and Casella, 2008). Previous reports based on simulations demonstrated that Bayesian models have a higher predictive ability than BLUP models (Solberg et al., 2008; Clark et al., 2011). Nevertheless, recent reports based on real data sets in livestock and plant breeding showed that the performance of these models can change depending on factors such as genetic architecture among traits and marker density (Moser et al., 2009; Heslot et al., 2012; Gao et al., 2013; Liu et al., 2014; Neves et al., 2014). Therefore, it is imperative to assess the performance of diverse models using real data sets, and thereafter, to identify the models that can provide improved prediction accuracy in a given population.

Although the application of GS in aquatic animals has been documented recently (Nielsen et al., 2009; Sonesson and Meuwissen, 2009; Nirea et al., 2012; Ødegård et al., 2014; Castillo-Juárez et al., 2015; Dou et al., 2016), no report appears to have employed real data sets for genomic prediction in shrimp. Considering the superior performance of GS when applied to livestock and crops (Zhang et al., 2011; Desta and Ortiz, 2014), GS is likely to be similarly applicable to shrimp species. Based on this proposal, we employed three statistical models (RR-BLUP, BayesA and Bayesian LASSO) to predict the GEBVs for two growth traits, body length and body weight to evaluate the performance of GS models in L. vannamei.

2 MATERIAL AND METHOD 2.1 Sample collection and genotyping

One full-sib family of L. vannamei was created in the breeding center of the Guangxi Fishery Institute (Fangchenggang, Guangxi Province, China). After about 100 days of culture from the post-larvae stage, 205 individuals were selected randomly and the body length and body weight of each individual were measured. The muscle tissue of each individual was collected and stored at -80℃ for genomic DNA extraction. The animals and genotyping method used in this study were the same as those reported previously (Yu et al., 2015). Briefly, all the individuals were genotyped based on the specific length amplified fragment sequencing (SLAF-seq) technique (Sun et al., 2013). A total of 6 359 SLAF markers which were genotyped successfully were collected. Those biallele SLAF markers could be considered as SNP markers and were chosen for further analysis. SNPs were then filtered based on either percentage of missing genotypes more than 20% or minor allele frequency less than 0.01 across the samples. Individual screening for the percentage of missing genotypes was also performed and any observations indicating more than 5% missing genotypes were discarded. Beagle 3.3.2 was used to impute missing genotypes (Browning and Browning, 2007).

2.2 Population stratification

A multidimensional scaling analysis was performed using the "cmdscale" function in the R software (R Development Core Team, 2014) to verify the genetic homogeneity of the data set. First, the matrix of genomic kinship (identity by state) was calculated with the GenABEL package based on the filtered markers (Aulchenko et al., 2007). Then the matrix of genomic kinship was input into the R software (R Development Core Team, 2014) and the multidimensional scaling was calculated using the "cmdscale" function.

2.3 Statistical models

Three statistical GE models were used to predict the GEBVs, namely RR-BLUP, BayesA, and Bayesian LASSO (Meuwissen et al., 2001; Park and Casella, 2008). The general linear BLUP model can be defined as:

where y is the vector of phenotypic observations; μ is the overall mean; 1_n is a vector of n (n=number of records); g is the vector of additive effects of SNPs with a variance of σ_g²; X is the corresponding design matrix with elements of X_ij=0, 1, and 2 for genotypes AA, AB, and BB, respectively for the i^th animal and j^th SNP; and e is the vector of residuals with a variance of σ_e² (Luan et al., 2009).

After the SNP effects were estimated by a model (RR-BLUP, BayesA, or Bayesian LASSO), the GEBV of each individual was computed for the validation of populations as:

RR-BLUP assumes that the marker effects are random (g~N(0, σ_g²)) (Meuwissen et al., 2001) and have equal variances . The variance parameters were unknown and were estimated using the restricted maximum likelihood method. The analysis was performed using the rrBLUP package (Endelman, 2011) in the R software (R Development Core Team, 2014).

The BayesA model assumes prior normal distribution with specific variance for each marker effect. The prior distribution of SNP variances follows a scaled inverse chi-square distribution with degrees of freedom df and scale parameter S. In this study, both the S and df of the distribution were set as described by Meuwissen et al. (2001). The Gibbs sampler was used for 40 000 iterations and the first 10 000 iterations were taken as burn-in (Ostersen et al., 2011). The analysis was implemented in the R software (R Development Core Team, 2014).

The Bayesian LASSO model assumes that the prior distribution assigned to the SNP effects is an identical, independent, double exponential density (DE) defined as , where λ is a sharpness parameter and λ² follows a gamma prior distribution . The residual variance (σ_e²) was assigned a distribution of scaled inverse chisquare with the degrees of freedom (df_e) and the scale parameter (S_e). In this study, we assumed that the residual variance accounted for one half of the sample variance of phenotype, and the parameters were set according to the method of de los Campos et al. (2013). This method was implemented using the BGLR package (de los Campos et al., 2009) in the R software (R Development Core Team, 2014) with 40 000 Markov Chain Monte Carlo (MCMC) iterations and with the first 10 000 iterations excluded as the burn-in (Ostersen et al., 2011).

2.4 Validation of the models

To systematically assess the performance of the GS models, a five-fold cross-validation approach was used in this study. Briefly, the data for each trait were divided randomly into five subsets, each containing 20% of the data. For each cross-validation experiment, one of the five subsets was retained as the validation set, and the other four subsets were combined and served as the training set. The process was repeated five times, each time with one subset as the validation set. Accordingly, each individual occurred only once in the validation set and had only one predicted GEBV (Luan et al., 2009).

2.5 Predictive ability of the GS models

The predictive ability of each GS model was assessed using the reliability of the GEBV and the bias of the predicted phenotypes. The reliability of the GEBV was evaluated by calculating the Pearson correlation between the GEBV of each shrimp sample in the validation population and their observed phenotypes. Bias was assessed using the regression coefficients of the observed phenotypes on the predicted phenotypes. A model was considered unbiased when the regression coefficient was close to one (Resende et al., 2012). The reliability of the GEBV and coefficient of regressions were calculated within each cross-validation analysis.

3 RESULT 3.1 Details of the phenotypes and markers

The average body length and body weight of the individuals used for the phenotype analyses were 117.13±7.36 mm and 20.32±3.72 g, respectively. After filtration, a total of 3 960 SNP markers in 44 linkage groups (LGs) were selected for the genomic prediction (Table 1). The average distance between adjacent SNP pairs within one LG was 1.04 cM, and the maximum distance between adjacent SNP pairs in each LG varied from 5.57 cM in LG14 to 36.80 cM in LG6. The number of SNPs in each LG also varied from 15 in LG14 to 194 in LG1

3.2 Population stratification

The multidimensional scaling analysis of an identity-by-state matrix for the 205 L. vannamei samples was shown in Fig. 1. The analysis revealed that all the points were not far from each other. Comparison of the multidimensional scaling variables (MDS1 and MDS2) showed that all the distributed points were well captured within the range -0.2 to 0.2. The points close to the borders were not considered as outliers, because they were relatively close to other neighboring points.

3.3 Predictive ability of the GS models

The reliability of the GEBV and the regression bias of the predicted phenotypes using three different GS models were shown in Table 2. Although these models used different prior assumptions of the marker effects, their predictive ability was very similar for the two growth traits in the analyzed L. vannamei population. The mean reliability of the GEBVs for body length and body weight using the three models was 0.296 and 0.411, respectively. For body length, BayesA performed slightly better than RR-BLUP and Bayesian LASSO, and the largest difference of 0.02 was between BayesA and Bayesian LASSO. For body weight, BayesA and Bayesian LASSO produced similar results that were 0.006 better than RR-BLUP. For both traits, the regression coefficients from the three models were close to one, indicating no obvious bias in the predictions. Additionally, no significant bias was detected among the different models.

4 DISCUSSION

In selective breeding projects for L. vannamei, body weight and body length are considered to be important selection indexes for growth traits. In this study, the average body length of the L. vannamei samples (about 117 mm) is comparable with the body lengths of 107 to 114 mm reported recently by Okpala and Bono (2016). However, the average body weight of the L. vannamei samples (about 20 g) is higher than the body weight of about 10 g reported in the Okpala and Bono (2016) study. Possibly, the body length and body weight are affected by sex, different culture environments, and different developmental stages (Okpala and Bono, 2016; Tan et al., 2016). It is worth noting that the selection strategy used in selective breeding projects may be a factor in the reported variations for growth parameters among different studies. It has been reported that the selection response (per generation) for the growth of L. vannamei can range from 2.3% to 25.0% owing to different selection strategies (Argue et al., 2002; Sui et al., 2015). Although traditional selective breeding methods can help to improve the genetic gains of such growth traits, they are still costly and time-consuming. Thus, there is a great need for improved breeding methods that can accelerate the genetic improvement of shrimp and decrease the costs.

GS has been shown to produce significant gains per unit cost and time in animal and crop breeding projects (Tayeh et al., 2015). The assessment of different GS models in various species has provided valuable information for the future application of GS. So far, no single model is guaranteed to produce best result across species or traits, because the performance of GS models depends on many factors, such as the architecture and heritability of the trait, the size of the training data set, the relatedness between the training population and validation population, and the marker density (Daetwyler et al., 2013).

In the present study, the three models that were tested showed similar prediction ability for two growth traits in shrimp. The regression coefficients of all the models indicated there was no obvious bias for the predictions. It has been reported that Bayesian models often have higher accuracy than BLUP models in the cases that a few quantitative trait loci (QTLs) account for a large proportion of the additive genetic variance (Verbyla et al., 2009; Erbe et al., 2012). BLUP-based models have been reported to have similar or higher accuracy than Bayesian models for traits that no known genes showing large effects (Luan et al., 2009; Silva et al., 2013). Therefore, we assumed that the growth traits of L. vannamei might be controlled by a large number of QTLs but with small effects. In a situation where the number of markers greatly outnumbers the number of genotyped animals, Bayesian models can take advantage of high-density markers and are expected to perform much better than BLUP-based models (Neves et al., 2014). In the present study, the low density of the markers might be a factor that contributed to the similar performance of the three models. Compared with BLUP-based models, Bayesian models are also able to efficiently capture the linkage disequilibrium information for a multi-breed reference population, but this advantage is reduced when there is a close relationship among the genotyped animals (Hayes et al., 2009b; Ostersen et al., 2011). For the L. vannamei population in the present study, all the individuals were from a full-sib family (Fig. 1) and this might potentially have affected the degree of performance of the Bayesian models.

In general, there appears to be increasing consensus that the heritability of a trait can affect the accuracy of the genomic prediction, and traits with high heritability can be predicted more accurately than those with low heritability (Daetwyler et al., 2010; Pszczola et al., 2012; Resende et al., 2012). In this study, we found the prediction reliability specific to body weight was significantly higher than the prediction reliability specific to body length for all three models, which might imply that body weight had higher heritability than body length in L. vannamei.

We believe that the current predictions produced by the three models could be improved by optimizing factors such as the density and distribution of markers, the size and structure of the training population, and the relation between training population and candidate population. Indeed, SNP density and distribution are crucial for the performance of statistical models. Solberg et al. (2008) and Weigel et al. (2009) found that significant gains in predictive ability were achieved by increasing the number of SNPs (Solberg et al., 2008; Weigel et al., 2009). Meuwissen (2009) suggested that 10N_eL markers were required to ensure the accuracy of genomic prediction, where N_e is the effective population size and L the genome size in Morgans. Because the L. vannamei genome is large (Yu et al., 2015), the current number of SNPs is too low to capture high level genetic variations. It has been reported that SNPs chosen randomly across the genome might affect the stability of genomic prediction (Spindel et al., 2015). Thus, an appropriate increase of marker density and a more even distribution of SNPs should improve the accuracy and stability of genomic prediction.

The accuracy of genomic predictions can also be affected by the size and structure of the training population (Hayes et al., 2009a; Moser et al., 2009). The size of a training population depends mainly on the N_e of a given population, and a population with high N_e will require a large number of individuals to build the training set (Clark et al., 2012). In this study, all individuals were from a full-sib family. As a result, the N_e was very small so the size of the training population could be low. Thus, the size of the training population in this study might have relatively little effect on the performance of GS models, even if it is lower than ideal. When all the samples are from a full-sib family, both linkage and linkage disequilibrium information can be useful for detecting gene effects, and polygenic effects can be captured (Lee et al., 2008). Intra-family predictions generally produce higher accuracy than predictions across families. Furthermore, the relation between training population and candidate population also can greatly affect the performance of GS models, and close relationships between the two populations can result in more accurate genomic predictions (van der Werf et al., 2010; Clark et al., 2012; Pszczola et al., 2012). We propose that the close relation among all the samples in this study resulted in the good prediction performances obtained for all three models. However, some limitations still exist in the use of a single fullsib family for the analysis. This study focuses on the evaluation of GS models based on a homogeneous population, and cannot reflect the effect of population structure, which plays a comprehensive role in model construction. In practice, GS has usually been applied to a broad-based population that may present some degree of population stratification. Therefore, further analyses based on a broad-based population would be required to systematically assess the performance of genomic prediction for L. vannamei.

5 CONCLUSION

Genomic prediction promises accurate selection in L. vannamei breeding programs. Three prediction models, RR-BLUP, BayesA and Bayesian LASSO, appear to be applicable for GS in shrimp. The direction of future studies should involve efforts to increase the number of training population and mining more SNPs that are evenly distributed across genome to improve the performances of GS models.