The diversity panel consisted of 200 genotypes of modern Canadian (CD) cultivars (n = 59), modern Chinese (CH) cultivars (n = 53), and Canadian × Chinese (CD-CH) progeny lines (n = 88) belonging to maturity groups 0, 1, and 2 (Supplementary Table S1). The diversity panel was evaluated in yield trials at the Elora Research Station (43°64′104.4” N; 80°40′567.4′′ W), Elora ON, and Woodstock Research Station (43°08′44.8′′ N 80°47′02.5′′ W), Woodstock ON during 2019 and 2020 field seasons. Two replications were evaluated per environment in a nearest neighbor Randomized Complete Block Design (nn-RCBD) with soybean lines randomly assigned.

Seedling emergence score was recorded for each plot 3 weeks after planting, based on the plot-wise number of plants observed. A scale of 0–10 was used where 0 corresponded to no emergence and 10 corresponding to 100% emergence. Pubescence color, flower color, and leaf morphology were recorded subsequently. Flower color was recorded at the R1 stage (one flower at any node), full maturity date was recorded at R8 stage where 95–100% of pods have turned brown; lodging: scored at maturity on a scale ranging from 1 to 5, where 1 = plants fully upright and 5 = plants fully prostrate; and height: as the distance between the terminal node and the ground, measured in cm (Ernpig and Fehr, 1971). All field observations were recorded on an iPad and exported as an MS Excel file.

Seed quality traits were measured using a Perten Diode Array 7,250 Near Infra-Red Spectroscopy (Springfield, United States) machine following manufacturer’s guidelines. Seeds were screened to remove off-types, dirt, and other impurities. A 100 seed weight was measured with a regular commercial scale. Hilum color, 100 seed weight, plot number, entry numbers, and experiment number were entered into the NIR machine for each entry. NIR results were exported as an excel file and screened for errors. Randomly selected genotypes were re-run to ensure that the readings were consistent. Within the soybean seed quality traits, only the protein and oil concentration (expressed as % on a dry seed basis) and 100 seed weight (g) were retained for analysis. One entry each from Elora 2019 and Woodstock 2019 was removed from analysis due to machine error (Supplementary Table S2).

Analysis of variance of seed yield, seed quality, and agronomic traits was conducted using the PROC GLIMMIX procedure in Statistical Analysis Systems (SAS) version 9.4 (SAS Institute Inc., Cary, NC, United States) for RCBD. The GLIMMIX procedure allows the use of generalized mixed linear model—a standard in agricultural research (Camp et al., 2018). “Genotype,” “environment,” and “genotype-by-environment” were considered fixed effects and “block (environment)” was considered random effect.

Due to the unbalanced number of genotypes in 2019 at Elora and Woodstock, data were sorted by year and environment and separated into three sets: 2019 (with 147 genotypes), 2020 (200 genotypes), and combined years (147 genotypes). Using PROC GLIMMIX procedure, the least squared means (LSMEANS) values were calculated for seed yield, protein concentration, and oil concentration for both combined environments and individual environments. Shapiro–Wilk test was conducted using PROC UNIVARIATE to determine the distribution of residuals. PROC PLOT was used to examine the normality of residual distribution. Homogeneity of error variance was tested by conducting Levene’s test on the absolute residuals.

Comparisons were made between environments, genotypes, and genotypes by environments. Tukey–Kramer multiple comparison test was invoked along with the LINES statement to generate statistically significant differences between comparison groups. CONTRAST statements were used along with ESTIMATE statements to test the statistical differences, if any, between the three different genotypic groups.

Leaf tissue was collected into labeled 10 ml plant-tissue collection tubes. One to two young leaves were collected into each tube. These tubes were then transported on ice back to the Soybean Research Laboratory at the University of Guelph in Guelph, ON. Leaf tissue samples were freeze dried with a Labonco FreeZone^® freeze dry system (Savant Moduly, Kansas City, MO, United States) for a period of 24 h and stored at −4°C.

Genomic DNA was extracted from samples of freeze-dried leaf tissue by using NucleoSpin^® Plant II DNA extraction kit by Macherey-Nagel following the manufacturer’s guidelines. Extracted DNA samples were spot tested with a NanoDrop 8,000 machine (Thermo Fisher Scientific, Waltham, MA, United States) to check for protein/RNA contamination and to verify the quality of genomic DNA. DNA concentration was established with the QuBit 4 DNA Analyzer (Thermo Fisher Scientific, Waltham, MA, United States) and was standardized to 10 ng/μl. A precise volume of 10 μl was pipetted out to two 96-well semi-skirted PCR plates, which were sent to Plateforme d’analyses génomiques [Institut de Biologie Intégrative et des Systèmes (IBIS)], Université Laval (Quebec, QC, Canada) for Genotyping-by-Sequencing (GBS) and SNP calling.

