Genome Wide Association Study and Genomic Selection of Amino Acid Concentrations in Soybean Seeds

Soybean is a major source of protein for human consumption and animal feed. Releasing new cultivars with high nutritional value is one of the major goals in soybean breeding. To achieve this goal, genome-wide association studies of seed amino acid contents were conducted based on 249 soybean accessions from China, US, Japan, and South Korea. The accessions were evaluated for 15 amino acids and genotyped by sequencing. Significant genetic variation was observed for amino acids among the accessions. Among the 231 single nucleotide polymorphisms (SNPs) significantly associated with variations in amino acid contents, fifteen SNPs localized near 14 candidate genes involving in amino acid metabolism. The amino acids were classified into two groups with five in one group and seven amino acids in the other. Correlation coefficients among the amino acids within each group were high and positive, but the correlation coefficients of amino acids between the two groups were negative. Twenty-five SNP markers associated with multiple amino acids can be used to simultaneously improve multi-amino acid concentration in soybean. Genomic selection analysis of amino acid concentration showed that selection efficiency of amino acids based on the markers significantly associated with all 15 amino acids was higher than that based on random markers or markers only associated with individual amino acid. The identified markers could facilitate selection of soybean varieties with improved seed quality.


INTRODUCTION
Soybean [Glycine max (L.) Merr.] is a major source of protein for humans and livestock in the world. For the past several decades, soybean meal has been the leading protein feed source for the animal and poultry production operations because of its high concentration of protein. Poultry and livestock industries use about 68 and 77% of the soybean meal consumed in the European Union and United States, respectively 1,2 . A major function of proteins in nutrition is to supply adequate amounts of required amino acids (Friedman and Brandon, 2001). Thus, genetic improvement of amino acid composition and balance is an important goal in soybean breeding. Developing new molecular markers for marker assisted selection (MAS) and genomic selection (GS) of amino acid composition in soybean will help to achieve this goal.
Quantitative trait loci (QTL) mapping of amino acids have been reported in soybean. Panthee et al. (2006) identified 32 simple sequence repeat (SSR) markers associated with 16 amino acids in soybean seeds based on 101 F6-derived recombinant inbred lines (RIL) from a cross of N87-984-16 × TN93-99. Fallen et al. (2013) reported ten QTLs associated with 17 amino acids and three genomic regions on chromosome 13 (4.89, 21.51, 40.69 cM) controlled multiple amino acids in 282 F5:9 RILs derived from a cross of Essex × Williams 82. As a sole dietary source of protein, soybean is deficient in lysine (Lys), threonine (Thr), methionine (Met), and cysteine (Cys) for poultry and swine. Warrington et al. (2015) conducted QTL analysis for the four amino acids in the Benning × Danbaekkong soybean population with 98 SSRs and 323 single nucleotide polymorphism (SNP) markers, and detected two QTLs on chr 8 and 20 for Lys; three on chr 9, 17, and 20 for Thr; four on Chr 6, 9, 10, and 20 for Met; and one on chr 10 for Cys ( Van Warrington, 2011;Warrington et al., 2015). Khandaker et al. (2015) analyzed MD96-5722" × "Spencer" RIL population and identified 13 QTLs associated with amino acids. However, reports of genetic diversity of amino acids and mapping of QTLs controlling amino acid in soybean germplasm are limited.
Because SSR, SNPs, and indels are abundant in plants and can be assayed with high-throughput technology, the markers have been widely used for genetic linkage mapping, association studies, diversity analysis, and tagging of genes controlling important traits (Liang et al., 2010;Lehne et al., 2011;Li et al., 2014;Shi et al., 2016;Taranto et al., 2016;Zatybekov et al., 2017;Qin et al., 2017a;Qin et al., 2017b;Chang et al., 2018). Genotyping by sequencing (GBS) takes advantage of the next-generation sequencing platforms and utilizes a highly-multiplexed system to assay DNA variants from reduced representation DNA libraries of plant materials (Elshire et al., 2011;Sonah et al., 2013). As a cost-effective technique, GBS has been successfully used in implementing genome wide association study (GWAS), genomic diversity study, genetic linkage analysis, molecular marker discovery and GS in plant breeding programs He et al., 2014;Qin et al., 2016;Shi et al., 2017).
With the decreased genotyping cost and improved statistical methods, GWAS and GS offer new approaches for genetic improvement of complex traits in crop species (Bernardo and Yu, 2007;Li et al., 2013;Morris et al., 2013;Yano et al., 2016;Zhang et al., 2017). GWAS is one of the powerful tools to overcome limitations in traditional QTL mapping (Luo et al., 2019). To date, it has been used to identify molecular markers for a broad range of complex traits in different plant species including Arabidopsis (Angelovici et al., 2017), wheat (Peng et al., 2018), maize (Li et al., 2013;Deng et al., 2017), rice (Huang et al., 2010;Yano et al., 2016), soybean (Fang et al., 2017); sorghum (Morris et al., 2013). In soybean research, GWAS were used in agronomic traits (Zatybekov et al., 2017;Chang et al., 2018), seed quality , seed traits (Xia et al., 2018), phosphorus efficiency (Lü et al., 2018), disease resistance (Qin et al., 2017b;Hanson et al., 2018) etc. As soybean is globally cultivated primarily for its protein and oil, and soybean protein is a complete protein as it contains all the essential amino acid that are required for human health. Numerous studies have reported on the QTL mapping and GWAS for protein Li et al., 2019). GS is to select desired individual within a population based on genomic estimated breeding values (GEBVs) (Hayes et al., 2009), GS has been shown more efficient than the traditional MAS for the improvement of traits controlled by QTL with minor effects (Bernardo and Yu, 2007;Heffner et al., 2009;Shikha et al., 2017;Zhang et al., 2017). GS has been applied to various agronomic traits and disease resistance in maize (Bernardo, 1996;Piepho, 2009;Albrecht et al., 2011;Technow et al., 2013;Shikha et al., 2017), rice (Onogi et al., 2015;Spindel et al., 2015;Duhnen et al., 2017), soybean (Jarquin et al., 2016;Xavier et al., 2016), and wheat Rutkoski et al., 2011;Poland et al., 2012;Battenfield et al., 2016), etc. Previous studies reported the efficiency of GS prediction by cross-validation approach (Dawson et al., 2013;Michel et al., 2016) and suggested that the size of the training population was critical (Xavier et al., 2016). Zhang et al. (2018) conducted GWAS for seed composition, including protein, oil, fatty acids, and amino acids, using 313 diverse soybean germplasm accessions genotyped with a high-density SNP array of the Illumina Infinium SoySNP50K BeadChip (Song et al., 2013). After filtered, a total of 31,850 SNPs with minor allele frequency (MAF) ≥5% were used for GWAS in their analysis and 87 chromosomal regions were identified to be associated with seed composition, explaining 8-89% of genetic variances.
However, little GWAS and no GS for amino acid concentrations in soybean has been reported so far. The main objectives of this study were to (1) evaluate amino acid compositions in soybean germplasm from China, Korea, Japan and U.S. (2) identify SNP markers associated with amino acid concentrations of soybean via GWAS, and (3) explore efficiency of GS for amino acids in soybean breeding. The newly identified markers are anticipated to facilitate MAS and GS of nutritional traits in soybean, and the soybean accessions with high concentrations of amino acids will be potential parents for soybean breeding.

