Epistatic Association Mapping for Alkaline and Salinity Tolerance Traits in the Soybean Germination Stage

Soil salinity and alkalinity are important abiotic components that frequently have critical effects on crop growth, productivity and quality. Developing soybean cultivars with high salt tolerance is recognized as an efficient way to maintain sustainable soybean production in a salt stress environment. However, the genetic mechanism of the tolerance must first be elucidated. In this study, 257 soybean cultivars with 135 SSR markers were used to perform epistatic association mapping for salt tolerance. Tolerance was evaluated by assessing the main root length (RL), the fresh and dry weights of roots (FWR and DWR), the biomass of seedlings (BS) and the length of hypocotyls (LH) of healthy seedlings after treatments with control, 100 mM NaCl or 10 mM Na2CO3 solutions for approximately one week under greenhouse conditions. A total of 83 QTL-by-environment (QE) interactions for salt tolerance index were detected: 24 for LR, 12 for FWR, 11 for DWR, 15 for LH and 21 for BS, as well as one epistatic QTL for FWR. Furthermore, 86 QE interactions for alkaline tolerance index were found: 17 for LR, 16 for FWR, 17 for DWR, 18 for LH and 18 for BS. A total of 77 QE interactions for the original trait indicator were detected: 17 for LR, 14 for FWR, 4 for DWR, 21 for LH and 21 for BS, as well as 3 epistatic QTL for BS. Small-effect QTL were frequently observed. Several soybean genes with homology to Arabidopsis thaliana and soybean salt tolerance genes were found in close proximity to the above QTL. Using the novel alleles of the QTL detected above, some elite parental combinations were designed, although these QTL need to be further confirmed. The above results provide a valuable foundation for fine mapping, cloning and molecular breeding by design for soybean alkaline and salt tolerance.


Introduction
Soil salinity and alkalinity are important abiotic stresses that adversely affect crop productivity and quality [1]. Approximately 20% of irrigated agricultural land is adversely affected by salinity and alkalinity [2], and salt-affected agricultural areas are continuously increasing. The salinity threat to agriculture exists in more than 100 countries [3]. Salinity inhibits seed germination and seedling growth; reduces nodulation; causes severe leaf chlorosis, bleaching and necrosis; and decreases total biomass and seed yield [4][5][6]. In China, there are 6.7 million ha of saline irrigated land, and 52.5-61.0% of soybean production loss is due to alkaline and salinity stresses [7]. The development of soybean cultivars with high salt tolerance is recognized as an efficient way to maintain sustainable soybean production in a salt stress environment [8][9][10]. However, the genetic architecture of the tolerance must first be elucidated. Therefore, the importance of alkaline and salt tolerance necessitates association mapping for these traits in soybean.
Although there have been some classic inheritance studies of salt tolerance in soybean [28,29], the molecular inheritance of this tolerance needs to be addressed. Lee et al. [30] identified one major QTL for salt tolerance, which was associated with markers sat_091, satt237 and satt339 on linkage group N. Based on the assumption that this QTL was identical to the Ncl locus identified by Abel [28], this locus was further confirmed in F 2 and recombinant inbred lines (RIL) [10,31]. Additional loci have been reported; e.g., Chen et al. [16] detected 11 QTL in RIL; Lee et al. [32] found a different single dominant gene in F 2 ; and Tuyen et al. [33] mapped a major QTL on linkage group D2 in F 2 and F 6 . In addition, Cho et al. [34] detected two markers, satt285 and satt079, that were significantly associated with foliar TRG accumulation, which is correlated with NaCl stress in soybean. It should be noted that all the above results in soybean were obtained from biparental populations. However, biparental population mapping strategies have several drawbacks [35]. For example, if the two parental lines do not segregate at a particular QTL, the QTL cannot be detected regardless of how many offspring are sampled in the mapping population. To overcome such shortcomings, association mapping strategies are recommended [36,37], and many such studies have been conducted [38][39][40][41][42][43][44][45][46][47]. However, very little is known about the detection of both QTL-by-environment (QE) and QTL-by-QTL (QQ) interactions.
The objective of this study was to mine novel QTL for alkaline and salt tolerances in soybean using epistatic association mapping (EAM) [36]. Elite alleles from the detected QTL were used to design parental combinations for cultivar improvement.

