Patterns of Genetic Variation in a Soybean Germplasm Collection as Characterized with Genotyping-by-Sequencing

Genomic characterization is playing an increasing role in plant germplasm conservation and utilization, as it can provide higher resolution with genome-wide SNP markers than before to identify and analyze genetic variation. A genotyping-by-sequencing technique was applied to genotype 541 soybean accessions conserved at Plant Gene Resources of Canada and 30 soybean cultivars and breeding lines developed by the Ottawa soybean breeding program of Agriculture and Agri-Food Canada. The sequencing generated an average of 952,074 raw sequence reads per sample. SNP calling identified 43,891 SNPs across 20 soybean chromosomes and 69 scaffolds with variable levels of missing values. Based on 19,898 SNPs with up to 50% missing values, three distinct genetic groups were found in the assayed samples. These groups were a mixture of the samples that originated from different countries and the samples of known maturity groups. The samples that originated from Canada were clustered into all three distinct groups, but 30 Ottawa breeding lines fell into two groups only. Based on the average pairwise dissimilarity estimates, 40 samples with the most genetic distinctness were identified from three genetic groups with diverse sample origin and known maturity. Additionally, 40 samples with the highest genetic redundancy were detected and they consisted of different sample origins and maturity groups, largely from one genetic group. Moreover, some genetically duplicated samples were identified, but the overall level of genetic duplication was relatively low in the collection. These findings are useful for soybean germplasm management and utilization.


Introduction
Genomic characterization of ex situ plant germplasm is playing an increasingly important role in germplasm management and utilization [1][2][3][4]. Characterizing germplasm using genomic tools has been technically and practically more feasible than before, thanks to the recent development of next generation sequencing technologies and bioinformatics tools, and drastically decreasing costs [5]. Such characterization can provide high resolution with genome-wide SNP markers to analyze genetic diversity and structure (e.g., see Milner et al. [3]; Sansaloni et al. [4]), group germplasm for the development of core collections (e.g., Jeong et al. [6]), and can support identification of genetic duplicates (e.g., see Ellis et al. [7]; Singh et al. [8]) for better germplasm management. It can also enhance the search for unique germplasm with traits of breeding targets for better varietal development (e.g., see Sansaloni et al. [4]; Mascher et al. [9]). Considerable efforts were made using genomic tools to characterize ex situ soybean [Glycine spp.] germplasm (e.g., see Jeong et al. [6]; Song et al. [10]) and to develop different core subsets of soybean accessions in different genebanks (e.g., see Wang et al. [11]; Cho et al. [12]; Oliveira et al. [13]; Kaga et al. [14]; Priolli et al. [15]). These characterizations not only allow for a better understanding of the genetic variation present in soybean germplasm, but also support enhanced soybean germplasm management and utilization. For example, Song et al. [10] reported considerable redundant germplasm in the USDA soybean collection with 23% of G. max accessions and 30% of G. soja accessions being at least 99.9% identical, based on 42,509 SNPs, to another accession in the collection. Moellers et al. [16] demonstrated the effectiveness of the USDA soybean core collection in the identification of multiple soybean Sclerotinia sclerotiorum resistances. However, some challenges remain in the genomic characterization of all 176,000 soybean accessions that are currently conserved in more than 24 genebanks worldwide [17,18]. Little is known about the extent of genetic diversity and unique accessions in each collection.
Plant Gene Resources of Canada (PGRC; the Canadian national seed genebank at Saskatoon) maintains a soybean germplasm collection of 1031 accessions. These accessions were acquired largely from Canadian soybean breeding programs over the last 50 years and from the USDA-ARS soybean collection and N.I. Vavilov All-Russian Institute of Plant Genetic Resources for accessions with known early maturity over the last 15 years. Little is known about the genetic diversity and structure of the PGRC collection. Thus, a genomic characterization of the collection was initiated in 2017 with the objective of collecting information to enhance soybean germplasm management and utilization. To make the characterization more relevant to soybean germplasm utilization, we also included 30 soybean cultivars and breeding lines developed and released from the Ottawa soybean breeding program of Agriculture and Agri-Food Canada (AAFC). The specific objectives of the characterization were to: (1) apply a genotyping-by-sequencing technique to genotype 541 PGRC soybean accessions and 30 Ottawa breeding lines; (2) analyze the patterns of genetic diversity in the assayed samples; and (3) identify accessions with the most genetic distinctness and redundancy.