DNA extraction, GBS, and SNP Discovery
Genomic DNA was extracted from freeze-dried fresh leaves of soybean plants using the CTAB (hexadecyltrimethyl ammonium bromide) method (Kisha et al., 1997). DNA library was prepared using the fragment digested by restriction enzyme ApeKI following the GBS protocol described by Elshire et al. (2011) and DNA sequencing was performed using GBS method (Elshire et al., 2011;Sonah et al., 2013). The 90 bp pair-end sequencing was obtained from each soybean genotype at the Institute of Genetics and Developmental Biology, Chinese Academy of Sciences, Beijing, China. The GBS dataset contained 3.26 M short-reads or 283.74 Mbp of sequence for each accession. The short reads were aligned to soybean whole genome sequence (Wm82.a1.v1) 3,4 using SOAPaligner/soap2 and SOAPsnp v. 1.05 was used for SNP calling (Li et al., 2009;Li, 2011).
Approximately a half million SNPs were discovered from the 249 soybean germplasm accessions. SNPs were eliminated if MAF was less than 5%, or missing and ambiguous alleles larger than 15%. After filtering, 23,279 SNPs remained for genetic diversity and association analyses.

Amino Acid Content Determination and Phenotypic Data Analysis
Soybean germplasm was grown at three locations, Shijiazhuang (114°83′E, 38°03′N), Cangzhou (116°7′E, 38°03′N), and Handan (114°48′E, 36°62′N) in Hebei province in a randomized complete block design with three replications in June 2012. Each plot consisted of six rows with a row length of three meters and raw space of 50 cm in all trials. The density was 225,000 plants per ha. The soil at Shijiazhuang was cinnamon. The organic matter, available P and available K concentration were 1.74% 29.9 mg/kg, 94.3 mg/kg, respectively. The soil at Cangzhou was light loamy. The organic matter, available P and available K concentration were 1.0-1.2%, 15 mg/kg, and 100 mg/kg, respectively. The soil at Handan was fluviatile loamy and the organic matter, available P and available K concentration were 1.6%, 19.3 mg/kg, 156.2 mg/kg, respectively. The plots were irrigated once at seed-filling stage. Plants were harvested after 95% leaves had fallen off. Ten plants were randomly chosen from the middle of a plot for seed traits analysis.
A total of 15 amino acids, Ala, Arg, Asp, Glu, Gly, His, Ile, Leu, Lys, Phe, Pro, Ser, Thr, Tyr, and Val in soybean seeds were measured by Biochrom 30 amino acid analyzer (Biochrom Ltd, Cambridge, UK) using the acid hydrolysis method (Davies and Thomas, 1973;Tsugita and Scheffler, 1982). Analysis was carried out by ion exchange chromatography under the experimental conditions recommended for protein hydrolysates. Each sample containing 0.1 g soybean seed powder was acid hydrolyzed with 10 ml of 6 N HCl at 110°C for 22 h in a 15 ml vacuum-sealed glass tube. The top hydrolysate in the tube was filtered into another 50 3 https://www.soybase.org/GlycineBlastPages/archives/Gma1.01.20140304.fasta.zip 4 https://www.soybase.org/GlycineBlastPages/index.php?db_select=Gma1.01 ml tube, and water was added to the tube. A total of 1 ml liquid from the 50 ml tube was transferred to a 1.5 ml tube and dried at 55°C, re-dissolved with 1 ml loading buffer and measured in the analyzer. The amino acid composition was calculated from the standard area obtained from the integrator and expressed as a percentage of the total weight.
Statistical analyses of the 15 amino acids were performed by JMP Genomics 7 (SAS Institute, Cary, NC, USA) 5 (Sall et al., 2012). The mean, range, standard deviation (SD), standard error (SE) and coefficient of variation (CV) were estimated for each amino acid concentration using 'Tabulate'; the distributions of amino acid concentrations were drawn using 'Distribution' in JMP Genomics 7.
Population Structure, Genetic Diversity, and Association Analysis STRUCTURE, a program that uses Bayesian method to analyze multi-loci data in population genetics (Pritchard et al., 2000) 6 , was used to analyze population structure and to create Q-matrix for association analysis. We used the default parameters of STRUCTURE 2.0 software: Admixture Model; Allele Frequencies Correlated; and Compute Probability of the Data (Kaeuffer et al., 2007). The number of subpopulation (K) was assumed to be between 1 and 12. Thus, each K was run 10 times, the Markov Chain Monte Carlo (MCMC) length of the burn-in period was 20,000 and the number of MCMC iterations after the burn-in was 20,000. For each simulated K, the statistical value delta K was calculated using the formula described by Evanno et al. (2005). The optimal K was determined using STRUCTURE HARVESTER 7 (Earl, 2012). After optimal K was determined, a Q-matrix was obtained and used in TASSEL 5 (Bradbury et al., 2007) for association analysis. Each soybean accession was then assigned to a cluster (Q) based on the probability that the genotype belonged to that cluster. The cut-off probability for the assignment to a cluster was 0.5. Based on the optimum K, a bar plot with 'Sort by Q' was obtained to visualize the population structure among the 249 accessions. Genetic diversity was also assessed and the phylogenic tree was drawn using MEGA 6 (Tamura et al., 2013) based on the Maximum Likelihood (ML) tree method  with the following parameters. Test of phylogeny: bootstrap method with No. of Bootstrap replications 500; Model/Method: General Time Reversible model, Rates among Sites: Gamma distributed with Invariant sites (G/I), Number of Discrete Gamma Categories: 6, Gaps/Missing Data Treatment: Use all sites, ML Heuristic Method: Subtree-Pruning-Regrafting-Ex-tensive (SPR level 5), Initial Tree for ML: Make initial tree automatically (Neighbor Joining), and Branch Swap Filter: Moderate. The population structure and the cluster information were imported to MEGA 6 for combined analysis of genetic diversity. For sub-tree of each Q (cluster), the shape of 'Node/ Subtree Marker' and the 'Branch Line' was drawn using the same color scheme of the STRUCTURE analysis.
Association mapping for the 15 amino acids was conducted separately based on the mixed linear model (MLM-Q+K) in TASSEL 5 8 (Bradbury et al., 2007) The SNP markers were considered significantly associated with amino acids if logarithm of the odds (LOD) value ≥3.0 based on MLM-Q+K models.
linkage Disequilibrium Analysis and SNP-Based haplotype Blocks TASSEL 5.0 (Bradbury et al., 2007) was used to calculate the linkage disequilibrium (LD) (r2) for all pairwise loci within a window of 1MB of each chromosome. Haplotype blocks (HAP) were constructed in Haploview (Barrett et al., 2004) with a cutoff of 1% (Contreras-Soto et al., 2017). The LD (r2) for all marker pairs was performed using the R script LDit 9 .