Evaluation of phenotypic values for STI and ATI
Alkaline and salt tolerances were measured in the LR, FWR, DWR, BS and LH of 257 soybean cultivars in 2009 and 2010, and the mean values, standard deviations, ranges, skewness and kurtosis were calculated (Tables 1 and S1). All the traits exhibited continuous distribution, and most showed a normal distribution. It should be noted that 2 germplasm accessions, Fengzitianandou and Baiqiu 1, were found to be highly resistant to salt; furthermore, 8 germplasm accessions, Linanbayuebai, Shengxiantiangengdou, Zunyizongzidou, Beichuanwuyanwo, Fengzitianandou, Guangxidalidou, Hedou 6 and Jidou 13, were highly resistant to alkaline condition.
The analysis of variance showed significant differences among all the cultivars for all the traits (Table 2), indicating that genetic variation exists among all the cultivars. Most of the other studied factors were also significant ( Table 2), suggesting that treatment, year and interactions should be considered in joint association mapping.

Mapping QTL for STI traits
A total of 83 QE interactions (24 for LR-STI, 12 for FWR-STI, 11 for DWR-STI, 15 for LH-STI and 21 for BS-STI) and one epistatic QTL of FWR-STI were detected using EAM implemented with an empirical Bayes algorithm. Among these QTL, 19 were confirmed using an enriched compressed mixed linear model (E-cMLM) method. Most of the detected QTL showed small effects on these traits, except for one LR-STI QTL associated with satt022 and one FWR-STI epistatic QTL between markers satt656 and sat_256. A summary of all the detected QTL, including the marker name, linkage group, position, variance, r 2 value and QTL type, is presented in Table 3. If two linked QTL were separated by less than 5 cM, they were considered a single QTL.
A total of 24 LR-STI QTL, with total phenotypic variance explained (PVE) values of 0.06-7.43%, were identified and mapped to linkage groups A1, A2, B1, C1, C2, D1b, E-H and J-N. Among these QTL, 6 were further identified by E-cMLM. It should be noted that the QTL associated with satt022 had a PVE greater than 5%. Furthermore, the two QTL associated with markers satt615 and sat_354 should be considered a single QTL because the genetic distance between the two markers is 2.5 cM. In addition, the QTL associated with markers satt226, satt615 and sat_354 were simultaneously identified in 2009 and 2010 using E-cMLM.
A total of 13 FWR-STI QE, with 0.17-3.10% PVE, were identified and mapped to linkage groups C2, D1b, E, G, I and M-O. One epistatic QTL with 15.80% PVE was identified between markers satt656 and sat_256, and sat_256 was also found to exhibit environmental interaction. Among these QTL, 3 were also identified using E-cMLM.
A total of 11 DWR-STI QE, with 0.13-2.05% PVE, were found and mapped to linkage groups E-G, I, K, M and N. Of these QTL, 3 were also identified using E-cMLM.
A total of 21 BS-STI QE, with 0.20-3.56% PVE, were identified and mapped to linkage groups A1-G, I, J, M and N. Of these QTL, 4 were also identified using E-cMLM. Closely linked QTL, e.g., satt452 and satt045 (1.56 cM), were considered a single QTL.

