Genomic selection for heterobothriosis resistance concurrent with body size in the tiger pufferfish, Takifugu rubripes

Parasite resistance traits in aquaculture species often have moderate heritability, indicating the potential for genetic improvements by selective breeding. However, parasite resistance is often synonymous with an undesirable negative correlation with body size. In this study, we first tested the feasibility of genomic selection (GS) on resistance to heterobothriosis, caused by the monogenean parasite Heterobothrium okamotoi, which leads to huge economic losses in aquaculture of the tiger pufferfish Takifugu rubripes. Then, using a simulation study, we tested the possibility of simultaneous improvement of parasite resistance, assessed by parasite counts on host fish (HC), and standard length (SL). Each trait showed moderate heritability (square-root transformed HC: h2 = 0.308 ± 0.123, S.E.; SL: h2 = 0.405 ± 0.131). The predictive abilities of genomic prediction among 12 models, including genomic Best Linear Unbiased Predictor (GBLUP), Bayesian regressions, and machine learning procedures, were also moderate for both transformed HC (0.248‒0.344) and SL (0.340‒0.481). These results confirmed the feasibility of GS for this trait. Although an undesirable genetic correlation was suggested between transformed HC and SL (rg = 0.228), the simulation study suggested the desired gains index can help achieve simultaneous genetic improvements in both traits.

Selective breeding is potentially able to boost aquaculture efficiency, now the fastest-growing food production industry 1 . In particular, pedigree-based breeding methods have contributed to aquaculture development by improving economically important traits, as seen in the salmonids and tilapias [2][3][4] . However, pedigree-based methods have innate drawbacks where it is assumed estimated breeding values (EBVs) of target traits for candidate individuals are the average breeding values of parents, ignoring Mendelian segregation within families 5 . Thus, pedigree-based methods can not differentiate EBVs among full sibs and large-scale pedigrees were needed to evaluate breeding values. On the other hand, molecular markers can be effective in handling Mendelian sampling by capturing genetic variance at DNA levels, i.e., full sibs could have different EBVs. By harnessing whole-genome high-density markers and advanced regression methods, Meuwissen et al. 6 proposed genomic selection (GS) to estimate the genomic estimated breeding values (GEBVs) of selection candidates. Thanks to the recent advances in genotyping by sequencing technologies, it is now affordable to genotype genome-wide single nucleotide polymorphisms (SNPs) for GS aquaculture breeding programs 7 . As expected, the greater performance of GS over the pedigree-based method in prediction and inbreeding control has been demonstrated by empirical studies using cultured fish populations 8,9 .
The tiger pufferfish Takifugu rubripes is a delicacy in Japan and is one of the most valuable marine fish species in Japanese aquaculture, ranking fourth in production value among cultured finfish 10 . To fulfill the growing demand for this species, selective breeding will be a practical approach to boost farming efficiency; however, tiger pufferfish aquaculture has not yet fully applied this technology 11,12 . Apart from growth-related traits, disease resistance should be highly valued in the breeding program, as disease outbreaks easily hamper the aquaculture industry. For instance, heterobothriosis, the gill disease caused by a monogenean parasite Heterobothrium okamotoi, severely threatens tiger pufferfish productivity and welfare 13 . The most severe infectious occurs at early OPEN 1 Fisheries Laboratory, University of Tokyo, Hamamatsu, Shizuoka 431-0214, Japan. 2 Veterinary Research Center, Nihon University, Fujisawa, Kanagawa 252-0880, Japan. * email: ahosoya@mail.ecc.u-tokyo.ac.jp