Candidate Gene Selection
Two databases including the annotations for genes at Soybase at https://www.soybase.org/dlpages/ 10 and the plant metabolic network (PMN) database 11 , were used for searching candidate genes related to amino acids in soybean.
Currently, three Williams 82 genome sequence assemblies are available at Soybase (Glyma1.1, and Glyma 2.0) 10 . However, we used Glyma1.1 as the reference because the SNP data were provided by Institute of Genetics and Developmental Biology, Chinese Academy of Sciences, at the time, Glyma1.1 was the best assembly available. We downloaded gene annotation of Glyma1.1 from Soybase and the corresponding gene positions in the Glyma 2.0 were obtained from https://www.soybase. org/correspondence/index.php 12 . For each SNP significantly associated with amino acids, we searched candidate genes within 10 kb of the SNP position. We also downloaded gene annotation from PMN for candidate gene discovery, because the metabolic pathway in PMN is updated with newer version of the genome (Phytozome v12: Gmax_275_Wm82.a2.v1.protein.fa).

Method 1: Ridge Regression Best Linear Unbiased Prediction
Ridge regression best linear unbiased prediction (RR-BLUP) was used to predict genomic estimated breeding value (GEBV) in GS and performed in the rrBLUP package (Endelman, 2011) with the R software Version 3.5.0 (Thuiller et al., 2009). The rr-BLUP is an effective and accurate prediction method as demonstrated in a wide range of traits and crops (Heslot et al., 2012;Jarquín et al., 2014;Lipka et al., 2014;Zhang et al., 2016).
We used 4:1 size ratio of training set and validation set randomly selected from the 249 accessions, which is a fourfold cross-validation, and repeated 100 times. Each training population subset consisting of 199 accessions was randomly selected from the association panel, and the remaining 50 accessions as the validation set (Resende et al., 2012;Shikha et al., 2017).
Two sets of SNPs were used to predict GEBV for each amino acid concentration in each accession: (1) all 23,279 high quality SNPs from GBS, and (2) all 231 SNP markers associated with 15 amino acid concentrations with LOD ≥3.0 from GWAS. In addition, we predicted GEBV for each amino acid concentration based on the SNP markers associated with the amino acid.
The prediction accuracy was estimated using the average Pearson's correlation coefficient (r) between the GEBVs and observed values for each amino acid concentration in the validation set (Zhang et al., 2010;Resende et al., 2012;Shikha et al., 2017). The training and validation sets were randomly created 100 times and the r value was estimated each time. The average r value was calculated for each amino acid. The r value indicates the prediction accuracy and the selection efficiency of GS.