Soybean Germplasm and DNA Extraction
The soybean collection conserved at PGRC has 1031 accessions, but only 743 accessions were available for seed distribution in 2016. For this study, we selected 541 accessions from the 743 accessions representing 28 countries of origin and one group of unknown origin, based on the log proportion of the group size with respect to the 29 groups involved. Fortunately, the selection also included 269 accessions with three known maturity groups (MG) (38 for MG 0, 172 for MG 00 and 59 for MG 000) [19]. Note that soybean cultivars are commonly classified into 13 MGs, depending on their adaptation to photoperiod and seasonal temperature. In Canada, the earliest soybean cultivars are MG 000, which are adapted to a region north of 50 • N. MG 00 and MG 0 cultivars are adapted to Manitoba, Eastern and Central Ontario and Western Quebec, while MG I and later are adapted to regions in Southern Ontario and Southern Quebec. To enhance our analysis of early maturity germplasm, we also acquired 30 short-season soybean accessions representing released cultivars and breeding lines from the active AAFC soybean breeding program at Ottawa Research and Development Centre over the last 30 years.
Approximately five seeds were randomly selected from each accession and planted in the greenhouse at the Saskatoon Research and Development Centre. Leaf tissue was collected at the 2-4 leaf stage separately for each plant, freeze-dried for 1-2 days in a Labconco (Kansas City, MO, USA) Freeze Dryer System, and stored at −20 • C. DNA was extracted from 12-15 mg of freeze-dried leaf tissue from one randomly selected plant per accession using the Qiagen BioSprint 96 DNA Plant Kit (Qiagen Inc., Toronto, ON, Canada), following the product handbook, with the exception of DNA elution in 100 µL rather than 200 µL of water. DNA quality was assessed using a 260/280-nm ratio from the Thermo Scientific Nanodrop 8000, and DNA was quantified using the Invitrogen Quant-iT TM PicoGreen ® dsDNA Assay Kit (Life Technologies, Burlington, ON, Canada) and adjusted to 20 ng/µL with water.

Genotyping-by-Sequencing
We screened genome-wide genetic variability of soybean germplasm using the genetic diversity-focused GBS (gd-GBS) protocol described by Peterson et al. [20]. This protocol modified the original GBS method [21] in combination with the ddRAD procedure [22] to employ long-read Illumina sequencing and indexing for sample deconvolution. Six sequencing libraries, each having 95 samples and one technical replicate (i.e., one sample was present in all six libraries), were prepared following the methods described in Peterson et al. [20] with the following modifications. First, four new adapters with balanced (A/C, G/T) barcodes (Table S1) were applied evenly across the samples to increase sample diversity and avoid phasing read errors on the Illumina HiSeq. Second, ligated fragments were cleaned using only 1× Agencourt AMPure XP Beads (Beckman Coulter, Mississauga, ON) rather than 1.5×, and eluted in 40 µL rather than 30 µL reagent grade water. After preparation, libraries were quantified using a ddPCR TM Library Quantification Kit for Illumina TruSeq (Bio-Rad Laboratories, Inc., Hercules, CA, USA) and diluted to 6 pM. Libraries were sequenced by NRC (National Research Council, Saskatoon, SK, Canada) on a HiSeq 2500 using Rapid Mode, 1 × 250 bp (dual indexed). The sequencing generated 1152 sequence read files in FASTQ format and these raw sequences were deposited into the SRA database of the National Centre for Biotechnology Information (NCBI) under BioProject ID of PRJNA705793.

