Genetic parameters and selection response for the harvest body weight of the giant freshwater prawn (Macrobrachium rosenbergii) in a breeding program in China

A multi-trait selective breeding program of Macrobrachium rosenbergii was initiated in China in 2015. In this program, the M. rosenbergii resources were widely collected from four countries, the origin of the founders was verified with 16 microsatellites and the pedigree was reconstructed, and the optimum contribution selection was used to make the mating design. In this study, we evaluated the genetic parameters and selection response for the harvest body weight (HBW) of M. rosenbergii after being communally reared for 95–109 days. The data were collected from two generations that comprised 25,212 progenies from 150 sires and 198 dams. The residual maximum-likelihood methodology was employed to evaluate the variance components, by fitting an animal model. The accuracy of estimated breeding values increased by 0.38% after pedigree reconstruction using microsatellite markers. The estimated heritability (h2) for HBW was moderate (0.212 ± 0.049) and the common environmental coefficient (c2) was low (0.063 ± 0.017) when all the data were used for the analysis. Within generations, h2 was moderate to high (0.198 ± 0.080 to 0.338 ± 0.049). c2 could only be estimated in G1, which was 0.055 ± 0.030. The average HBW of males was significantly larger than that of females (P < 0.01). h2 estimated for female HBWs were higher than that for males within generations, while h2 estimated for female HBWs were lower than that for males across generations. But they were not significantly different (P > 0.05). The genetic correlations between sexes were moderate to high within each generation (0.529 to 0.763). Two methods were used to estimate the realized response. One method was calculated from the differences between the least squares means of the selected population HBW and that of control population HBW, which was 14.01%. The other method was calculated from the differences between the EBVs of the selected population HBW and that of control population HBW, which was 11.52%. The predicted responses derived from two sets of genetic parameters acquired from within- and across- generation datasets were 11.68% and 10.67%, respectively. The present study provides valuable information for breeding programs of M. rosenbergii.

for seed industry; Special funds for major science and technology of breeding new agriculture (aquatic) varieties in Zhejiang province (2016C02055-2-1). The funder provided support in the form of research materials, but did not have any additional role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript. The commercial company: Jiangsu Shufeng Prawn Breeding Co. LTD did not provided financial support but did play a role in this study. The specific roles of all the authors are articulated in the 'author contributions' section. and high pond survival. There was likely some genetic connection among these individuals. The second strain was commercial larvae from Thailand Charoen Pokphand Group, characterized by fast growth in late stages. The last two strains were the progeny of wild populations from Burma and Bengal. A total of 134 prawns (56 females and 78 males) were used to construct the base population. All these strains were collected in 2015 and bred in Jiangsu Shufeng Prawn Breeding Industry Co., Ltd.