Method 2: Genomic Best Linear Unbiased Prediction
GS was also performed with the genomic best linear unbiased prediction (gBLUP) and the method was extended to compressed best linear unbiased prediction (cBLUP) by using the Compressed Mixed Linear Model (CMLM) approach in GAPIT (Lipka et al., 2012;Tang et al., 2016; http://www.zzlab.net/GAPIT/gapit_help_ document.pdf). In order to conduct a four-fold cross-validation for estimating prediction efficiency, we randomly selected 199 accessions as the training set and the remaining 50 accessions as the validation set to predict GEBV for each accession. GEBV was calculated using the cBLUP in GAPIT using the SNP markers which were associated with the 15 amino acid concentrations with LOD ≥3.0 from GWAS. The Pearson's correlation coefficient (r) between GEBV and observed value of the amino acid concentrations in both training and validation sets were calculated based on the 249 accessions. A total of 100 replications were used to calculate the r values and the average r value for each amino acid was used as the indicator of prediction accuracy.

Phenotypic Variation and Association of Amino Acids in Soybean Seeds
The concentration of 15 amino acids, Ala, Arg, Asp, Glu, Gly, His, Ile, Leu, Lys, Phe, Pro, Ser, Thr, Tyr, and Val varied widely among the 249 accessions (Supplementary Table 2). Concentration distribution of all amino acids except for Val, Ile and Gly in the accessions was near normal, indicating the amino acids are complex traits (Supplementary Figure 1). Glu and Asp were the main components of soybean seeds, which consisted of 20.1% and 13.3% of the total 15 amino acids, respectively. Glu had the highest concentration (74.42 ppm) among the 15 amino acids, followed by Asp (49.15 ppm). Two to five times of difference were observed between the accessions with the lowest and the highest concentration of Arg, Gly, Ile, Leu, Pro, Thr, and Val (Supplementary Table 2). The large variations of the amino acids were also indicated by the high CV values (Supplementary Table 2).
Most of the correlation coefficients among the 15 amino acids were greater than the threshold of 0.124 at P = 0.05 significant level (Table 1). Significant and negative coefficients were also observed between Asp and Ile, Asp and Val, Ile and Gly, Ile and Ser, etc. (Table 1). Based on the correlation coefficient values, the 15 amino acids except for Arg, His, and Pro could be divided into two groups ( Table 1). Group one consisted of five amino acids: Ala, Asp, Glu, Gly and Ser, their pairwise correlation coefficients were greater than 0.75 except for the pair between Glu and Gly (r = 0.6) ( Table 1). Group two contained seven amino acids: Ile, Leu, Lys, Phe, Thr, Tye, and Val with r values greater than 0.48 for all pairs. However, most correlation coefficients of amino acids between the two groups were negative (Table 1). Since the content of amino acids within each group were all significantly and highly correlated, they could have practical application in breeding program, e.g. breeders don't need to improve amino acid individually, they can simultaneously improve multiple amino acids within the same group.

Association Mapping and SNP Marker Identification
The population structure of the 249 soybean accessions was initially inferred using STRUCTURE 2.3.4 (Pritchard et al., 2000) and the peak of delta K was observed at K = 6, indicating the presence of six sub-populations (clusters, Q1-Q6) (Figure 2A).
In total, 51 of the 249 accessions were assigned to Q1 subpopulation with 50 accessions from China; 65 assigned to Q2 with 42 from U.S., 21 from China and two from Korea; 55 assigned to Q3 with 54 cultivars from China; 42 assigned to Q4 with 27 cultivars from China and 12 accessions from U.S.; 21 assigned to Q5 with 16 from U.S.; and 15 to Q6 with all 15 from China ( Figure 2B, and Supplementary Table 1). Phylogenetic analysis of the 249 soybean accessions using MEGA 6 also showed that the clustering of accessions was consistent with that inferred by STRUCTURE ( Figure 2C).
A total of 318 SNP markers consisted of 231 SNPs were associated with the 15 individual amino acid at LOD ≥3 (Supplementary Table 3 and Supplementary Figure 2). Because some SNPs were associated with two or more amino acids as pleiotropic association, the number of SNPs was only 231 ( Table 2). Of the 318 SNPs, 11 were associated with Ala, 29 with Arg, 9 with Asp, 34 with Glu, 29 with Gly, 19 with His, 51 with Ile (Figure 3), 20 with Leu, 14 with Lys, 9 with Phe, 24 with Pro, 11 with Ser, 21 with Thr, 13 with Tyr, and 24 with Val (Supplementary Table 3 and Supplementary Figure 3).
The total number of haplotype blocks was 3,458 based on 23,279 SNPs, the 231 SNPs were positioned in 85 of these haplotype blocks (Supplementary Table 3). Many haplotype blocks contained more than two SNP markers. For example, Gm12_4525341 and Gm12_4525326 were in the same haplotype block and associated with Arg; Gm06_289575, Gm06_399885, and Gm06_582930 were in the same haplotype block on Chr 6 and were associated with Gly (Supplementary Table 3).
The number of the haplotype blocks varied among chromosomes, e.g. 12 of the 85 haplotype blocks were on Chr 16; 11 haplotype blocks on Chr 18; 1 on Chrs 6 and 9. Twenty of the 85 haplotype blocks had significant association with more than one amino acids, e.g. Gm20_42531505 on the Chr. 20_Block 2 was significantly associated with Thr, Gly, Ile, Tyr, Leu, Phe; Two SNP markers, Gm04_43207248 and Gm04_43207187 in the Chr.  Table 4).
Based on phenotypic patterns of the amino acid concentration among accessions, the 15 amino acids could be divided into two groups which were showed in phenotypic variance section. SNP markers associated with amino acids in each group were also found. Twenty-five SNP markers were associated with five amino acids, Ala, Asp, Glu, Gly, and Ser in group one (Table 3), and 28 SNP markers with seven amino acids, Ile, Leu, Lys, Phe, Thr, Tyr, and Val in group two ( Table 4). The SNP markers in each group can be used to simultaneously select multiple amino acids within the group. Such as Gm10_48103776 was associated with five amino acids, Ala, Asp, Glu, Gly, and Ser in group one with LOD values of 2.93, 3.15, 3.51, 2.35, and 3.60, respectively (Table 3) and it can be used to simultaneously select soybean lines with higher contents of the five amino acids in soybean breeding progress. For group two, such as Gm20_42531505 was associated with seven amino acids, Ile, Leu, Lys, Phe, Thr, Tye, and Val with LOD values of 3.53, 4.55, 2.89, 4.79, 5.04, 3.87, and 2.10, respectively ( Table 4), indicating that it can be used to simultaneously select the soybean lines with higher contents of seven amino acids. Meanwhile, both phenotypic and genetic data supported there were two groups of amino acids existed in soybean.