Bioinformatics Analysis
All FASTQ files were first assessed with FastQC [23] to determine if there was any 3 adapter read-through; if so, such files were trimmed with Trimmomatic version 0.32 [24] using the following trim settings: ILLUMINACLIP:TruSeq3-PE-2.fa; SLIDING-WINDOW:10:24; and MINLEN:100 variables to remove any adapter sequence bases, trim bases where the average Phred score over a 10 base window was less than 24 and remove any sequences shorter than 100 bases, respectively. The statistics of raw and filtered sequence reads were generated using an in-house Perl script.
SNP calling was made for each sample from the trimmed FASTQ files against the soybean genome assembly of Wm82.a2.v1 [25,26]. The reference FASTA chromosome and scaffold names were modified with a custom Perl script to reflect the formatting required for our custom pipeline. Specifically, clean sequences of each sample were aligned to the reference using Bowtie2 version 2.2.6 [27] with the output in SAM format. Variant calling was done using SAMTOOLS version 0.1.18 [28] as follows: converted to BAM format with VIEW; sorted with SORT; variant identification with MPILEUP and converted to VCF format with BCFTOOLS VIEW. VCF file processing and final data output were generated using a collection of in-house Perl and Bash shell scripts referred to as the "referenceGeno" pipeline, which is available upon request to the corresponding author. These scripts helped to (1) filter the resulting SAMTOOLS VCF files to generate a list of potential SNPs in haplotype and genotype with 15%, 30%, and 50% missing data across all 576 accessions and (2) format them in different tab delimited tabular forms for different sequence analyses such as STRUCTURE and MEGA for further genetic analyses.

Genetic Diversity Analysis
SNP genotype data with up to 50% missing values were firstly cleaned with six control samples (one original sample plus its replicates in six HiSeq runs); only those markers having the same marker values for six control samples were kept. There was further removal of the markers with singletons, with the same genotype across the assayed samples and with more than three genotypes per sample. Minor allele frequency, the extent of missing SNPs, SNP distribution with respect to chromosome, country of origin and maturity were generated in a Microsoft Excel ® file.
The genetic structure of the 571 samples was analyzed based on the SNP data of 19,898 markers (or M50 dataset) using a model-based Bayesian method available in the program STRUCTURE version 2.2.3 [29]. The STRUCTURE program was run 30 times with 30-core parallel computing in a Linux server for each subgroup (K) value, ranging from 1 to 10, using the admixture model with 10,000 replicates for burn-in and 10,000 replicates during analysis. The final sample grouping was determined based on (1) the likelihood plot of these models, (2) the change in the second derivative (∆K) of the relationship between K and the log-likelihood [30], and (3) the stability of grouping patterns across 30 runs. For a given K with 30 runs, the run with the highest likelihood value was selected to assign the posterior membership coefficients to each sample. The posterior membership coefficients were displayed in a bar plot. The size and composition of each optimal cluster were analyzed with respect to sample origin and maturity.
The inferred genetic structure was further compared for consistency with the genetic relationships of individual samples obtained from two commonly applied approaches. A principal coordinate analysis (PCoA) of all 571 samples was performed using the R routine, AveDissR, for assessing genetic distinctness and redundancy [31,32] and plots of the first three resulting principal components were generated to assess the sample associations. A neighbor-joining (NJ) analysis of all 571 samples was also conducted using PAUP* [33] based on the SNP data of 19,898 markers and a radiation tree was displayed using MEGA 5.05 [34]. The resulting PCoA plots and NJ trees were individually labeled with respect to sample origin and maturity group.
An analysis of molecular variance (AMOVA) was performed with Arlequin version 3.01 [35] on 19,898 markers to quantify the genetic variation present among 24 groups of sample origin and among five maturity groups ( Table 1). The corresponding pairwise genetic distances among the origin and maturity groups were also generated. Additional AMOVA was also made to quantify the genetic variation among genetic groups identified from the STRUCTURE and PCoA analyses.

Identifying Distinct and Redundant Germplasm
We applied the same approach as proposed by Fu [36] that was equivalent to the simple matching coefficient of Sokal and Michener [37] to calculate the pairwise genotypic dissimilarity and to generate the average pairwise dissimilarity (APD) per sample for the assessment of germplasm distinctness and redundancy. The higher the APD value, the more genetically distinct the sample is in the collection. The lower the APD value, the more genetically redundant the sample is in the collection. Ranking the APD values of all the assayed samples provides a means of identifying the most distinct and most redundant samples [36]. Specifically, we applied AveDissR to generate the APD value per sample. To verify the redundant and distinct samples identified, the PCoA plot was also made with the labels of the identified groups to determine their genetic associations with the whole assayed samples.
Extra effort was also made to generate a pairwise dissimilarity matrix for all 571 samples based on 19,898 markers using AveDissR with modification to output the matrix file. This analysis helped identify the sample pairs with the lowest pairwise dissimilarity values for the assessment of genetically duplicated samples.