DNA extraction, amplification, and pedigree reconstruction
Pleopods were collected from 134 founder broodstocks and stored in 95% alcohol. Genomic DNA was extracted by the standard phenol and chloroform method [15]. Primer sets for 58 dinucleotide microsatellite loci from previously published literatures [16][17][18][19] or developed in our laboratory were used for amplification of alleles. PCR reactions were carried out in a 20 μl reaction mixture, containing 25 ng of template DNA, 5 pmol of each primer, 1.5 mM of MgCl 2 , 200 μM of each dNTPs, 1X Taq buffer, and 0.25 U of Taq DNA polymerase (TAKARA). PCR amplification was performed in a BIO-RAD T100™ thermal cycler with the following cyclic conditions: initial denaturation at 94˚C for 5 min, followed by 30 cycles at 94˚C for 30 s, 30 s for annealing at a specific temperature (Table 1), 45 s for elongation at 72˚C, and a final extension at 72˚C for 5 min. The PCR products were genotyped by ABI 3730. A total of 16 microsatellite markers were selected ( Table 1). The information of 16 loci in the 134 founders was shown in Table 2. The number of alleles per locus ranged from 7 to 22, with an average of 14 alleles per locus. The observed heterozygosity (H O ) varied from 0.120 to 0.806 while the expected heterozygosity (H E ) varied from 0.681 to 0.893. The polymorphic information content (PIC) varied from 0.630 to 0.879. The PIC measures were all larger than 0.5 in this study. Cervus 3.0 software was used to carry out the identity analysis for the 134 founder individuals with the likelihood-based method and the information of the 16 loci [20]. COLONY 2.0 [21][22][23] as implemented in R [24], was used to infer full-sib relationships according to the genotypes of candidate parents. COLONY uses a Bayesian approach to reconstruct the most likely full-sib groups. This analysis was carried out using the full-likelihood method. The 134 founders were assigned to 110 families, with 56 sires from 54 families, and 78 dams from 64 families.
Mating and production of families G 0 generation was established in 2016 by using an incomplete diallel crossing. The breeding design and the results of the family number for each breeding unit is shown in Table 3. Our previous experiments showed that the number of effective populations in the Bengal and Burma populations was too small, and the growth performance was always lower than that of the selected populations (Thailand and "Nantaihu No. 2"). Therefore, intra-population mating families and cross-breeding families between the two wild populations were no longer established. The survival rates of the two wild populations were better than that of the selected populations, so a certain number of hybrid families with other populations were designed. Similarly, the number of effective populations in Thailand is not large, so no intra-population mating families were established.
A selection index was used to rank all harvest shrimps from G 0 to G 1 . The relative weights were 70% for individual breeding values for harvest body weight (HBW), and 30% for family breeding values for pond survival. Herein, we focus on the genetic evaluation for the harvest body weights from G 0 to G 1 . Pond survival was estimated using a standard threshold (probit) and sire-dam model, which will be described in detail in next paper for publication. G 1 generation was established by using an optimal contribution selection based on a selection index, with a rate of inbreeding of 0.5% [25]. 72 mating groups were designed to generate G 1 . Within each mating group, one male was mated to 3-5 females from different families. The control population was constructed using 15 full-sib families, whose parents had the mean selection index in the G 0 generation and were from different full-sib families. Not all groups can be mated successfully according to the mating design. Sometimes only one pair could succeed in a group, which lead to the number of half-sib families being less than the number of candidate sires used in two generations.
Female and male candidate parents with a healthy appearance were selected, and each was reared in a cuboidal net cage with a side length of 0.25 m until fertilization. These cages were marked by a four-digit code, which was used as the individual ID. According to the mating design, 3-5 females and one male from different families were isolated and put in a separate 3 m 2 tank. These candidates were identified by unique color combinations of "Visible Implant Fluorescent Elastomers" (VIE) for each family. Each female prawn carrying fertilized eggs was put into a separate 170 L tank until hatching. Two thousand larvae were randomly selected from each tank and put into a 50 L culture tank. The family establishment periods ranged from 32 to 35 days in G 0 and G 1 (Table 4). A total of 78 families were produced in G 0 by  mating 56 males and 78 females, and 120 families were produced in G 1 by mating 90 males and 120 females (Table 5).

Rearing and family tagging
Larvae were fed twice a day with Artemia salina during the first 7 days after hatching. From the 8 th day, egg yolks was gradually added according to different developmental stages. Larvae were fed 4-8 times a day. Total amount and composition of the diet were adjusted daily according to different stages. The temperature was controlled at 30 ± 0.5˚C. Daily water exchanges were gradually increased, reaching 100% at the post-larval stage.
Six hundred post larvae with strong vitality were randomly selected from each family, and put into a 3 m 2 concrete tank until they grew to a body length of about 2.5 cm. Then the VIE tag was individually applied to identify each family. For G 0 , 150 prawns per family were randomly selected. For G 1 , 300 prawns per family were randomly selected. The selected prawns were assigned equally to each pond.
Pathogens, including Enterobacter cloacae, Citrobacter freundii, Aeromonas hydrophila, Enterobacter aerogenes, Nodavirus, Norovirus and Dicistrovirus, were detected by sampling and testing eggs and larvae at different stages from each family, using reverse transcriptasepolymerase chain reaction (RT-PCR). Families with pathogens were eliminated.