Candidate Gene Selection
The linkage disequilibrium (LD) of soybean genome was analyzed, the average distance of markers at half of the maximum LD decay rate was about 200kb. Considering the LD decay value may vary from genomic region to region, we used the 10kb windows as previously reported (Xie et al., 2018). We identified 704 genes with all or partial sequence within the 10 kb windows that flanked each of the 217 out of 231 unique SNPs associated with one or more amino acids (Supplementary Table 5) and the other 14 SNPs did not have any candidate genes at the 10 kb windows on the chromosomes. Based on gene annotations of the soybean whole genome assembly Gmax_275_Wm82.a2.v1 from Soybase and PMN (Phytozome v12: Gmax_275_Wm82.a2.v1.protein.fa), we found that 15 SNPs were in 14 genes related to amino acid metabolism in gene ontology annotation terms (Supplementary Table 6), e.g. in the region flanking the SNP Gm03_36417795, there was a candidate gene "Glyma03g28476 (Glyma 1.1)/Glyma.03g129100 (Glyma 2.0)" encoding for pyrroline-5-carboxylate reductase (Delauney and Verma, 1990) 13 (Supplementary Table 6). This enzyme catalyzes the last step of L-proline biosynthesis through the L-glutamate degradation pathway. In the region flanking the SNP Gm03_36465287, there was a gene Glyma03g28530 (Glyma 1.1)/Glyma.03g129700 (Glyma 2.0) encoding β L-selenocystathionase, a key enzyme catalyzing L-homocysteine and L-cysteine interconversion. L-homocysteine and L-cysteine interconversion is an intermediate step for conversion between methionine and cysteine (McCluskey et al., 1986) 14 (Supplementary Table 6).  SNP markers associated with 15 amino acid, and SNP markers associated with an individual amino acid.