GBS was conducted following the methods and recommendations outlined by Elshire et al. (2011), Sonah et al. (2013), and Torkamaneh et al. (2020a,c). The GBS library was created with ApeKI restriction enzyme digestion. A 158 million single-end reads were generated with an Ion Torrent Proton System (Thermo Fisher Scientific Inc., USA). These were processed using the Fast-GBS.v2 pipeline (Torkamaneh et al., 2020c). FASTQ files were demultiplexed, trimmed, and then mapped against the soybean reference genome (Williams82 (Gmax_275_Wm82.a2.v1); Schmutz et al., 2010) with an average success rate of 94.4%. SNPs were identified from the mapped reads and filtered out if (i) they were multi-allelic, (ii) the overall read quality (QUAL) score was <20, (iii) the mapping quality (MQ) score was <30, (iv) read depth was <2, and (v) missing data >80%. Missing data imputation was performed using BEAGLE v5.1 (Browning et al., 2018) following the protocol laid out by Torkamaneh and Belzile (2015).

GWAS was conducted using the rMVP package in R (Yin et al., 2021) utilizing Fixed and random model Circulating Probability Unification (FarmCPU; Liu et al., 2016) on the following traits: seed yield, seed protein concentration, seed oil concentration, plant height, 100 seed weight, days to maturity, and lodging score. Of the 200 lines included in the original panel, only 192 were included in GWAS (Supplementary Table S1). The genotypes that were excluded were Canadian cultivars and are listed in Supplementary Material. The FarmCPU model uses multiple loci linear mixed model (MLMM) and incorporates multiple markers simultaneously as covariates in a stepwise MLM to partially remove the confounding between testing markers and kinship (Liu et al., 2016). A genomic PCA matrix (P) and a genomic kinship (VanRaden) matrix (K) were used to capture the population structure and relatedness among individuals in the panel (Hyun et al., 2008; VanRaden, 2008; Li et al., 2013). Genomic Association and Prediction Integrated Tool (GAPIT; Lipka et al., 2012) was used to capture the LD decay of the SNP panel (Tang et al., 2016). An adjusted p value following methodology outlined by Brzyski et al. (2017) was used to ensure a false discovery rate (FDR) < 0.05 and to establish a significance threshold (Wang et al., 2012; Brzyski et al., 2017).