Results
Phenotypes. We produced test fishes by artificially crossing 11 wild males and 10 wild females, and subjected 240 4-month-old individuals to an artificial infection for 37 days. Heterobothriosis resistance was evaluated by counting the number of parasites attached to the branchial cavity walls (HC), and the standard length (cm) was measured on each fish (SL). The phenotypic mean was 15.85 (± 9.15 S.D.) for HC and 9.83 (± 0.78 S.D.) for SL ( Fig. 1 and Supplementary Table S1). As the plot shows, the distribution of HC was non-normal (Shapiro-Wilk test: p = 3.79 × 10 -6 , alpha level = 0.05) while SL approximated a normal distribution (Shapiro-Wilk test: p = 0.406, alpha level = 0.05). Therefore, we applied a square-root transformation on (HC + 1), approximating a normal distribution (Shapiro-Wilk test: p = 0.235, alpha level = 0.05). Transformed HC was used in the following genetic analysis. Weak but significant phenotypic correlation was observed between HC and SL (Pearson's r analysis: r = 0.157, p = 0.015; 95% confidence interval: 0.031 ≤ r ≤ 0.278).
Genotyping. We genotyped genome-wide SNPs of each individual using AmpliSeq 39 . The MiSeq sequencing generated an average of 174,870 (± 83,576 S.D.) raw reads per fish. After the quality-trimming step, the mean number of reads for each fish was 161,426 (± 83,576 S.D.) with the mean read length of 124 bp. The survived reads were mapped onto a reference fugu genome (FUGU5/fr3) for SNP calling. Following the quality filtration of SNPs, 6718 putative SNPs were yielded. Missing SNPs were imputed using LinkImputeR 40 . At this imputation step, 11 SNPs were discarded and 6707 imputed SNPs were called for each individual with the imputation accuracy of 0.888.    www.nature.com/scientificreports/ of the genetic correlation between the transformed HC and SL was also estimated. We detected a moderate antagonistic genetic correlation (r g = 0.228), where large individuals were suffering from higher parasitic loads. This genetic correlation could be due to the phenotypic correlation, although phenotypic correlation between HC and SL was weak as described above. Therefore, we tested correlation between GEBV for each trait using a univariate linear model (i.e. GBLUP); to predict GEBV for HC, SL was included as the covariate to minimize non-genetic effects from SL. If genetic correlation exists between the two phenotypes, the GEBVs would also show a correlation. We found positive correlation (Pearson's r = 0.252, p = 7.67 × 10 -5 ). This supports genetic correlation between the two traits.

Model comparison of genomic prediction.
To examine the availability of GS for HC and SL, we applied 12 regression models: i.e., GBLUP, Bayes A, Bayes B, Bayes C, Bayes Ridge, Bayes LASSO, Bayesian reproducing kernel Hilbert space (Bayesian RKHS), support vector machine regression with a linear kernel (SVR-linear), SVR with a poly kernel (SVR-poly), SVR with a radial basis function kernel (SVR-rbf), feedforward neural networks (FNN), and multi-task feedforward neural networks (multi-task FNN). We compared predictive ability defined as Pearson's r between the GEBVs and observed phenotypes by means of a tenfold cross-validation scheme. Predictive abilities for transformed HC ranged from 0.248 to 0.344 under 12 models (Table 1). Among these models, SVR-poly and SVR-rbf models were inferior, while two deep learning models were slightly better. On the contrary, the two SVR based models ranked at the top for prediction of SL, and deep learning models were inferior. Bayes RKHS and GBLUP models showed good performance in both traits. In this study, we tested the availability of LGSI methods for simultaneous improvements of HC and SL, assuming a genetic correlation estimated above (r g = 0.228), using simulation studies. We simulated six scenarios each different in selection schemes, i.e. random mating (RAND), GS on HC only (GS HC ), GS on SL only (GS SL ), selection applying Smith-Hazel index 28,29 (S1 SHI and S2 SHI , different in economic weights) and desired gains index 31 (S DGI ), for 10 generations with 50 replications (Fig. 4). In short, RAND was based on random mating while GS HC and GS SL were based on GS on either of the traits. GEBV was estimated by GBLUP. In S1 SHI , selection candidates were ranked based on the Smith-Hazel index. Since economic importance

Recurrent selection Initiation
Random sampling (F i , n = 2,000) Ten generations (i = 1 to 10)  The founder population (n = 10,000) was constructed, and 20 sires and 20 dams were randomly sampled to produce 8000 progenies. Then, 2000 fish were randomly sampled from the progeny pool as the broodstock population (F 0 ). (b) The workflow of recurrent selection schemes. Parents (20 sires and 20 dams) were selected from F 0 according to the scenario-specific selection criteria and 8000 progenies were generated. The selection scenarios were: RAND, random selection; GS HC , selection on Heterobothrium okamotoi counts (HC); GS SL , selection on standard length (SL); S1 SHI and S2 SHI , selection based on genomic Smith-Hazel index (SHI); S DGI , selection based on the desired gains index (DGI). S1 SHI has the same economic weights for both traits, and S2 SHI uses the similar vector of economic weights as the vector of desired gains in S DGI . Then, random sampling was applied to select 2000 progenies as the broodstock population for the next generation. A total of 10 generations (F 1 to F 10 ) of this process were replicated 50 times.

Scientific Reports
| (2020) 10:19976 | https://doi.org/10.1038/s41598-020-77069-z www.nature.com/scientificreports/ for each trait has not been evaluated in the tiger pufferfish aquaculture industry, we assume both traits have equal economic weights, which is w = [− 1, 1] for HC and SL for S1 SHI (HC is expected to decrease by selection). For S DGI , d was set as [− 3, 0.3] for HC and SL, so that SL can be improved preferentially while HC can be reduced by 30% after 10 generations (− 3 * 10/100 = − 30%). To compare the two selection index methods, we also ran an additional scenario (S2 SHI ) based on Smith-Hazel index, where the economic weight for each trait was set the same as the designed weights of S DGI (w = [− 3, 0.3]). As expected, only S DGI could improve the two traits simultaneously, where true breeding values (TBVs) of parasite load (HC) decreased while SL increased in each generation (Fig. 5).

Discussion
In this study, we tested the possibility of GS for genetic improvements in heterobothriosis resistance of the tiger pufferfish from empirical data and conducted a simulation study to design a GS breeding strategy that could improve the resistance trait concurrent with growth-related traits. Overall, our results suggest GS for the parasite resistance trait is feasible (predictive ability = 0.248-0.344) and breeding strategy incorporating the DGI method can simultaneously improve both HC and SL, even though an unfavorable antagonistic genetic correlation was suggested (r g = 0.228).
With 6707 SNP makers, moderate estimated heritability of transformed HC (h 2 = 0.308, SE = 0.123) and SL (h 2 = 0.405, SE = 0.131) were obtained, indicating selective breeding for those traits is feasible. The estimated heritability was comparable to those estimated for resistance against sea lice in Atlantic salmon (h 2 = 0.22 to 0.33 with 35 k SNPs) 22 and bacterial cold water disease resistance (survival days) in farmed rainbow trout (h 2 = 0.33     9 . This suggested our small SNP panel could successfully capture the genetic variance for HC in the tiger pufferfish. In this study, we could not detect significant SNPs from GWAS. Even with the small SNP panel and small sample size, strong effect SNPs (the sex-determining SNP) could be detected in a cultured population of the tiger pufferfish 39 . Therefore, our GWAS result suggests the parasitic resistance is controlled by a large number of quantitative trait locus (QTL) with small or moderate effects, and marker-assisted selection is not feasible. This result is consistent with the previous QTL analysis using the interspecies hybrid system of pufferfishes 19 . Although the effect was not significant, genes neighboring the SNP with highest -log 10 (p) values on chromosome 9 (12,024,615 bp) deserve further investigation, because the genomic region including this site was reported to have a small QTL effect on host specificity of H. okamotoi 19 . The predictive abilities for HC estimated under 12 models were moderate (0.248-0.344), and within the range observed for other disease resistance traits examined in other fish species 21,23,24,42 . The predictive abilities of Bayesian hierarchical linear models (i.e. Bayes A, B, C, LASSO, and Ridge) were similar (0.303-0.312) and scarcely higher than the GBLUP model (0.307 ± 0.018) for HC. This suggests that these linear models did not greatly differ regarding the predictive ability and the assumptions of the prior distribution of genetic effects have a limited impact on this trait. Bayes RKHS showed slightly better performance in HC compared to these linear models. For SVR-poly and SVR-rbf models, relatively low abilities for HC were observed, however, high abilities were found for SL. Since the default hyperparameters were used in the SVR models, hyperparameter tuning may aid achievement of better performance for HC as in the case of the previous study 43 . The architectures of FNN and multi-task FNN were tuned to achieve high predictive ability of GS for HC, however, the same architecture was applied to calculate the predictive ability of GS for SL. As expected, these models resulted in high predictive ability for HC but low for SL. This indicates that a deep learning model is task-specific and high accuracy can be obtained with careful optimization as described previously 44 . However, a great improvement in predictive ability was not achieved by FNN methods compared to GBLUP and Bayesian models even with the model complexity.
Our simulation study showed the availability of DGI for simultaneous genetic improvement in HC and SL even when the unfavorable antagonistic genetic correlation was assumed. The two scenarios incorporating the Smith-Hazel index showed the undesired consequences, where the average TBV for both SL and HC increased (smaller HC is favored). This happened because the breeding scheme only selected the individuals with the top LGSI values, but the high LGSI calculated by the Smith-Hazel index method does not guarantee the selected individuals are superior in both of the traits 45 , especially when target traits show a negative correlation. On the other hand, DGI, a variation of the selection index methods, allows selection with restrictions on multiple traits via the desired gains vector (d). In this study, the d vector was set with intending to reduce HC by 30% during 10 generations while maximizing SL. The desired gains vector (d) can be further optimized by comparing simulation scenarios with various d to achieve the self-defined breeding goal. Unfavorable genetic correlation between body size and disease resistance is commonly observed in aquaculture species, e.g. vibriosis in Atlantic cod 46 , bacterial cold water disease in rainbow trout 47 , and piscirickettsiosis in coho salmon 48 . Therefore, it is expected that DGI or the similar LGSI method can be widely applied for the simultaneous improvement of disease resistance trait and growth-related traits, which are the primary targets of most breeding programs.
In summary, the availability of GS for HC and SL was confirmed in this study. Moderate heritability for both traits suggests the genetic return from GS is high. GBLUP and Bayesian linear regression models showed similar predictive abilities for these traits. Although an unfavorable antagonistic genetic correlation was suggested between the two traits, the GS breeding strategy incorporating DGI can be a solution for the simultaneous genetic improvement.

Methods
Sample fish. The empirical experiments were performed in the Fisheries Laboratory, University of Tokyo (Hamamatsu, Shizuoka Prefecture, Japan). All samples (n = 240) were generated by a full-factorial mating among 10 wild males and 11 wild females, which were commercially caught from Wakasa Bay (Fukui Prefecture, Japan). For the mating, artificial fertilization was applied following the previous study 12 with minor modification. In brief, females were anesthetized with 200 mg/l of 2-phenoxyethanol and then ripened by injection of 150 µg/kg of luteinizing hormone-releasing hormone (LHRH, Sigma-Aldrich, St. Louis, MP, USA). Gametes were stripped from each individual and fertilized per male-female pair (110 pairs in total). Fertilized eggs of each maternal half-sib family were mixed and kept in a hatching jar. After hatching, each maternal half-sib was kept in a holding tank for 1 month and then all families were mixed and cultured in a three-ton communal tank. Rearing and feeding conditions were set as previously described 10 . At 4 months age, 240 fish were randomly collected and subjected to the artificial challenge test.
Artificial infection and phenotyping. Artificial infection was done following previous studies 12, 49 . A day before the infection, fish were equally distributed into three identical one-ton experimental tanks (80 individuals/tank) supplied with H. okamotoi-free fresh seawater (UV treated and filtered). Meanwhile, eggs of H. okamotoi were collected from tanks containing infected fish and kept in a glass jar containing fresh seawater until infection. Hatching was induced by physical stimulation (shaking at 140 rpm for 10 min) and the density of oncomiracidia, the free-living larvae of H. okamotoi, in the suspensions was determined under the microscope just before the infection. At infection, the water depth of experimental tanks was adjusted to 15 cm, and approximately 4000 oncomiracidia was introduced into each tank. At 3 h post-exposure, fish were transferred into three, newly-setup one-ton holding tanks and reared for 32 days, when H. okamotoi reaches maturation and moves to the branchial cavity walls (BCW) 13 . At the 32-day mark, fish were euthanized, measured for SL and the BCWs dissected from both sides. For each fish, the caudal fin was clipped and kept in 600 µl TNE8U buffer 50 ( Genotyping. Genomic DNA was extracted using a Gentra Puregene Tissue Kit (QIAGEN, Hilden, Germany) following the manufacture's instruction and applied for AmpliSeq genotyping as previously described 39 .
In short, 3187 genome-wide target regions were amplified by the first PCR with the custom AmpliSeq primer pools. After the adapter ligation and purification steps, PCR products were barcoded by a second PCR with 8-base index oligo sequences (Supplementary Table S1) for individual identification. The libraries of 240 individuals were pooled and sequenced on Illumina MiSeq System using the MiSeq reagent kit v2 (300 cycles) from both ends. The raw FASTQ reads were quality-trimmed using trimmomatic-0.36 51 40 . At first, SNPs which did not fulfill the maximum missingness per SNP and sample of 0.9 were filtered out, and then the missing genotypes were imputed. All samples were retained but 11 out of 6718 SNPs were discarded. Subsequently, the imputation accuracy was accessed with numbermasked option (set as 500). PLINK (v1.0.7) 57 -recodeA option was used to generate the allele coding matrix from the imputed VCF files.
Population structure. The population structure was analyzed with a nonlinear dimension reduction technique, t-distributed stochastic neighbor embedding (t-SNE) 41,58 . At first, t-SNE transforms the genomic data into conditional probabilities that represent pairwise similarity in the high-dimensional space. Then, transformed data were applied to a heavy-tailed Student's t-distribution that measures pairwise similarities of corresponding samples in the low-dimensional embedding space. Finally, t-SNE minimized the sum of the Kullback-Leibler divergence between those two probability distributions (Kullback-Leibler divergence is the measure of the difference between two probability distributions). The t-SNE analysis was implemented in sklearn.manifold. TSNE class of Python/scikit-learn-0.20.3. The perplexity was set as 20 and default values were used for the other parameters.
Heritability and genetic correlation. Heritability and genetic correlation were calculated by a multivariate linear mixed model as follows: where y i is a vector of phenotypes for trait i (i = 1 for transformed HC and 2 for SL); X i and Z i are incidence matrices for fixed effects β i and random effects u i , respectively. The model assumes the random effects ( u ) follow a multivariate normal distribution as u = u where σ 2 g i and σ 2 e i are the genetic variance and the residual variance for trait i, respectively. Then, the genetic correlation ( r g ) was computed as: where σ g 1 ,g 2 is the genetic covariance between two traits.
The genetic correlation estimated by means of the multivariate model could be biased from the phenotypic correlation. Therefore, we further tested correlation between GEBV for each trait using GBLUP model. In the prediction model for HC, SL was included as the covariate to minimize non-genetic effects from SL. The prediction models are described as following: where y is a vector of phenotypes; X is an incidence matrix for the fixed effect β (for the prediction of HC, SL was added as a covariate); Z is an identity matrix for the random effects u , which models the breeding values; ε

Genome-wide association study (GWAS). To investigate the associated markers with transformed HC
and SL, GWAS was performed based on the linear mixed model: where y is the vector of the phenotypes; β is a vector of fixed effects other than SNP effects; g is the vector of random effects that models the polygene background; τ is a vector of fixed effects which represent the additive SNP effects; X , Z , and S are incidence matrices relating to β , g , and τ , respectively. ε is a vector of normal residuals. g and ε follow the normal distributions as g ∼ N 0, K σ 2 g and ε ∼ N 0, Iσ 2 ε , respectively; where K is the realized relationship matrix described above. A restricted maximum likelihood (REML) algorithm was performed to solve the linear mixed model using GWAS function of R/rrBLUP-4.6 with the parameter: n.PC = 10. The p-values were calculated for each SNP marker. The Bonferroni-corrected significant threshold was set as α = 7.454 × 10 -6 (0.05 divided by the number of SNPs: 6707).
Model comparison of genomic prediction. Predictive abilities were obtained under 12 regression models described below. The tenfold cross-validation scheme was applied for the ability calculation following the procedure described by Hosoya et al. 63 . Samples were randomly and equally divided into ten subsets: one for testing and the remaining for training. The phenotypic values of the test set were masked, and the regression model was trained using the training set. GEBVs of the test set were predicted and predictive ability was calculated as the correlation (Pearson's r) between GEBVs and observed values of the test set. This step was repeated with rotating the test sets among the ten subsets, and the average of Pearson's r was obtained. This cross-validation process was repeated 10 times to obtain the mean and the standard error of the mean (S.E.) for the predictive abilities. Transformed HC instead of the original phenotype was used in GBLUP and Bayesian models. To generate the consistent random state for sampling among 12 models, we fixed seeds for random sampling among the models.
Genomic best linear unbiased prediction (GBLUP). The same model as the Eq. (4) was used for GBLUP.
Bayesian models. The models of Bayes A, B, C, Ridge, and LASSO 64,65 are expressed as follows: where y, X, β and ε are same as GBLUP; µ is an intercept; 1 n is a vector of one; p is the total number of genotypes for one individual; z i is a vector of genotypes at SNP i; g i is a vector of random effects that represent the genetic effects for SNP i. following a specific prior distribution. Bayes A assumes a scaled-t distribution for the prior while Bayes B assumes a mixture of gaussian distribution and a point mass at zero. The prior of Bayes C is a mixture of scaled-t distribution and a point of mass at zero. The prior of Bayes Ridge and Bayes LASSO is a normal distribution and a double exponential distribution, respectively. These models were solved using R/ BGLR-1.0.8 66 with nIter = 10,000 and burnIn = 2000.

Bayesian reproducing kernel Hilbert spaces regression (Bayesian RKHS). Bayesian RKHS is a
Bayesian approach of semi-parametric regression 67 structured as: where each parameter is the same as the Bayesian models, while u and ε follow the normal distribution as u ∼ N(0, K σ 2 u ) and ε ∼ N(0, Iσ 2 ε ) , respectively; where K is a reproducing kernel which approximates the marker-based relationship matrix and I is an identity matrix. The model was solved using R/BGLR-1.0.8 with nIter = 10,000 and burnIn = 2000.

Support vector machine regression (SVR).
The SVR method can be viewed as a convex optimization problem that finds a function from genotypes to phenotypes at most ε-deviation from all observed phenotypes while balancing the model complexity and prediction error 68,69 . The method of Lagrange multipliers is used to solve the optimization problem, and the derived approximate function follows: where the input x is a vector of all genotypes for a single sample; N is the sample size; a * i and a i are Lagrange multipliers; k(x, x i ) is a kernel function; x i is a vector of genotypes for individual i ; b is a residual. SVR-linear, SVR-poly, and SVR-rbf using linear, polynomial, and radial basis for kernel function, respectively. The SVR models were implemented by the sklearn.svm.SVR function in Python/scikit-learn-0.20.3 70 . The gamma parameter was set to 'auto' and the default setting was used for the other parameters.
Neural networks. Feedforward neural networks (FNN), inspired by the biological neural network, can model genotype-phenotype regression 71 . Neural cells were modeled by non-linear functions (or activation func-(5) y = Xβ + Zg + Sτ + ε, y = µ1 n + u + ε, Scientific Reports | (2020) 10:19976 | https://doi.org/10.1038/s41598-020-77069-z www.nature.com/scientificreports/ tion) and the network was mimicked by the chain structure. Our FNN had one input layer, two hidden layers, and a regression output. The number of input units was 6707, equivalent to the number of SNPs. The first hidden layer has 200 hidden units and the second 20. The rectified linear unit was used as an activation function in hidden layers. FNN was trained by minimizing the loss function, that is, the mean squared error in this case: where n is the sample size of the training group; y i and y i are observed value and the predicted value of individual i , respectively. Multi-task deep neural network (Multi-task FNN) is based on an assumption where HC and SL share underlying genetic architecture to some extent. These models can improve the accuracy of estimation of the main output using the related auxiliary task as an inductive bias to the main task in reproducing kernel Hilbert space 72,73 . The model has one input layer, two hidden layers, one main regression output, and one auxiliary regression output. The first hidden layer has 200 units, which is a sharing layer for both tasks. Both outputs have separated second hidden layers that differ in the number of the hidden units (20 for the main trait and 100 units for the auxiliary trait). For the estimation of HC, the main regression output is for HC and auxiliary regression output is for SL. The model setting of the main regression and auxiliary regression output was reversed for SL estimation. The activation function and the loss function were the same as the FNN model described above. FNN and multitask FNN were implemented in Python/keras-2.4.3 package 74 with tensorflow-gpu-2.2.1 backend 75 . FNN used "Adam" optimizer, and multi-task FNN used "RMSprop" optimizer, both with the default parameters. Both models were trained by model.fit method in Python/keras with the parameters as epochs = 30, batch_size = 128, and others followed the default. Many combinations of model architecture, loss function, activate the function, and optimizer was arbitrarily chosen and tested to find the models here that have a high accuracy of GP for HC.
Simulation study. To investigate the breeding strategy that can improve SL and HC simultaneously, six scenarios different in recurrent selection schemes were simulated for ten generations with 50 replicates using R/ AlphaSimR-0.11.0 package 76 . The tested scenarios were named for the selection scheme applied: random mating (RAND), GS on HC only (GS HC ), GS on SL only (GS SL ), Smith-Hazel index (S1 SHI and S2 SHI ), and Desired gains index (S DGI ). The workflow of the simulation study is illustrated in Fig. 4 and the details are described in the Supplementary material S2. LGSI forS1 SHI , S2 SHI and S DGI was constructed as: where b is a vector of index coefficients; y is a vector of GEBVs. b for S1 SHI and S2 SHI was computed as: where P and A are phenotypic and genetic variance-covariance matrices, respectively; w is the economic importance of both traits and set as [− 1, 1] assuming both traits have equal economic weights (HC is expected to be decreased by selection) for S1 SHI . P was obtained by varP function of R/AlphaSimR, and A was obtained by mmer function of R/sommer-4.0.1 following the same procedure for the estimation of genetic correlation except that original HC value was used. On the other hand, b for S DGI was computed as: where P and A are the same as S SHI while d is a vector of desired gains. The combination of d can be chosen arbitrarily depending on the breeding goal of the program. In this study, we set d at [− 3, 0.3] for HC and SL so that SL can be improved preferentially while HC can be reduced by 30% after 10 generations. To compare different selection index methods, the vector of economic importance, [− 3, 0.3], which is the same as d in S DGI was assigned to w in S2 SHI .

Ethics statement.
All experiments done in this study were approved by and carried out in accordance with the IACUC (Institutional Animal Care and Use Committee) of the Graduate School of Agricultural and Life Sciences, University of Tokyo (P17-161). All methods were performed in accordance with the IACUC guidelines and regulations.

Data availability
Amplicon sequence reads have been deposited in the DDBJ Sequence Read Archive (DRA Accession: DRA010341). Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.