Growth and harvest
Three earthen ponds were used for communally rearing two generations (Table 4), P11 (2000 m 2 ) and P14 (1933 m 2 ) were utilized in the G 0 generation, and P10 (1933 m 2 ) and P11 (2000 m 2 ) were utilized in the G 1 generation. The temperature was 25 ± 6˚C during the grow-out stage. Juveniles in grow-out ponds were fed with a commercial pellet containing 37% crude protein, 6% crude fat, 2% crude fiber, 2% calcium, 1% total phosphorus, 1.8% lysine and 0.6% methionine. Feeding volume accounted for 2-3% of body weight, which were fed twice a day. After 95-109 days communally rearing the two generations (Table 4), a total of 25,212 individuals were harvested over the two generations (S1 Table).
All the surviving prawns were sampled, and for each of them the body weight, sex, family VIE code, pond, and harvest date were recorded. The phenotypes of the male prawns were classified. The claw of the male prawns contributes an important proportion of body weight. So, the male prawns with zero, one, or two claws were recorded as M0C, M1C, and M2C, respectively. Reproductive status is sometimes used as a common classification criterion for female prawns. While, the results of Hung et al. [26] suggested that body weights of female reproductive statuses were essentially the same traits and they could be analyzed together. This analytical model has also been used in many studies of M. rosenbergii [3,7]. So the females were not classified in this study. Individuals which were above average body weight and healthy, were transferred to a cuboidal net cage with a side length of 0.25 m until they were transferred to the 3 m 2 tanks for fertilization. A four-digit code label was placed above the cage to distinguish different individuals.

Data analysis
The results of the significance test of fixed effects and covariate (age) are shown in Table 6. All the main effects, the two-way interactions and covariates among them were statistically significant for all individuals (P < 0.001), except for interaction between generation and pond (P > 0.05), which was statistically significant for each sex (P < 0.001). Within each sex, different effects showed different degrees of influence, which was mostly because the classification criterions of the two sexes were different. When estimating variance components and genetic correlations between sexes, the reconstructed and unreconstructed pedigrees were used at the same time for comparisons. When calculating selective responses and making the mating design, only reconstructed pedigree analysis results were used.

Variance components and heritability estimation
Age at harvest was linearly related to HBW. Therefore, age was used as a covariate in the model. The variance components of HBW were evaluated with a restricted maximum likelihood (REML) approach in ASReml 4 [27]. The animal model was as follows: where y ijklmn is the observed HBW of the mth individual of the breeding population; μ is the overall mean; Gen i (the ith generation), Sex j (the jth gender), their interactions with Spec k (the kth sex specific phenotype, including M0C, M1C, M2C in males), and Pon l (the lth pond) is a linear covariate nested within the interaction among generation, sex and pond; a m is the additive genetic effect of the mth prawn, a( 0,As 2 a ), where A is the additive genetic relationship matrix among all prawns, including the parent prawns of the G 0 generation; f n is the random effect of the nth full-sib family, which is caused by the separate rearing of the full-sib families before communal rearing, f~(0,Is 2 c ), where I is an identity matrix; e ijklmn is the random residual error of the mth individual, e~(0, Is 2 e ). For within-generation analysis, the generation effect was excluded from the model. In the G 0 generation, the common environmental effect was removed from the model because of the convergence problem.
The genetic parameters before and after the pedigree reconstruction were estimated. The accuracy (r) of estimated breeding values before and after pedigree reconstruction was Where s m is the standard error of the mth individual, f m is the inbreeding coefficient of the mth individual, and s 2 A is the additive genetic variance obtained from the mixed model.