Genomic Selection for Amino Acid Concentration Based on RR-BlUP in rrBlUP
The correlation coefficients between GEBV and observed value varied among amino acids based on all 23,279 SNPs (column-2 in Table 5), the r value was 0.61 for Arg; 0.50 for Phe; between 0.35 and 0.50 for His, Lys, Thr and Tyr; between 0.25 and 0.35 for Ala, Glu, Ile, Leu, Pro, and Val; and less than 0.25 for Asp, Gly, and Ser. The r values for most amino acids were less than 0.5, suggesting GS prediction accuracy for most amino acids was low based on genome-wide random SNPs.
The correlation coefficients between GEBV and observed value of the 15 amino acids were equal or higher from 231 SNPs than those from the 23,279 SNPs (column-3 vs column-2 in Table 5). The r value was larger than 0.6 for Arg, Ile, Lys, Phe, and Thr, and between 0.5 and 0.6 for Asp, Gly, His, Leu, Tyr, and Val, indicating that associated markers were more efficient to predict amino acids for soybean lines than all the SNPs (Figure 4 and column-3 in Table 5).
Of the 231 SNPs, a total of 171, 42, 12, 4, 1 and 1 SNPs were associated with only one, two, three, four, five, and six amino acids, respectively. A total of 11,29,9,34,29,19,51,20,14,9,24,11,21,13, and 24 SNP markers were associated with Ala, Arg, Asp, Glu, Gly, His, Ile, Leu, Lys, Phe, Pro, Ser, Thr, Tyr, and Val, respectively (Supplementary Table 3). We used the SNP markers only associated with individual amino acid to predict the GEBV for each amino acid, the r values for the 14 amino acids were higher than those from the 23,279 SNPs except for Phe, but equal or lower than those from the 231 SNP markers except for Val (Table 5).
T-test was conducted to compare the r values from the 231 SNPs and from the all 23,279 SNPs and found that the r value from FIGURe 3 | The QQ plot between the expected LOD (-log(P-value)) value and the estimated LOD (log(P-value)) value of amino acid Ile based on 23,279 SNPs as an example (all 15 QQ-plot for the 15 amino acids showed in Supplementary Figure 3). the 231 SNPs in column-3 for each amino acids was significantly higher than that in column-2 from all SNPs with P = 0.01 level in Table 5, indicating that using the associated SNPs had better prediction for GS than using all randomly SNPs (Table 5).