Mapping QTL for ATI traits
A total of 86 QE interactions (17 for LR-ATI, 16 for FWR-ATI, 17 for DWR-ATI, 18 for LH-ATI and 18 for BS-ATI) were identified using EAM implemented with an empirical Bayes algorithm. Among these QTL, 29 were confirmed using an E-cMLM method. Most of the detected QTL showed small effects on these traits. A summary of all the detected QTL is presented in Table 3.
A total of 17 LR-ATI QTL, with 0.41-4.20% PVE, were detected and mapped to linkage groups A2, B2, C2, D2, F and J-N. Of these QTL, 8 were also identified using E-cMLM. It should be noted that several markers were associated with both LR-ATI and LR-STI, i.e., satt615-sat_354, satt160, satt247, satt245 and satt022. Fourteen closely linked marker pairs were found to be associated with both LR-STI and LR-ATI, e.g., satt289 (LR-STI) and satt277 (LR-ATI). In addition, using E-cMLM, two markers, satt683 and sat_280, were found to be associated with LR-ATI in 2009 and 2010.
A total of 17 DWR-ATI QE, with 0.10-3.12% PVE, were identified and mapped to linkage groups A1, B1, B2, C2, D1b, D2, E, F, H, J, K, M and O. Of these QTL, 6 were also detected using E-cMLM. Several closely linked QTL were considered a single QTL, e.g., satt509 and sat_261 (0.45 cM). It should be noted that two markers were associated with both DWR-ATI and DWR-STI, i.e., sct_190 and satt220, and nine closely linked marker pairs were associated with both DWR-STI and DWR-ATI, e.g., satt441 (DWR-STI) and satt417 (DWR-ATI). In addition, using E-cMLM, marker satt132 was found to be associated with DWR-ATI in 2009 and 2010.
A total of 18 LH-ATI QE, with 0.10-3.12% PVE, were found and mapped to linkage groups A1, B2, C2, D1b, D2, E-I, K and N. Of these QTL, 5 were also identified using E-cMLM. It should be noted that three markers were associated with both LH-ATI and LH-STI, i.e., staga001, satt309 and sat_372. Three closely linked marker pairs were associated with both LH-STI and LH-ATI, e.g., satt045 (LH-STI) and satt452 (LH-ATI). In addition, using E-cMLM, marker satt452 was found to be associated with LH-ATI in 2009 and 2010.
In summary, the above mapping results produced three particularly noteworthy observations. First, almost all the detected QTL were small (1-5%). Second, almost all the detected QTL were found to exhibit environmental interactions. Finally, some markers were found to be associated with the STI and ATI indicators of multiple traits. For example, marker sat_354 was associated with LR-ATI, LR-STI, FWR-ATI, FWR-STI, BS-ATI and BS-STI, and markers satt160 and satt022 were associated with LR-STI, LR-ATI, BS-STI and BS-ATI. These data indicate that correlations exist among the above five traits ( Table 2).

Mapping QTL from original traits
Original trait observations and marker information were used to perform EAM. A total of 78 QE (17 for LR, 14 for FWR, 4 for DWR, 21 for LH and 21 for BS) and 3 epistatic QTL for BS were detected and are listed in Table 4. Among these QTL, 74 were also identified using E-cMLM. Most of the detected QTL showed small effects on these traits, except for two epistatic QTL: sat_390 6 sat_254 and satt244 6 sat_160.
A total of 17 QE interactions for LR, with 0.07-2.38% PVE, were detected and mapped to linkage groups C2, D1a-F, H, J-L, N and O. Of these QTL, 16 were also identified using E-cMLM. Note that seven markers, satt289, sat_252, satt672, satt263, satt149, sat_224 and satt244, were found to be associated with LR in more than two environments. When all the LR datasets were jointly analyzed by E-cMLM, 15 QTL-by-year interactions were detected, which was consistent with the identification of 15 QE interactions using EAM.
A total of 14 FWR QE interactions, with 0.06-1.06% PVE, were identified and mapped to linkage groups A1-B1, C2, D2-F and J-N. Of these QTL, 13 were also identified using E-cMLM. Note that nine QTL were found to be associated with LR in more than two environments. When all the FWR datasets were jointly analyzed using E-cMLM, 13 main-effect QTL, one QTL-bytreatment interaction and QTL-by-year interaction were detected, which was consistent with the identification of 13 QE interactions using EAM.
Four DWR QE interactions, with 0.08-4.38% PVE, were identified and mapped to linkage groups B1, D2, H and K. Of these QTL, one was also identified using E-cMLM.
A total of 21 LH QE interactions, with 0.07-2.00% PVE, were identified and mapped to linkage groups B1, B2, C2-L and N. All of these QTL were also identified using E-cMLM. Note that nine QTL were found to be associated with LH in more than two environments. When all the LH datasets were jointly analyzed using E-cMLM, 20 main-effect QTL and five QTL-by-year interactions were detected. Among these QTL, twenty were found to be consistent with the QE interactions for LH index identified using EAM.
A total of 21 QE and 3 QQ interactions for BS, with 0.08-11.41% PVE, were detected and mapped to linkage groups A2, B1, C2-F, H and J-O. Of these QTL, 23 were also identified using E-cMLM. Several closely linked QTL were considered a single QTL, e.g., satt413 and satt672 (1.36 cM), and thirteen QTL were found to be associated with BS in more than two environments. When all the BS datasets were jointly analyzed using E-cMLM, 26 main-effect QTL and two QTL-by-year interactions were detected. Among these QTL, 23-were found to be consistent with the QE and QQ interactions for BS identified using EAM.
A comparison of the QTL for ATI and STI indicators revealed 56 common QTL: 14 for LR, 11 for FWR, 3 for DWR, 15 for LH and 13 for BS. For example, markers satt289, satt222 and satt022    were associated with both LR and LR-STI, and markers satt683 and satt022 were associated with both LR and LR-ATI. Among these markers, marker satt022 is common; in other words, marker satt022 was associated with LR, LR-STI and LR-ATI.