Genetic correlations between sexes
SPSS version 19.0 software was used to analyze the significant difference in HBW between males and females. Means of HBW between males and females within each generation were examined by an independent sample t-test. A bivariate analysis model was used to evaluate the genetic correlations between males and females within each generation, where the male and female body weights were regarded as two different traits. The bivariate analysis model was the same as the univariate model, the year and sex effects were excluded: where y klmn is the observed HBW of the mth individual with male or female gender; other fixed and random effects were the same as those in formula 1.
The phenotypic variance (s 2 p ) is obtained by adding additive variance (s 2 a ), common environmental variance (s 2 c ), and residual variance (s 2 e ) together. All variance components were estimated using within-and across-generation datasets. A complete pedigree from G 0 is available and used in all variance estimation. Heritability (h 2 ) was obtained from the s 2 a divided by s 2 p . The common environmental coefficient (c 2 ) was obtained from the s 2 c divided by s 2 p . Z-score was utilized to detect if there was significant difference in the heritability between males and females, and if the genetic correlations between different sexes were significantly deviated from one [29]: z ¼ xi À xj ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ðs 2 i þ s 2 j Þ q Where x i and x j are the heritabilities estimated for males and females; σ i and σ j are their standard errors, respectively. When detecting whether an estimate is significantly deviated from one, x j and σ j are set to one and zero, respectively.

Realized genetic gain estimation
Two methods were used to estimate the realized response. One method was calculated from the differences between the least squares means of the selected population HBW and that of control population HBW. The other method was calculated from the differences between the EBVs of the selected population HBW and that of control population HBW (Formula 1). The least squares means of the selected and control populations were evaluated by the following linear model: Where y ijklmn is the mth individual HBW in G 1 ; μ is the overall mean; Pop i (selected or control population), Sex j (the jth gender), Pon l (the lth pond) and their interactions with Spec k (the kth sex specific phenotype) (Sex j � Spec k , Sex j � Spec k � Pon l ) are fitted as fixed effects in G 1 , Age m (Sex j � Spec k � Pon l ) is a linear covariate, which is nested in the interaction among sex, sex specific phenotype and pond; Pop j � Fam n is the random effect of the nth family which is nested in the jth population; e ijklmn is the random residual error of the mth individual.

Predicted genetic gain estimation
Formula 1 was used to estimate the EBVs (Estimated Breeding Value) of all animals in G 0 and G 1 based on the best linear unbiased prediction (BLUP) by ASReml 4 software. The predicted genetic gain of HBW after one generation selection was achieved by calculating the difference of the average EBVs between G 1 and G 0 (including the control families in G 1 ). The selection response was achieved by calculating the difference between the least squares means of the selected population HBW and that of the control population HBW in G 1 .

Descriptive statistics
The observation numbers, simple means, minimum and maximum, standard deviation, and coefficients of variation (CVs) values for HBW in each generation are shown in Table 7. The mean body weight of the selected population of G 1 was smaller than that of the base population, G 0 . The relatively lower stocking and harvest density and longer days for grow-out test in G 0 may be responsible for this trend (Table 4). In G 1 , the mean and standard deviation of the selected population HBW were higher than that of the control population. The mean HBW in males was higher than that in females in each generation, ranging from 22.03% to 53.47%. The HBW CVs in the two generations were 29.54% to 33.07%, where males were also higher than Table 7. The number, mean and standard deviation for HBW of individuals with different phenotypic characteristics, by sex, population, and generation. females. The differences of HBW between males and females were significant in each generation (P < 0.01). In G 0 , the percentages of male specific characteristics were 2.62% for M0C (males with no claw), 10.24% for M1C (males with one claw), and 87.11% for M2C (males with two claws). In G 1 , the percentages of male specific characteristics were 6.69-7.80% for M0C, 21.99-23.69% for M1C, and 68.03-71.09% for M2C. A total of 22 males and 6 females had no record. For males, M2C had the highest mean body weight, followed by M1C and M0C.