SNP Discovery
The HiSeq run of 576 soybean samples (Table S2) yielded approximately 548.4 million raw forward sequence reads from six libraries (Table S3). The number of raw forward sequence reads per sample ranged from 269,633 to 1,618,902 with an average of 952,074 and the filtered read counts varied from 255,484 to 1,521,313 with an average of 883,226 (Table S3). The SNP call detected 43,891 SNPs across the 576 samples on 20 chromosomes (43,161 SNPs) and 69 scaffolds (730 SNPs). As expected for the genotyping-by-sequencing technique, there were a large number of SNPs (or 19,318) having missing values for more than 50% assayed samples. Further removal of the same SNP genotype (2979 SNPs) or more than three genotypes (398 SNPs) across the assayed samples and inconsistent genotypes across the six repeat samples (1298 SNPs) generated three datasets of 19,898 SNPs with missing levels of 50% or less (M50 for short); 13,948 SNPs with 30% or less (M30); and 6861 SNPs with 15% or less (M15). The SNP distributions across the 20 chromosomes and 69 scaffolds for these three SNP datasets are shown in Figure 1A. The patterns of variation in SNP count per chromosome were similar for three missing value levels. On average, each chromosome had 985 SNPs at the M50 level, 691 SNPs at the M30 level, and 341 SNPs at the M15 level. Specific SNP distribution with respect to the level of missing value is illustrated in Figure 1B, and SNP counts generally decreased toward the increased level of missing value (up to 50%). Further characterization of the SNPs in the M50 dataset revealed that a majority of the SNPs had minor allelic frequency of 0.1 or lower ( Figure 1C).

Patterns of Genetic Variability
The Bayesian inference of genetic structure by STRUCTURE without consideration of accession information revealed three optimal clusters of the 571 soybean samples with a strong delta K support (Figure 2). Clusters 1, 2, and 3 consisted of 144, 285, and 142 samples, respectively, and each cluster had mixed memberships originating from different countries (see Table S2 for the detailed memberships). The PCoA plot also revealed three distinct genetic groups of the 571 soybean samples ( Figure 3). The membership of each group was the same as those inferred from STRUC-TURE analysis, with one exception (see Table S2 for the sample membership). One sample in PCoA Group 3 had a membership coefficient of 0.504 for STRUCTURE Cluster 2 (vs. 0.496 for Cluster 3). Also, it was clear that the samples from Canada, China, the Republic of Korea, Japan and Russia were widely spread into three groups (Figure 3, Table 1 and Table S2). The 30 Ottawa breeding lines were located in Group 1 and Group 3 only. The first and second PCoA components explained 2.59% and 1.83% variances, respectively.  The NJ clustering revealed the same three distinct groups of all 571 soybean samples ( Figure 4) as those inferred from STRUCTURE ( Figure 2) and PCoA (Figure 3). However, the NJ tree revealed more sub-groups in each group and the mixture of sample origins in each group. There were no groups specific to samples originating from a specific or single country only ( Figure 4A). Interestingly, all the samples of known maturity groups were also widely spread into three NJ groups, although the 30 Ottawa breeding lines (largely belonging to these three maturity groups) were located only in two NJ groups ( Figure 4B), the same as in the PCoA plot ( Figure 3). The NJ tree also helped to identify six genetically unique samples (CN32352, CN107502, CN45107, CN107642, CN36139, and CN107548). Quantifying the genetic variations of all 571 soybean samples representing 24 groups or countries of origin through AMOVA revealed 7.45% variance residing among these 24 groups and 92.55% variance present within groups ( Table 2). The pairwise genetic distances among 24 groups ranged from 0 to 0.403 with a mean of 0.107, but only 109 out of 276 group pairs showed statistically significant distances (Table S4). This result indicates that the country of sample origin may not always be informative to identify genetically distinct samples. The AMOVA analysis also revealed that 5.55% variance was found to reside among five maturity groups (including Ottawa breeding lines and an unknown origin group) and 94.45% within the maturity groups. The pairwise group distances were not statistically significant among three maturity groups, but these three groups had significant pairwise group distances with the Ottawa breeding lines (e.g., 0.235 between group 000 and Ottawa breeding lines). Further AMOVA analysis of three PCoA groups revealed 34.52% variance resided among the inferred three groups and 65.48% was present within the groups. The largest significant pairwise distance was 0.413 between Group 1 and Group 3, followed by those between Groups 1 and 2 (0.363) and between Groups 2 and 3 (0.282).