Predictions for novel parental combinations
The best way to improve a trait is to pyramid all the desirable elite alleles into one cultivar, if possible. By maximizing the number of elite alleles using a Monte Carlo simulation experiment, the ideal novel cultivar combination can be designed. Using this method, 21 elite alleles of 27 QTL for LR-STI were pyramided by combining cultivars Guangxibayuehuang, 0804, Shangqiu 832012, Qingyuanxiaoqingdou, Daheidou, Ganjiangnan, Shuangliuliuyuehuang, Baiqiu 1, Xu 0701 and Wenfeng 5.

Discussion
The mapping results from this study are reliable for three reasons. First, 19 QTL for STI, 29 QTL for ATI, and 74 QTL for original traits, which were detected using EAM, were confirmed by E-cMLM (Tables 3 and 4). Furthermore, 36 QTL for STI (11 for LR-STI, 7 for FWR-STI, 5 for DWR-STI, 4 for LH-STI and 9 for BS-STI), 6 QTL for ATI (2 for LR-ATI, 1 for FWR-ATI, 1 for DWR-ATI and 2 for BS-ATI) and 29 QTL for original traits (6 for LR, 6 for FWR, 10 for LH and 7 for BS) found in this study have also been identified by other researchers (Tables 3 and 4). For example, one major QTL on linkage group N plays an essential role in enhancing soybean salt tolerance in different genetic backgrounds [10,28,30], and the present study also confirmed several associations with salt tolerance, e.g., satt237 for BS-STI and satt234 for LR, LH, FWR and BS. Chen et al. [16] identified 11 QTL significantly associated with salt tolerance, and Cho et al. [34] detected 15 QTL on linkage groups B2, C2, D2, G, J and K for foliar TRG accumulation, which is postulated to function as a compatible solute and/or osmoprotectant under adverse NaCl stress conditions. These 26 QTL were all confirmed in this study. Finally, several salt tolerance genes from Arabidopsis thaliana and soybean (http://www.ncbi.nlm.nih.gov/ gene?term = salt%20tolerance) were found to be located in close proximity to tolerance-associated markers in this study.
The potential candidate genes for salt tolerance are summarized in Tables S2 and S3. Among these candidate genes, four soybean genes, Glyma13g41980, Glyma15g03400, Glyma17g37430 and Glyma17g35340, which have homologs in Arabidopsis, were associated with five markers examined in this study, with physical distances of 135.76-820.04 kb. Six soybean salt tolerance genes, LOC100808889, LOC100807827, LOC100800981, LOC100795117, LOC100814727 and LOC100797515, were closely linked to seven markers analyzed in this study, with physical distances of 135.41-619.74 kb. Although these consistent results were observed, one might raise doubts about our conclusions because we did not consider cultivar relatedness. In our opinion, the inclusion of all the main and epistatic QTL in the genetic model reduced the importance of controlling for genetic background.
Tuyen et al. [33] reported that tolerance to NaCl may not always accompany tolerance to alkaline stress. Their evidence was based on a discrepancy between the QTL for alkaline tolerance, located between markers satt669 and sat_300 on linkage group D2, and the QTL for saline tolerance, located on linkage group N, in wild soybean JWS156-1. However, the opposite phenomenon was observed in the present study. For example, satt237 on linkage group N was associated with both LR-ATI and BS-STI; satt615 on linkage group D2 was associated with both LR-STI and LR-ATI; and sat_354 on linkage group D2 was associated with LR-STI, LR-ATI, FWR-STI, FWR-ATI, BS-STI and BS-ATI.
Similar results for other linkage groups were also found in this study. More importantly, the cultivar Fengzitianandou was found to be highly resistant to NaCl and alkaline stresses. In addition, tolerance to salt and alkaline stresses is related to tolerance to Al in soybean.
Sharma et al. [48] identified two major QTL for tolerance to Al. These two QTL were also detected in the present study and were located in the marker intervals satt160-satt252 (F) and satt202-satt371 (C2), respectively. Note that satt160 was simultaneously associated with LR-STI, LR-ATI, FWR-ATI, DWR-STI, LH-STI, BS-STI, BS-ATI and LH. Bianchi-Hall et al. [49] identified six QTL associated with tolerance to Al stress. Among these QTL, Al tol1-1, Al tol1-3, Al tol1-5 and Al tol1-6 were found to be located in close proximity to the associated markers in this study. For example, satt329 (associated with LR-ATI) and satt233 (associated with BS) on linkage group A2 were close to Al tol1-1; satt509 (associated with LR-STI, DWR-ATI and FWR) on linkage group B1 was close to Al tol1-3; satt215 (associated with LR-STI, FWR-ATI and FWR), satt431 (associated with LR-ATI), satt132 (associated with DWR-ATI, LH-STI and BS), satt244 (associated with LR and BS) and sat_224 (associated with LR) on linkage group J were close to Al tol1-5; and satt238 (associated with LR-ATI, FWR-ATI, FWR and BS) on linkage group L was linked to Al tol1-6. These results demonstrate that the mechanisms of tolerance to alkaline/salt and Al conditions are similar.
Most of the detected QTL showed small effects, except for four QTL: one main-effect QTL associated with satt022 and three epistatic QTL: satt656 6 sat_256 for LR-STI, sat_390 6 sat_254 for BS and satt244 6 sat_160 for BS. The results of this study demonstrate that large differences in tolerance are caused not only by a few large-effect QTL but also by the cumulative effect of numerous small-effect QTL. This small-effect QTL phenomenon has been observed in Arabidopsis [50], rice [51] and maize [52], although the high PVE of the epistatic QTL is surprising because alkaline and salt tolerances in plants are derived from interactive molecular pathways [53]. In addition, almost all QTL exhibit significant environmental interactions, indicating that the genetic architecture across various salt and alkaline stresses is sensitive or variable.
Research into tolerance to NaCl and alkaline stresses usually employs two types of indicators: a salt and alkaline tolerance index and an original trait indicator. These two types of indicators were also analyzed in this study, producing complementary results: tolerance genes were found proximal to both types of associated markers.
The utilization of varieties with multiple resistance genes is an effective way to reduce the effects of adversity. As the numbers of both cloned genes and detected QTL increase, a Monte Carlo simulation experiment becomes a valuable means of parental combination prediction and selection strategy design [36,54]. In this study, several parental combinations were predicted. Certain cultivars are encountered repeatedly in these combinations and may be used to improve multiple traits. For an STI indicator, e.g., Sudou 1 may be chosen for improving FWR and BS; Daheidou for LR and LH; Qingyuanxiaoheidou for LR, LH and BS; Qinyan 1 for DWR and BS; and Wuhuabayuehuang for DWR and BS. These predictions may be valuable for several reasons. First, two of the selected parents, Zhejiangsiyuebai (DWR-STI) and Caoqing (BS-STI), are on the list of 348 ancestors of 651 soybean cultivars released during 1923-1995 in China [55]. Youbian 31, Sudou 1, Fengjiao 66-12, Ludou 1, Wenfeng 5, and Jindou 2 are on the list of 171 varieties that were bred earlier and were derived from 519 soybean cultivars [55]. Second, some accessions exhibit beneficial agronomic traits and strong resistance to adversity, i.e., Sudou 1, Ludou 1, Kexi 8, Daheidou, Edou 2, Niumaohuang, Wenfeng 5, Dianbaiheidou and Baiqiu 1. Finally, some selected cultivars have been widely grown in certain areas, i.e., Sudou 1, Ludou 1 and Wenfeng 5. In addition, we should consider the sizes of novel alleles when predicting parental combinations. Thus, the tolerant variety Fengzitianandou may be included in a parental combination because of its large effect. Of course, a prerequisite for the above prediction is that all the above QTL will be further confirmed.
The phenomenon of QTL clusters has been reported in soybean [56] and cotton [57]. Previous work has indicated that numerous disease resistance loci are clustered in various regions of the soybean genome, e.g., on chromosomes D1b and F. This phenomenon was also evident in the current results. For example, QTL linked to markers satt226, satt615, sat_354 and satt514 on chromosome D2 were clustered in an 8.55 cM interval region. Three salt tolerance genes, Glyma17g37430, Glyma17g35340 and LOC100807827, were found to be located near the toleranceassociated QTL clusters of satt256, satt413 and satt672 (Tables S2  and S3). The phenomenon of QTL clusters may be used to explain trait correlation [57].

Plant materials and DNA marker analysis
All the 257 soybean cultivars were obtained, by stratified random sampling, from 6 geographic ecotypes in China (Fig. 1), were planted in three-row plots in a completely randomized design, and were evaluated at the Jiangpu experimental station at Nanjing Agricultural University in 2009 and 2010. The plots were 1.5 m wide and 2 m long. These seeds were used to conduct plastic container experiments.
Approximately 0.3 g of fresh leaves obtained from each cultivar in 2009 was used to extract genomic DNA using the cetyltrimethylammonium bromide method, as described by Lipp et al. [58]. To screen for polymorphisms among all the cultivars, PCR was performed with 135 simple sequence repeat (SSR) primer pairs. The primer sequences were obtained from the soybean database Soybase (http://www.ncbi.nlm.nih.gov). PCR was performed as described by Xu et al. [56].

Evaluation of alkaline and salt tolerances
A salt-water flooding method [59] was used to evaluate the alkaline and salt tolerances of all the soybean cultivars. In brief, twelve soybean seeds for each cultivar were sown in a 30620615 cm plastic container with sand added to a height of 3.56cm. The seeds were then treated with control (CK, pH: 7.0), 1006mM NaCl (pH: 7.0) and 106mM Na 2 CO 3 (pH: 11.1) solutions, with two replications each. A 350-ml aliquot of the appropriate solution for each treatment was applied to each plastic container filled with sand. Twelve soybean seeds for each treatment were grown in a growth chamber under white fluorescent light (600 mmol m 22 s 21 ; 14 h light/10 h dark) at 2561uC. The length of the main root (LR), fresh and dry weights of roots (FWR and DWR), biomass of seedlings (BS) and length of hypocotyls (LH) of healthy seedlings from 5 plants in plastic containers under simulated alkaline and salt conditions were measured 7 days after sowing. The units used were centimeters for length and grams for weight. To measure the degree of salt and alkaline tolerances, the original trait observations may be transformed into salt and alkaline tolerance indices for each trait using the below equations: alkaline tolerance index ATI ð Þ~x CK {x Na 2 CO 3 .
x CK |100% where x CK , x NaCl and x Na 2 CO 3 stand for the phenotypic values exhibited following control, saline and alkaline treatments, respectively.

Population structure
Population structure plays an important role in association mapping. To investigate the population structures of all the selected cultivars, the STRUCTURE program [60] and the approach of Evanno et al. [61] were employed. The number of subpopulations (K) was set from 2 to 10. The number of replicates for each K was 20, and the total average of the mean log-likelihood at a fixed value of K was used. Using the ad hoc statistic DK, which is based on the rate of change in the log-  probability of data between successive K values, the DK value was much higher for the model parameter K = 4 than for other values of K (Fig. 2). By combining this high DK value with knowledge of the breeding history of these cultivars, a K value of 4 was chosen. The Q matrix was calculated based on information from 135 SSR markers and was incorporated into the association mapping.

Epistatic association mapping
The EAM approach suggested by Lü et al. [36] was also used to analyze all the five-trait datasets in this study; for technical details, the reader is referred to the original study by Lü et al. [36]. All the above analyses were performed using the SAS program.
The enriched compression mixed linear model (E-cMLM) approach, suggested by Li [37] and expanded by Zhang et al. [62], was used to confirm the results of the epistatic association mapping. In this analysis, the dataset for each year and each treatment was analyzed, and all the datasets for all the twoyear and three-treatment datasets were jointly analyzed.

Supporting Information
Table S1 Phenotypic variation in the length of the main root (LR), fresh and dry weights of roots (FWR and DWR), biomass of seedlings (BS) and length of hypocotyls (LH) among healthy seedlings measured in 257 soybean cultivars in 2009 and 2010. (DOC)