Variance components, heritability, and common environmental coefficient
The results based on reconstructed pedigrees were very close to the results of the unconstructed pedigree, and there was no significant difference ( Table 7). The most obvious difference is that the variance components of females in G 0 could be estimated with reconstructed pedigree but not by unconstructed pedigree. This was mostly because the pedigree reconstruction increased genetic linkages among some founders. The accuracy of estimated breeding values after and before pedigree reconstruction was 0.633 and 0.631, respectively. After pedigree reconstruction, the accuracy improved by 0.38%.
Variance components, heritability, and the common environmental effect for HBW were estimated within-and across-generations. As shown in Table 8, the estimated heritabilities within each generation (0.338 ± 0.049, 0.198 ± 0.080, respectively) were moderate to high, and significantly deviated from zero (P < 0.05). The common environmental coefficient of G 1 was 0.055 ± 0.030. The common environmental variance in G 0 was negligible. The estimated heritability for the cross-generation dataset was moderate (0.212 ± 0.049). The common environmental coefficient was a little larger (0.063 ± 0.017) than within generations.
The HBW heritabilities of male and female prawns were moderate to high, ranging from 0.199 ± 0.089 to 0.491 ± 0.170. For within generation datasets, the estimated heritabilities for female HBWs were higher than that for males. For across generation dataset, the estimated heritability for female HBWs were lower than that for males. No significant differences were obtained in the heritability between male and female prawns within-and across-generations (P > 0.05). The common environmental coefficients were larger in females than in males, from both within-and across-generations. The genetic correlations between sexes were moderate to high within each generation (0.529 to 0.763). In G 0 , it was significantly different from 1 (z = -5.17, P < 0.01). In G 1 , it was not significantly different from 1 (z = -1.68, P > 0.05).

Selection response
The least squares means and EBVs of HBW in selected and control populations were shown in Table 9. The realized genetic gains after carrying out the selection for one generation were 14.01% and 11.52% by the above two methods, respectively. The predicted genetic gains of HBW per generation were calculated from two sets of genetic parameters. Set 1 was acquired from the across-generation dataset, while Set 2 was acquired from the mean value of two within-generation datasets ( Table 10).
The mean breeding values in G 1 were significantly increased compared to that in G 0 . After performing the selection for one generation, the genetic gains of HBW were 3.17 g and 3.47 g with parameters from Set 1 and Set 2, respectively. The percentage increases were 10.67% and 11.68%, respectively. The predicted genetic gains using Set 1 were a little lower than that predicted using Set 2.

Heritability and common environmental coefficient
After pedigree reconstruction, the estimated variance components were very close to the results of the unconstructed pedigree and the accuracy of estimating breeding value has been improved 3.8%, probably because there was not enough genetic relationship among founders  (Table 9).
themselves. The meaning for pedigree reconstruction is it would play an important role in avoiding mating among full-sib individuals.
The heritability for HBW of M. rosenbergii across generations was 0.212 ± 0.049, which was similar to the estimate of 0.22 ± 0.056 in India [4]. However, it was higher than the estimates of 0.14-0.15 in Vietnam [6,7] and 0.11 ± 0.08 in Thailand [29]. A previous selection of M. rosenbergii in China only achieved an HBW heritability estimation of 0.056, which was considered to be caused by low genetic variations of the foundations [3]. In this study, the base population was composed of four different strains from four different countries (China selected strain, Thailand, Burma, and Bengal populations), which could assemble greater genetic variations. The present data supported the opinion that low genetic variations of the foundations led to a low HBW heritability estimation of M. rosenbergii [3]. Besides, many other factors, such as different geographical populations, environmental conditions, and statistical models, would also affect the heritability estimation of the same species [4,30].
In G 0 , the common environmental variances could not be successfully estimated. It was most probably because there were not enough genetic ties among families. Lack of the common environmental coefficient (c 2 ) in the linear mixed model led to a higher heritability estimation in G 0 . c 2 was small but significantly different from zero across generations, which showed that the period of family construction between generations had a certain effect on HBW (Table 8). Notably, the standard error of heritability and the common environmental coefficient were small, which was possibly because of the large number of testing individuals in each family. Another issue that needed to be addressed was that aquaculture conditions should be fixed annually to obtain more accurate testing and evaluation results. The density in G 0 was very low (3.19-3.29 ind/m 2 ), which was far below the production density. So we increased the density in G 1 (9.75-14 ind/m 2 ). This would affect the social structure to some extent which could also be seen from the proportion of different male morphology in the two generations, but the impact on the analysis results was very limited. We analyzed the correlation among the initial stocking body weight, harvested body weight and survival rate, and found the results were very close in the two generations. It showed that the differences of density had little effect on growth, survival and other traits (S1 Fig, S2 Table).