Genetic Distinctness and Redundancy
The APD estimates for the 571 assayed samples ranged from 0.0986 to 0.2522 with a mean of 0.1284 and standard deviation of 0.018 (Table S2). The frequency distribution of these APD values is given in Figure S1A, and there were 2, 19 and 106 samples with APD estimates greater than 3, 2 and 1 standard deviation, respectively, while there were only 72 samples with APD estimates lower than one standard deviation. The regression analyses of the APD estimates over the three levels of SNP missing values (50%, 30% and 15%) revealed that the correlations among those APD estimates were high, ranging from 0.813 to 1.027 with R2 greater than 0.902 ( Figure S1).
As the APD value of a particular sample measures the overall genetic difference of the sample against the remaining samples of the collection, we selected 40 more genetically distinct samples with APD estimates greater than 0.1565 and 40 more genetically redundant samples with APD estimates lower than 0.1076 for further germplasm management ( Table 3). Note that the 40 samples with the most genetic distinctness included the first three distinct samples identified from the NJ tree ( Figure 4). The genetic relationships of the selected 80 samples with the remaining samples are also displayed in Figure 3B. The selected genetically distinct samples spread over into three groups, but were mainly located in Group 1, while the selected genetically redundant samples mainly harbored in Group 2, with two samples placed in Group 3. The frequency distributions of these 80 selected samples with respect to origin and maturity are shown in Table 1.
We also identified 16 sample pairs with the lowest pairwise dissimilarity values that were roughly equivalent to 37 or fewer (out of 19,898) SNP markers showing genetic differences between two samples (Table 4). Such extremely low pairwise dissimilarity values clearly indicate the presence of genetic redundancy. The 16 (out of 162,735) sample pairs had only 22 samples, of which 13 were considered as genetically duplicated samples (Table 4). More specifically, each of the 13 samples was at least 99.8% identical, based on 19,898 SNP markers, to at least one other sample. Note that origin of country is given in the ISO 639-2 code; UNK is for unknown and OBL is for the AAFC Ottawa breeding line. CG stands for country or group (see Table 1). MG is for five maturity groups and PG for three PCoA-based groups (see Table S2 or Figure 3). Note that any three of the four samples in each group of sample pairs and any one sample in individual sample pairs could be considered as genetically duplicated samples to at least one other sample.