SNP markers significantly associated with a trait identified through GWAS were compared to previously reported markers and genes annotated in SoyBase Genome Browser (http://soybase.org) and NCBI RefSeq database following similar methodology to Zhang et al. (2015) to determine potential candidate genes. A length of 250 kb was added or removed from either end of the significant marker to locate potential regions for comparison based on the LD rate of the current population. In selecting candidate genes, the following criteria was used as: (i) genes of known function in soybean related to the trait under study, (ii) genes with function-known orthologs in Arabidopsis related to the trait under study, and (iii) genes pinpointed by the peak SNPs. Putative candidate genes were subsequently researched in the literature for verification.

Mean yield across environments was 2,590 ± 727.9 kg/ha, with a range of 126 kg/ha—4,805 kg/ha (Supplementary Table S3). Analysis of variance revealed that genotype and genotype-environment were the main sources of variation, with environment also showing significance, albeit of smaller magnitude (Supplementary Table S4).

The mean protein concentration observed across environments and years was 41.0% ± 1.95% (dry basis) with a range of 34–46.8%. Analysis of variance showed that genotype, environment, and genotype-by-environment effects were all significant at explaining the variation in the observed protein concentration.

The mean oil concentration across all environments was 19.8% ± 1.23% (dry basis), with a range of 14.9–23.1%. For oil concentration, genotype, environment, and genotype-by-environment effects were all significant.

Correlations were calculated to evaluate the relationships between seed yield, protein concentration, oil concentration, seed weight (g), height (cm), days to maturity, emergence score (1–10, %), and lodging score (1–5, %). Yield was found to be positively correlated with height (r = 0.47; p < 0.0001) and lodging score (r = 0.28; <0.0001), emergence sh (r = 0.56; p < 0.0001). Yield and oil concentration (r = −0.12; p < 0.0001), as well as yield and protein concentration (r = −0.08; p = 0.0035) showed significant negative correlations (Supplementary Table S5).

Protein concentration was negatively correlated with yield (r = −0.08; 0.0035), oil concentration (r = −0.40; <0.0001), height (r = −0.23; p < 0.0001), and days to maturity (r = −0.15; p < 0.0001). A significant positive relationship was observed between protein concentration and seed weight (r = 0.18; p < 0.0001). Protein concentration was not correlated with emergence nor lodging.

Oil content showed significant negative relationships with seed yield (r = −0.12; p < 0.0001), protein content (r = −0.40; p < 0.0001), height (r = −0.05; p = 0.0411), seed weight (r = −0.12; p < 0.0001), days to maturity (r = −0.26; p < 0.0001), and lodging score (r = −0.24; p < 0.0001). There was a significant positive relationship observed between oil and emergence (r = 0.07; p = 0.0088).

Correlation analysis between each location-year for seed yield, seed protein content, and seed oil contents revealed that Elora 2019 showed a significant positive relationship with Woodstock 2019 (r = 0.46; p < 0.0001); however, Elora 2019 was not correlated to either Elora 2020 or Woodstock 2020 for this trait. Yield at Elora 2020 was correlated with yield at both Woodstock 2019 (r = 0.17; p = 0.0045) and Woodstock 2020 (r = 0.16; p = 0.0017). For both protein and oil concentration, all environments were found to be correlated with each other (Supplementary Table S6).

A total of 158 million single-end reads were generated by GBS using the Ion Torrent Proton system. These reads were then mapped against the soybean reference genome (Williams82 (Gmax_275_Wm82.a2.v1); Schmutz et al., 2010) with an average success rate of 94.4%. From a total of 119,065 SNPs identified from mapping, 31, 931 SNPs remained after filtering as described in M&M. A final number of 27,911 SNP markers with minor allele frequency (MAF) > 0.05 were retained for GWAS.

Two major subpopulations, presumably corresponding to the Canadian and Chinese dichotomy, were identified (Figure 1A). The kinship matrix revealed a low level of genetic relatedness among the 200 genotypes (Figure 1B). The LD decay (r²) of the population was observed to decline to half its maximum value at 250 Kb (Figure 1C). Genomic SNP coverage for the panel of 200 soybean genotypes is depicted in Figure 1D. Evidently, the extent of LD decay varied among the different chromosome regions, resulting in uneven coverage and some regions with no SNPs identified.

GWAS was carried out using combined environment LSMEANS generated from the analysis reported in above with a total of 27,911 SNP markers used for the following traits: soybean seed yield, seed protein concentration, seed oil concentration, plant height, 100 seed weight, days to maturity, and lodging score using FarmCPU model where P + K values were used as covariates to minimize false discovery rate.

SNP markers that were significantly associated with the traits of interest are listed in Tables 1 and 2. The Manhattan plots and the corresponding Q-Q plots for these traits are depicted in Figures 2, 3. In total, 14 significant marker-trait associations were identified. Of these, only the SNPs significantly associated with soybean seed yield (three SNP), seed protein concentration (one SNP), and seed oil concentrations (one SNP) were selected for candidate gene searching. The significant SNP markers associated with the agronomic traits were excluded from candidate gene search due to resource limitations and details of those traits being out-of-scope for the current study. For the agronomic traits, a total of four SNPs were identified for 100 seed weight, two SNPs for days to maturity, two for lodging score, and one SNP for plant height.

The effect magnitudes of the minor allele on seed yield ranged from −111.2 to 162.9 (Table 1; Figure 2). One SNP was significantly associated with seed yield on Chr 5 (S05_27,809,193) while two were detected on Chr 14 (S14_5,870,227 and S14_5,884,688). However, a candidate gene was identified for only S14_5,884,688, with an effect magnitude of −108.8 (Table 1). Glyma.14 g072200 was reported to encode inositol-pentakisphosphate 2-kinase 1 (IPK1) and was identified as a potential candidate gene based on its function and proximity to S14_5,884,688. For seed protein concentration, a single SNP on Chr 5 (S05_3,040,140) was identified as significantly associated, with an effect magnitude of 0.86. Glyma.05 g03760, which encodes proprotein convertase subtilisin/kexin, was identified as a potential candidate gene for S05_3,040,140 (Table 1; Figure 2). For seed oil concentration, a single SNP on Chr 19 (S19_43,240,106; Table 1; Figure 2), with an effect magnitude of −0.40, was identified through FarmCPU. Glyma.19 g171000, which encodes zinc finger FYVE domain containing protein, was identified as the potential candidate gene for S19_43,240,106.

Allelic effect for the significant marker-trait associations for seed yield, oil, and protein concentration QTL was also measured. Based on allele frequency, phenotypic data, and effect magnitudes of the minor allele for the significant marker-trait associations identified from FarmCPU, it is likely that both the Canadian and Chinese genotypic groups may have potentially contributed the favorable allele to the seed yield QTL S14_5,870,227, S14_5,884,688, and S05_27,809,193 and the favorable allele for protein concentration QTL (S05_3,040,140) in the CD-CH group (Table 3). The Canadian group was identified as likely to have been the major contributor of the favorable allele for the seed oil QTL (S19_43,240,106) in the CD-CH group (Table 3).

For 100 seed weight, significant SNP-trait associations were detected in Chr 17 (S17_14,271,552) and Chr 18 (S18_2,625,222, S18_3,536,348, and S18_3,820,958; Table 2; Figure 3). The effect magnitudes of the minor allele ranged from −0.85 to −1.27 for this trait. Chr 18 also contained one of the two significant SNP associations for days to maturity (S18_17,449,562), with the other SNP located on Chr 15 (S15_46,719,323) for days to maturity (Table 2; Figure 3). The effect magnitudes ranged from 4.87 to 6.38 for this trait (Table 2). Only a single significant SNP was identified for plant height on Chr 5 (S05_4,738,203) with an effect of −4.19. Lastly, two SNPs were identified on Chr 9 (S09_38,461,706) and Chr 19 (S19_39,376,171) for lodging, with effect magnitudes of 0.12 (Table 2; Figure 3).