Sexual differences of HBW
The mean HBWs of male and female M. rosenbergii were significantly different at each generation (Table 7). Males were 53.46% heavier than the females in G 0 and 19.39% heavier in G 1 (an average result from both selected and control populations). Previous studies also reported that male M. rosenbergii were much heavier, that is 62% [4] or 93% [6] higher than female prawns. Sexual size dimorphism is a common phenomenon in crustaceans, like Penaeus monodon, Fenneropenaeus chinensis, and Litopenaeus vannamei [31][32][33]. In M. rosenbergii, the average HBW of males is higher than that of females, which might be because females would allocate lots of energy to ovarian development and incubation [34,35]. The HBW CVs were 27.54-35.40% in males and 17.32-23.40% in females, which were similar to values reported by Luan,.13% in males and 6.57-34.54% in females) [3], though lower than those reported by Pillai, et al. (51% in males and 71% in females) [4]. The relatively lower HBW CVs is due to a relatively longer grow out period of M. rosenbergii used in the present study.
The HBW heritability of female M. rosenbergii (0.278 ± 0.067) was lower than that of males (0.309 ± 0.067) although there were no significant difference. This is different from previous studies in M. rosenbergii [3,4,6,7,29], which was most probably because the male prawns were classified in the genetic analysis, while female prawns did not. When the male and female HBWs were treated as separate traits, the genetic correlation was positive (0.529 ± 0.093 to 0.763 ± 0.142) and significantly different from one in the additional bivariate model, indicating that gender-related genes may have an effect on HBW. In addition, other factors such as behavioral factors, social interactions, and food deprivation might also affect the sexual dimorphism of HBW heritability [36].
The classification of sex related characteristics in M. rosenbergii has been reported in many studies. For males, the most common classification is divided into five distinct morphotypes: small claw, orange claw, blue claw, old blue claw, and no claw males, according to the claw color and the body size [37][38][39][40][41]. Genetic parameter evaluations of M. rosenbergii breeding populations were often based on the above classification [4,6,7,41]. Notably, according to our preliminary statistics, the claw weight of male M. rosenbergii was very high, and one claw accounted for about 10% of HBW. Males often face intensive social interactions in the process of communally rearing, which would lead to the loss of one or two claws. This will greatly affect the final harvest weights. In this study, one claw and no claw males accounted for 23.53% of the total number (Table 7). Therefore, the claw number was used as a classification method for males in this study.

Selection response
Although Gjedrem [42] summarized that the genetic gain for one generation was about 10% to 20% for aquatic animal species, the reported selection response of HBW of M. rosenbergii was not that high. In Vietnam, the selection response of HBW of M. rosenbergii reached an average rate of about 7% per generation [6,7]. In China, it was 6.56% per generation [3]. In the present breeding program, the realized response reached 14.01% or 11.52% after performing one generation selection, almost two times that of previous breeding programs. The predicted response was also over 10% after performing one selection. The obvious genetic progress is likely to benefit from the great genetic variation, which is the most important factor to determine the genetic progress of breeding objectives. Moreover, the optimum contribution selection method might also play an important role. Selecting and mating parents is quite important for a breeding program. The optimum contribution selection method provides a powerful tool to establish equilibrium between the genetic gains of the next generation and limit the inbreeding rate by restricting the increase in average co-ancestry [43]. This is as optimal contribution selection puts more selection pressure on Mendelian sampling term over truncated selection [44][45][46]. Additionally, selection theory points out that sustained genetic progress was from the creation of a covariance between Mendelian sampling terms and the long-term genetic progress of candidate parents [47,48]. The optimum contribution selection is more applied in terrestrial animals [49][50][51][52]. In aquatic animals, most studies have been carried out based on simulation and practical breeding projects [53][54][55].
In the present breeding program, increased density in G 1 would result in more social interaction, and due to the suppression of growth via social dominance, estimates of realized genetic gain based on least squares mean was a little higher than that based on breeding value and predicted genetic gain. The good consistency between them indicated that the current genetic information could exactly partition the additive genetic variance. Our results show that the optimal genetic contribution selection can also be used as an effective breeding strategy in aquatic animals.