Soybean is one of the most important oilseeds as well as protein source crop worldwide. In Kazakhstan, soybean planting area has been increased from 45 thousands hectares in 2006 to 120 thousands in 2016 [1]. This massive enlargement is under crop diversification policy adopted by KZ government with the final goal to reach up to 400 thousand hectares by 2020 occupied by commercial soybean annually. Designated areas of the expansion is land devoted to agriculture in south-eastern (SEK), eastern (EK) and northern (NK) regions of KZ [2, 3].

Soybean is a new crop for Kazakhstan. It dictates the necessity of preliminary studies such as evaluation of a number of diverse soybean varieties grown in the new environment and their potential use as an acceptor germplasm for the specific conventional breeding purposes. Recent findings regarding non-random associations between certain alleles of flowering genes and yield in soybean lines grown in these targeted regions [1] were a promising start in this direction. In this study, the time flowering span was assessed in all three regions in association with variation in E genes and yield performance. Specific allele combinations of the four E genes and respective optimal ranges of flowering and maturity time were identified for each experimental site [1]. However, more comprehensive research-based support is needed for successful development and implementation of specific conventional soybean breeding programs in KZ.

Several genomic oriented tools were generated in soybean community in recent years to facilitate breeding programs worldwide (soybase.org). The list of these genomic tools includes assembly of Williams 82 genomic sequence (http://www.soybase.org /SequenceIntro.php), Affymetrix SoyChip annotation [4], searching engines for Glycine max and Glycine soja sequences (http://plants.ensembl.org/Glycine_max/Info/Index), SoySNP50K iSelect BeanChip from Illumina [5], and etc. The 50,000 (50 K) SNP iSelect array was found to be particularly instrumental in genetic mapping of QTL (quantitative trait loci) for complex agronomic traits, including abiotic [6, 7] and biotic stress resistances [8] and seed quality [9, 10].

Genome-wide association studies (GWAS) is considered to be one of the most promising approaches in the identification of QTL of agronomic traits [11, 12]. In soybean GWAS experiments indicated a strong bias towards environmental factors in MTA discoveries [13–15]. For instance, obtained results from three different GWAS studies related to the identification of QTLs for yield performance in Canada [16], USA [13, 15], and Brazil [11] showed different responses and QTLs for yield components were identified in different parts of the genome. Presumably, in order to generate a reliable data for regionally running breeding projects, separate MTA experiments are required. The main objective of this study was to identify non-random MTAs in soybean field trials in three different environments in KZ. This is the first attempt based on association mapping approach for identification of soybean productivity with related QTLs in Kazakhstan.

113 soybean accessions from 5 different breeding origins were compared for yield performance in three soybean growing regions of Kazakhstan. North American accessions demonstrated the highest yield in the majority of trials (Fig. 1), supporting superior adaptability of this genetic material to the KZ environment. It is interesting to note that the pedigree of a few local high-performance cultivars also included in this study are heavily based on the US bred soybean background. On the other hand, principal coordinate analysis (PCA) showed that the most of local varieties genetically closer to East European germplasm, although North American and East European accessions were not far apart from each other (Fig. 3).

The size and level of genetic variation in studied genetic panels appear to be critical for the positive outcome of GWAS based projects. It was demonstrated that experiments with sample size less than 384 accessions [17] and large LD blocks [18, 19] may lead to the identification of false positive associations. On the other hand, in the study by Turner et al. 2016, there was shown that smaller panels may allow detection of false negative associations that would not have been detected in the larger panels [20]. These study results on relatively small soybean sample size (n = 113) in fact confirm the above mentioned findings. Initially, 64 MTAs were detected using TASSEL MLM model with the application of K + Q matrices. However, further validation of the revealed associations using t-test suggested that only 42 of them are presumably true associations. GWAS led to the identification of 30 SNP markers involved in 42 confirmed MTAs. Comparison of physical positions of identified SNPs with positions of previously mapped QTLs (https://soybase.org/search/qtllist_by_symbol.php) showed that 32 out of 42 identified MTAs have been found to be at the same locations as known QTLs. Remaining 10 MTAs identified in this study might be denoted as presumably novel QTLs. Also, a high percentage of positive matches regarding physical locations of identified SNPs and known QTLs support GWAS as a useful research tool for the search of non-random MTAs in soybean in Kazakhstan.

The other causative factor perhaps relevant to the identification of novel QTLs in this study is the impact of specific environmental conditions listed in Additional file 5. Here significant environmental influence (81.9%) on yield performance was detected by ANOVA in the genotype x environment interaction (GEI) assessment. Four MTAs were significant in two or more field trials, and this can be apparently explained by large environmental effect causing a relatively low stability of the identified loci.

In this study, 24 MTAs related to plant adaptation traits including plant height were identified. For the length of flowering time, two MTAs were revealed on chromosomes Gm06 and Gm19 (Table 3). The Gm06_20,370,075 was located in the close proximity to the GmPhyA1 gene, which was annotated as a phytochrome receptor corresponding to the flowering E1 locus [21]. Gm19_49,964,637 has located approximately 2.5 Mb from the GmPhyA3 gene, another phytochrome receptor corresponding to the maturity E3 locus [22]. Four identified SNPs, including those on chromosomes Gm10 and Gm19, were associated with maturity time, and their map locations coincided with the genetic positions of E2 and E3 maturity loci [21, 22]. Gm10_48,586,134 in chromosome Gm10 was approximately 3.9 Mb apart from GmGIa (Glyma10g36600), a gene that involved regulation of soybean maturity and flowering time [23]. The Gm19_48,168,077 in chromosome Gm19, as in case of flowering MTA, was located approximately 1.5 Mb from the GmPhyA3 gene (Glyma19g41210). It is interesting to note, that location of another SNPs Gm07_10009107 was in proximity to the Phytochrome B (Glyma07g11790) gene, which commonly encodes inhibitory effect on plant germination [24]. SNP location of the Gm09_42,241,644, which is related to the MTA for PH, was coincided with the physical position of previously mapped QTL for internode length (soybase.org) (Additional file 3).

The second group of identified MTAs belongs to yield component traits, which includes NFN, NSP, TSW, and YPP. In total, GWAS allowed the identification of 20 MTAs for yield components spread on ten different chromosomes (Table 3).

Overall thirty SNP markers associated with time to flowering and maturation, plant height, number of fertile nodes, seeds per plant and yield were identified. Physical locations of 32 identified out of 42 total MTAs coincide well with positions of known analogous QTLs (www.soybase.org). This result indicates importance of revealed MTAs for soybean growing regions in Kazakhstan. The other 10 MTAs might lay claim to some novelty until additional proof is obtained. Obtained results would serve as required prerequisite for forming and realization of specific breeding programs towards effective adaptation and increased productivity of soybean in three different regions of Kazakhstan.