Genomic Selection for Amino Acid Concentration Based on CMlM in GAPIT
Based on cBLUP method using CMLM in GAPIT, the average r was estimated ( Table 5 and Figure 5). The average correlation coefficient in the training set was greater than 0.7 and was higher than those in validation set. The average values in validation set were greater than 0.5 for amino acids except for Pro. Two comparisons were tested to validate the stability of GS using different estimate methods and approaches: (1) RR-BLUP in rrBLUP vs cBLUP in Gapit, and (2) selfvalidation (training set by itself) vs cross-validation (training set). For the first comparison, the 15 r values in column-3 ("231 SNPs in 249 accessions") was compared to those in column-6 ("231 SNPs in validation set") in Table 5 and we found a strong association between the average r values from RR-BLUP in rrBLUP and from cBLUP in Gapit (r = 0.85) based on the 231 associated SNPs. For the second comparison, the 15 r values in column-5 ("231 SNPs in training set") was compared to those in column-6 ("231 SNPs in validation set") in Table 5 and we found a strong association between the average r values from cBLUP in Gapit (r = 0.84) based on the 231 associated SNPs. The strong association with high r value >0.8 between different methods and approach showed that we can use the 231 SNPs to select high amino acid content in soybean through GS.

Application of Marker-Assisted Selection to Genetic Improvement of Soybean
Previous studies using bi-parental segregating populations have identified QTLs controlling 15 amino acids in soybean seeds (Panthee et al., 2006;Fallen et al., 2013;Khandaker et al., 2015;Warrington et al., 2015). The QTL were associated with 84 molecular markers on 14 chromosomes (Supplementary Table 7). In this study, we identified 231 unique SNP markers significantly associated with 15 amino acids (Supplementary Table 3). Eight SNPs were in the same regions of SSR markers that were associated with amino acid concentrations reported by Panthee et al. (2006) (Panthee et al., 2006). Two SNP markers,  Table 7). In addition, Gm09_43488824 at 43,488,824 bp on chr 9 associated with Asp was in the regions reported by Panthee et al. (2006) and Fallen et al. (2013). As GWAS for amino acid concentrations in soybean, Zhang et al. (2018) reported that 54 SNPs, as 92 markers were associated with 18 amino acids; 38 of the 54 SNPs associated with only one amino acid; and 11 SNPs associated with 2 to 12 amino acids. The SNP markers for each amino acid were located at one chromosome such as Pro or Ser, nine chromosomes such as Arg or Asp, up to 11 chromosomes such as Try (Supplementary Table 7). Comparisons with the SNP markers associated with amino acids reported by Zhang et al. (2018), most of SNP markers were located at different regions of soybean chromosomes. However, there were four regions similar to our results: (1) 3.71-3.82 Mb of chr 7 for Arg; (2) 33.85-35.73 Mb of chr 16 for Arg; (3) 16.28-17.65 Mb of chr13 for Asp; and (4) 8.27-9.33 Mb of Chr 8 for Gly. From our study, the SNP marker Gm07_3811476 was associated with Arg at 3,811,476 bp on chr 7, which was near with around 98 kb to the SNP markers ss715597475 at 3,713,267 bp on chr7 for Arg reported by Zhang et al. (2018). Another SNP, Gm16_33853366 close to ss715624781 with 1.87 M distance on chr 16 was also associated with Arg; Gm16_33853366 was at 33,853,366 bp and ss715624781 at 35,721,993 bp on chr 16. For Asp, the Gm13_17646967 at 17,646,967 bp was close to ss715616790 at 16,286,313 bp with a distance 1,360,654 bp on chr 13. The SNP markers, Gm08_8480396 and Gm08_8538031 associated with Gly from this study were close to the two SNP markers, ss715602750 and ss715602851 with Gly  and the four markers are located at a region with one FIGURe 4 | The correlation coefficient (r) among 15 amino acids between the observed values (each amino acid concentration) and the GEBVs predicted from the 231 SNP markers using RR-BLUP in rrBLUP software.
Mb distance on chr 8. Thus, the four regions were validated to be associated with one of the amino acid, Arg, Asp, or Gly. These SNPs identified for 15 amino acids in this study can be used as molecular markers to select lines with high amino acids content through marker-assisted selection (MAS). PCR-based KASP SNP genotyping can be used in soybean breeding program to select high amino acids through MAS. Targeted region sequencing such as tGBS (targeted genotyping-by-sequencing) (Simko et al., 2018) can also be used for MAS and GS based on the sequences flanking these SNPs (Ott et al., 2017).
From this study, 14 candidate genes were found to be related to amino acid metabolism based on gene annotations from Soybase and PMN with gene ontology annotation terms using the DNA sequences in the 15 regions with the 15 SNPs in column-B of Supplementary Table 7 significantly associated with amino acids (Supplementary Table 6). Our further research will develop the molecular markers such as PCR-based assays or targeted region sequencing to validate these candidate genes in our association panel and others. Gene-silence through CRISPR/Cas9 may be used as an approach to validate these candidate genes.