Discussion
The germplasm characterization presented here revealed some interesting patterns of genetic variation in the soybean collection held at Plant Gene Resources of Canada. First, there were three distinct genetic groups present in the assayed samples, and each group was a mixture of the samples originating from different countries and with different known maturity groups (Table 1 and Figure 3). Second, the samples originating from Canada were clustered in all three genetic groups, but the 30 Ottawa breeding lines were clustered only in two groups (Table 1 and Figure 4). Third, the 40 samples with the most genetic distinctness represented three genetic groups with diverse sample origin and known maturity ( Table  3). In contrast, the 40 samples with the highest genetic redundancy consisted of different sample origins and maturity groups, but were largely from one genetic group. Fourth, the extent of genetic duplication was relatively low in the collection and only 13 samples were identified as genetically duplicated samples (Table 4). These findings are useful for soybean germplasm management and utilization.
Three different diversity analyses (Figures 2-4) revealed the presence of three distinct genetic groups of soybean germplasm in the PGRC soybean collection, but each of these groups was well mixed with diverse sample origins and not unique to germplasm from China, Korea, Japan or Russia. This was unexpected, as soybean is thought to have been domesticated in China around the eleventh century BCE and then disseminated to surrounding countries around the first century CE [38]. Also, it is not consistent with those distinct clusters reported by Song et al. [10] of soybean wild and landrace genotypes from different countries, in which genetic clusters were well aligned with the sample origins (see Figure 1 of Song et al. [10]). One possible explanation is that the inferred genetic groups reflects the unique gene pool generated by the Canadian soybean breeding programs over the last 80 years [39,40], as the PGRC soybean collection was largely acquired from those Canadian breeding programs aimed for improving short-season soybean. Also, the known maturity groups were well presented in the inferred genetic groups (Table 3). However, the diversity analyses of the Canadian soybean cultivars released over the last 80 years did not reveal marked genetic groups [36,41]. Clearly, more research is needed to understand the origin of these distinct genetic groups.
We applied pairwise dissimilarity to identify genetic duplicates and average pairwise dissimilarity to identify the samples with the most genetic redundancy. Overall, the extent of genetic duplication in the PGRC soybean collection was relatively low (Tables 3 and 4). This finding, however, is inconsistent with those reported for the USDA soybean collection [10], but it is not surprising either, given the collection history and acquisition sources. However, such duplication identification was not without limitations [36]. It depends on the number of samples assayed and the SNP markers used. More samples would enhance the power of identification. Fewer genome-wide SNP markers with missing values could also affect the reliability of estimating average pairwise dissimilarity (see Figure S1). Also, the most genetically redundant samples were relative to the whole set of samples assayed, and there is no solid threshold to group samples with and without genetic redundancy. Moreover, we did not examine within-accession variation, as such variation was expected to be low [42] for soybean with an outcrossing rate of 1-2% [43,44]. In spite of this, we cannot rule out the effects, if any, of the within-accession variation on the estimation of pairwise dissimilarity and duplication identification. Thus, these identifications may not necessarily identify the true duplicated accessions per se [36], but suggest the potential duplicate candidates for germplasm management.
Our characterization also generated a unique set of genomic resources for genetic analyses of soybean germplasm. The acquired SNP genotype data had an adequate genome sampling across the 20 chromosomes ( Figure 1A) and can be applied through a genomewide association analysis to identify genetic regions associated with various traits of breeding targets such as agronomic traits, early maturity, and quality traits, if those assayed accessions are phenotypically evaluated. These efforts will facilitate the search for genetic variants of breeding importance from the conserved germplasm. Also, genotyping-by-sequencing is known to generate SNP genotypes with large amounts of missing data [45], as variation could occur in the restriction site during library preparation, PCR bias during library amplification, and/or flow-cell sequencing bias due to increased multiplexing. However, our preliminary analysis seemed to suggest that the SNP genotypes with up to 50% level of missing values across the samples were still informative, at least for estimating genetic diversity ( Figure S1). Further diversity analyses based on the M30 and M15 datasets revealed the same three distinct genetic groups, although the group distinctness was slightly reduced (or more dispersed) for the M15 dataset (results not shown). The patterns of variation with respect to sample origin and maturity group remained essentially the same as for the M50 dataset.
The revealed patterns of genetic variation have implications for managing soybean germplasm. The genetic variation was relatively low with respect to sample origin and maturity group, but considerably larger among three distinct genetic groups (Table 2; Figure 3). Thus, some attention should be paid to manage these three distinct genetic groups. The genetic distinctness measured with average pairwise dissimilarity can be informative to the development of core subsets of soybean germplasm from the collection for germplasm utilization, as the identified samples can serve as the candidates for further consideration, including the field evaluation on the traits of breeding targets. This could be achieved following the same integrated approach for the development of the flax core collection [46]. Also, the revealed genetic distinctness can be useful for the selection of a set of distinct soybean germplasm for safety backup in other genebanks. As mentioned above, the genetic redundancy in the collection was relatively low, as the lowest average pairwise dissimilarity values had just one standard deviation apart from the mean, and those redundancies mainly resided in the genetic group 2. The identified genetic duplicates can be further verified in the field to determine their accession duplication.
The findings presented here also have some implications for germplasm utilization in soybean breeding. Breeders searching for germplasm with early maturity need to screen all the soybean accessions in the collection, as the inferred genetic groups were not unique to specific maturity groups, nor to specific sample origins. A field evaluation of maturity traits is needed. Identification and exploration of genes associated with early maturity traits can proceed with any genetic group, but would be more fruitful with the group in which the majority of the 30 Ottawa breeding lines resided.

Conclusions
The genotyping-by-sequencing generated abundant SNP markers across the 20 soybean chromosomes. The diversity analysis revealed three distinct genetic groups present in the soybean collection held in Plant Gene Resources of Canada. These groups were a mixture of the samples originating from different countries and samples of known maturity groups. The 30 Ottawa breeding lines were clustered with two of the three groups. The analysis also identified 40 samples with the most genetic distinctness and 40 samples with the highest genetic redundancy and showed that the extent of genetic duplication in the collection was relatively low. These findings are useful for soybean germplasm management and utilization.