Genomic Selection
Prediction accuracy is the main parameter to measure the performance of GS (Jarquin et al., 2016;Zhang et al., 2016;Duhnen et al., 2017). Prediction accuracy is affected by several factors including GS models, marker density, level of LD, QTL number, the population size specially the training population size, relationship between training population and validation population, and trait heritability (Jarquin et al., 2016). Zhang et al. (2016) estimated prediction accuracy (r value) of seed size based on 309 soybean accessions and reported r = 0.85 when 2000 SNPs or 31,045 SNPs were included, r = 0.8 when 1000 SNPs or 500 SNPs were used. They also identified 48 SNPs on 12 chromosomes associated with seed size based on GWAS. The r value ranged from 0.64 to 0.74 when 5, 10, and 15 of the 48 SNP markers were used, which were 25% higher than those calculated from the same number of randomly selected SNPs. Our results showed that the highest r value (0.56) was obtained based on the model including 231 SNPs significantly associated with one or multiple amino acids, followed by the model including SNPs significantly associated with individual amino acid (r = 0.45), and the least was the model including all SNPs (r = 0.34). A t-test showed r values were significantly different among the sets.
We also estimated the GEBV and r values using the cBLUP in GAPIT. Based on the set of 231 SNPs, the correlation coefficient was greater than 0.7 in the training population and greater than 0.5 in TABle 5 | The averaged correlation coefficient (r) among 15 amino acids between the observed values (each amino acid content) and the GEBVs predicted from (1) all 23,279 SNPs, (2) the 231 SNP markers, and (3) only the associated SNP markers with the specific amino acid content using RR-BLUP in rrBLUP software, and from (4) the 231 SNP markers in reference set (training set) and inference set (validation set) using CBLUPin Gapit.

RR-BlUP in rrBlUP
CBlUP in Gapit *Associated SNPs signifies that the average correlation coefficient (r) for each amino acid in column-4 was calculated with the SNP markers only associated with the individual amino acid to predict the GEBV for each amino acid, such as for r = 0.33 for Ala, which was calculated from 11 associated SNPs.
FIGURe 5 | The average correlation coefficient (r) among 15 amino acids between the observed values (each amino acid concentration) and the GEBVs predicted in both training set and validation set from the 231 SNP markers using cBLUP method in GAPIT.
Frontiers in Plant Science | www.frontiersin.org November 2019 | Volume 10 | Article 1445 validation population. The high correlation between the reference and inference (0.84) based on 15 amino acids, further confirmed the reliability of the GS. A high correlation (0.85) of the prediction accuracy between rrBLUP and GAPIT based on 231 SNPs, indicated that both RR-BLUP in rrBLUP or cBLUP in GAPIT were consistent.

CONClUSION
In this study, soybean accessions with high concentrations of amino acids in seeds, and molecular markers associated with individual and groups of amino acids were identified. These soybean accessions with high amino acid concentrations could be used as parents in soybean breeding programs. The SNP markers strongly associated with the concentrations of the amino acids could be used to improve the nutritional quality of soybean through marker-assisted selection. In addition, fourteen candidate genes that were related to amino acid metabolism were also identified. These candidate genes will lead to a better understanding of the molecular mechanisms that control amino acids metabolism in soybean seeds. Genomic selection analysis of amino acid concentration showed that the selection efficiency of amino acids based on the markers significantly associated with 15 amino acids was higher than that based on genome-wide random markers or markers only associated with an individual amino acid. These results suggest that including a set of markers significantly associated with multiple amino acids in genomic selection is likely to help breeders to efficiently select soybean varieties with improved amino acid content.

DATA AVAIlABIlITY STATeMeNT
SNP data can be found in the ENA using accession number PRJEB34546 (https://www.ebi.ac.uk/ena/data/view/PRJEB34546).

eThICS STATeMeNT
All data and materials are not related to human and animals. This research is not related to any plant specimens to be deposited as vouchers or any other association for this section.

AUThOR CONTRIBUTIONS
JQ, AS, YC, FW, and WR carried out phenotyping and genotyping. AS, JQ, SL, and QS analyzed the data. JQ composed the draft of the manuscript. MZ and CY directed and managed this research. AS and QJS reviewed and edited the manuscript. All authors have read, made corrections, and approved the final manuscript.

FUNDING
The authors would like to thank Prof. Lijuan Qiu (Chinese Academy of Agricultural Sciences) for providing seeds of 249 soybean accessions. This study was supported by the Hebei Province Natural Science Foundation for Distinguished Young Scholars (C2014301035), National Natural Science Foundation of China (31100880), and Key Project of the Natural Science Foundation of Hebei Province (C2012301020).