Unrevealing the Genetic Basis of Alkaline Tolerance in Soybean Using Genome-Wide Association Study

1. Introduction

Soybean [Glycine max (L.) Merr.] is cultivated and consumed globally because its seeds possess abundant quantity of edible oil and protein [1]. Moreover, the importance of soybean in the improvement of soil fertility, as well as its significance in human health has been well recognized [2]. Globally, the soybean consumption is maximum in China; however, the country is heavily reliant on soybean imports [3,4]. However, there are many environmental challenges including both biotic and abiotic stresses which limit the soybean production in China. Among, these challenges, alkaline stress has been considered as a major detrimental factor for soybean productivity [5,6]. It has been documented that China possesses around 100 million ha of salt-alkaline soil, and there is a need to bring this land under crop cultivation [7]. Hence, production of the alkaline-tolerant cultivars of soybean is an important breeding objective. Studies have documented the significant correlation of the alkaline stress with the shoot and root biomass related traits in soybean [8]. However, the production of alkaline-tolerant cultivars of soybean needs the detailed information of genetic basis that regulates the alkaline tolerance in soybean. Comprehensive understanding of the genetic elements and genes controlling the alkaline tolerance in soybean will allow their use in the production of alkaline-tolerant cultivars of soybean through molecular breeding [7,9].

Previously, the efforts put forward on unraveling the genetic elements and genes governing the alkaline tolerance in soybean was negligible, and only four QTLs associated with alkaline stress tolerance have been documented in SoyBase (https://www.soybase.org/). However, the traditional QTL mapping involving low-throughput markers were utilized to identify these QTLs [5]. The traditional genetic mapping possessing the demerit of low-resolution has restricted the use of these genetic elements in the breeding of alkaline tolerance in soybean [10]. In this context, the association mapping approach possessing high resolution and precision has emerged the practical approach in marker-assisted breeding applications [11,12]. Recent breakthrough seen in the genome sequencing and high-throughput genotyping has allowed the routine application of association mapping in crop plants [13,14]. Multiple studies have documented the use of GWAS for different abiotic stresses in multiple plant species, such as drought tolerance in rice [15], salinity tolerance in wheat [16], alkaline tolerance in maize [17], drought tolerance in chickpea [18], salinity tolerance in mungbean [19], heat stress tolerance in Faba bean [20] and pea [21].

In the crop germplasm, every gene has multiple alleles and these all alleles cannot be identified by using the biallelic SNP markers. Some of these alleles can be rare and superior resulting their loss in the GWAS identification process [12]. In this context, the emergence of haplotype markers that are multi-allelic can capture the rare and superior alleles in the crop germplasm. The performance of haplotype markers has been observed better than the SNP markers by both empirical and simulation methods [22,23,24]. The better efficiency of haplotype markers compared to SNPs in association mapping and genomic selection (GS) analyses can be attributed to the multi-allelic nature of the haplotype markers [12]. The desirable haplotypes for abiotic stress tolerance have been identified in different crop species salinity tolerance in soybean [25], salt tolerance in wheat [26], drought tolerance in wheat [27].

In soybean, GWAS has been performed for different traits, including abiotic stress related traits. For example, GWASs have been performed for salt tolerance using the 305 diverse soybean accessions that were genotyped with the Illumina Infinium SoySNP50K BeadChip [28]. Similarly, GWAS analysis of the salt tolerance trait was carried out in a global soybean panel using the SoySNP50K BeadChip for genotyping [29]. GWAS mapping has been carried out for drought tolerance using diverse of panel of 259 released Chinese soybean cultivars genotyping with the Illumina SoySNP6k iSelect BeadChip [30]. The GWAS analysis has been performed for seed-flooding tolerance in soybean by using the Yangtze-Huai soybean breeding germplasm (YHSBG) population consisting of 347 diverse soybean accessions genotyping with the RAD-seq method [31]. GWAS analysis was also carried out for cold tolerance in the soybean at the germination stage using 260 soybean accessions genotyped with the 30,799 SNPs [32]. However, very limited studies of GWAS have been reported for the alkaline stress tolerance in soybean [10,33]. Hence, concerted efforts are needed to identify highly consistent and stable genetic elements as well as true functional genes governing alkaline tolerance in soybean.

To discover the alkaline tolerance at the genetic level, the association mapping panel of 326 diverse soybean accessions were used to perform genetic mapping, haplotype, and candidate gene analyses. The outcome of the present study involves the major QTLs, haplotypes, and candidate genes linked to soybean alkaline tolerance. These genetic elements can be incorporated in the breeding of alkaline tolerance in soybean, and the candidate genes can be validated for their further deployment in soybean breeding.

2. Results

2.1. Determination of Optimum Alkaline Treatment

To determine the optimum level of alkaline treatment that can be used for the evaluation of alkaline stress tolerance in the GWAS panel of 326 soybean accessions, we tested the performance of the two soybean cultivars viz., DS118 and JY86 under the different alkaline treatments viz., 50 nmol/L, 100 nmol/L, 150 nmol/L and 200 nmol/L. The DS118 is a well-known alkaline-tolerant cultivar, whereas JY 86 is a local model soybean cultivar in the Jilin province, China. Our results revealed that the above-ground shoot of DS118 cultivar grows well under the alkaline treatments of 50 mmol/L, 100 mmol/L and 150 mmol/L (Figure 1). However, the plant height of DS118 cultivar gradually decreased with the increase in the alkaline treatment concentration from 50 mmol/L to 200 mmol/L (Figure 1); besides at the alkaline concentration of 200 mmol/L, the leaves of the plant turned yellow as well as are mildly wilted. The underground root growth was also inhibited to different degrees with the increase in the alkaline treatment concentration, and the inhibition was more obvious when the concentration of alkaline treatment was increased beyond 150 mmol/L (Figure 1). However, in the case of JY86 cultivar the phenotype of the shoot growth seems good at the alkaline treatment of 50 mmol/L, but the plant leaves began to turn yellow as the concentration of alkaline treatment reached 100 mmol/L (Figure 1). Moreover, at the alkaline treatment of 150 mmol/L the leaves were wrinkled and curled as well as leaf number was reduced (Figure 1). The underground root of JY86 cultivars was obviously inhibited under 150 mmol/L alkaline treatment, and the roots are dying at the alkaline treatment of 200 mmol/L. Based on the phenotypic characterization, it become evident that local soybean cultivar viz., JY86 are more sensitive to alkaline stress relative to the DS118 cultivar. By considering our results, the 150 mmol/L was considered as an optimum concentration for the evaluation of GWAS panel of soybean accessions for alkaline stress tolerance.

2.2. Phenotypic Analysis of Alkaline Tolerance-Related Traits

The alkaline tolerance of 326 soybean accessions of GWAS panel was evaluated using eight traits of alkaline tolerance-related traits viz., seedling fresh weight (SFW), seedling dry weight (SDW), root dry weight (RDW), root fresh weight (RFW), relative chlorophyll content (CC), root number (RN), root tips number (RTN) and root length (RL). The phenotypic appearance of the four representative cultivars under control (CK) and alkaline treatment (AT) are present in the Figure S1. The values of descriptive statistics paramaters for all the eight studied traits are presented in Table 1. The ANOVA (F-value) analysis showed that genotype (G), environment (E) and G × E interaction were highly significant (P<0.0001) for all the eight studied traits, expect that G and G × E varaition were non-significant for the RFW trait (Table S1).

2.3. Population Structure and LD Analysis

By applying the quality control measures viz., MAF<0.05 and a missing genotype at 0.01, we obtained a total of 3,311,166 high-quality SNPs (Table S2). These SNPs were distributed across all the 20 chromosomes in soybean; the lowest number of SNPs were located on the Chr.11 (62,554), whereas the highest number of SNPs were present on the Chr.15 (254,495) (Figure 2A; Table S2). Markers density across the different chromosomes varies considerably, and the highest marker density were observed on the Chr.15 (4916.83 SNPs/Mb), and lowest marker density were observed on the Chr.11 (1799.08 SNPs/Mb) (Figure 2B; Table S2).

The analysis of population structure also revealed no distinct structure with a continuous distribution (Figure 3A; Table S3). In addition, the kinship matrix represented by the dendrogram and heatmap does not showed the clear clustering of the soybean accessions (Figure 3B). The linkage disequilibrium (LD) characteristic was also studied for the association mapping panel of 326 accessions, and the results are shown in the graph (Figure 3C). The mean value of r² in the soybean genome was 0.68, and at the 0.68 the decay of LD started and at the 0.34 the half decay is reached. At the 71.6 kb the LD decay curve and half decay intersected, and this distance determines the linkage across the whole soybean genome. The genomic interval of ±71.6 kb in the upstream and downstream of the stable significant SNPs are considered as QTL.

2.4. Association Mapping Analysis of Alkaline Tolerance

In the current study, we identified a total of 28 SNPs significantly associated with the six alkaline tolerance-related traits viz., RDW, RFW, SFW, CC, SDW and RTN across the two environments (Table 2; Figure 4 & Figure 5). These significant SNPs include one, one, four, seventeen, four and one for the traits viz., RDW, RFW, RTN, SDW, SFW and CC, respectively. Besides, twenty and eight significant SNPs were identified in the AT and ACR environments, respectively. These significant SNPs identified in the AT and ACR environments for different alkaline tolerance-related traits are associated with alkaline tolerance in soybean. Among these 28 significant SNPs, the Chr04_51929177 was identified by all the seven GWAS models viz., BLINK, FarmCPU, MLMM, super, GLM, MLM and CMLM in the ACR environment. The significant SNPs viz., Chr20_25660093 was identified by the five GWAS models viz., MLMM, super, GLM,

MLM and CMLM as well as in both environments viz., AT and ACR. Moreover, the SNPs viz., Chr18_27591088, Chr04_52367033, Chr04_51934424 and Chr04_52131534 were identified by four out of total seven GWAS models in the ACR environment. The Chr01_38897254 was identified by three GWAS models viz., MLMM, super and CMLM in the ACR environment. The remaining 21 significant SNPs identified in the AT or ACR environments were detected by either two or one GWAS models.

Based on the above results, the significant SNPs viz., Chr01_38897254, Chr18_27591088, and Chr20_25660093 are identified through multiple GWAS models and environments (AT & ACR), and thus considered as stable marker-trait associations (MTAs) regulating the alkaline tolerance in soybean. By considering the upstream and downstream distances within the linkage disequilibrium (LD) decay (±71.6 kb) around these three stable significant SNPs, we delineated these three genomics regions as stable QTLs viz., qAT1, qAT18 and qAT20, respectively. Besides, 15 significant SNPs falling within the LD decay distance of ±71.6 kb are identified on the Chr.14 to be associated with SDW in the AT environment across two GWAS models, this region can be also considered as the stable genomic loci/QTL associated with alkaline tolerance in soybean, and are named as qAT14. Moreover, the four significant SNPs viz., Chr04_51929177, Chr04_52367033, Chr04_51934424 and Chr04_52131534 that are associated with the SFW in the ACR environment are identified within the LD decay distance of ±71.6 kb on Chr.04 across the all seven GWAS models, this genomic region can be also delineated as QTL named as qAT4. Overall, we identified a total of five stable QTLs viz., qAT1, qAT4, qAT14, qAT18 and qAT20 associated with the alkaline tolerance in soybean.

2.5. Candidate Genes Identification

A total of 39 model genes were identified within the physical genomic interval of five major QTLs viz., qAT1, qAT4, qAT14, qAT18 and qAT20, which consisted of four, 14, 11, one and nine genes, respectively (Table 3). Furthermore, based on the in-silico analysis, gene annotations and literature survey, we defined a total of 17 genes underlying the genomic interval of five QTLs, that include one, eight, three, one and four genes within the physical interval of qAT1, qAT4, qAT14, qAT18 and qAT20. The selection of these candidate genes was based on gene function annotation such as genes governing salt tolerance, drought tolerance, protein phosphorylation, heat shock proteins, basic helix-loop-helix (bHLH) transcription factors, heat shock transcription factor, S-ribonuclease binding protein 1, GRAS family transcription factor, pentatricopeptide repeat (PPR-like) superfamily protein, C2H2-like zinc finger protein, abiotic stress tolerance, auxin signaling, kinase activity, ethylene signaling and salicylic acid signaling. These 14 genes possess one of the above function, thus were considered as putative candidate genes regulating alkaline tolerance in soybean.

2.6. Haplotypes Identification for Alkaline Tolerance

Our results showed that the four significant SNPs identified on Chr.04 associated with the SFW in the ACR environment were present within the LD decay range of ±71.6 kb, and these SNPs thus form haplotype blocks viz., Hap4 on Chr.04 (Figure 6A). Besides, the 15 significant SNP identified on the Chr.14 are associated with SDW in the AT environment fall within the the LD decay range of ±71.6 kb, thus form another haplotype block viz., Hap14 (Figure 6B). The haplotype block Hap4 possesses three haplotype alleles viz., Hap4_1, Hap4_2 and Hap4_3 (Figure 6A). The Hap4_2 regulates the higher value of SFW under the ACR environment, thus regulates the higher alkaline tolerance; whereas the Hap4_1 regulates the lowest value of SFW under the ACR environment and regulates the lower alkaline toerlance, besides, the Hap4_3 regulates the intermediate value of SFW under the ACR environment, thus regulates the moderate alkaline tolerance in soybean. The haplotype block Hap14 possesses two haplotype alleles viz., Hap14_1 and Hap14_2 (Figure 6B). Haplotype allele viz., Hap14_1 regulates the lowest value of SDW under the AT environment, thus are regulating the lower alkaline tolerance in soybean; whereas the Hap14_2 regulates highest value of the SDW under the AT environment, thus are regulating the higher alkaline tolerance in the soybean.

Moreover, the haplotype analysis of the 17 candidate genes revealed that haplotype alleles of the four candidate genes viz., Glyma.04G252300, Glyma.04G253100, Glyma.14G083700 and Glyma.20G072500 showed signigicant difference in the regualation of alkaline tolerance-related traits viz., RDW, RFW, SDW, SFW, RTN & CC under the treatments of either AT or ACR or both (Figure 7). The haplotype alleles underlying these four candidate genes across the GWAS panel of accessions varies from two to four. These haplotype alleles governs the alkaline tolerance from the lowest to highest through moderate levels.

The immobile nature of plants has exposed them to environmental constraints; however, the alkaline stress has emerged one of important constraint that negatively affects the plant growth and development as well as productivity and overall production. Therefore, production of alkaline-tolerant cultivars of soybean has received increased attention among the soybean researchers. In the present study, eight traits related to alkaline stress tolerance viz., seedling fresh weight (SFW), seedling dry weight (SDW), root dry weight (RDW), root fresh weight (RFW), relative chlorophyll content (CC), root number (RN), root tips number (RTN) and root length (RL) were evaluated in the GWAS panel of 326 accessions of soybean under alkaline treatment. The highly significant correlation was observed among these eight alkaline tolerance-related traits in the GWAS panel under all three environments viz., CK, AT and ACR. Overall, the comprehensive analysis of RDW, RFW, SFW, CC, SDW, RN, RTN and RL is an effective way for the evaluation of alkaline tolerance at the seedling stage.

Almost negligible efforts have been put forward to elucidate the genetic basis of alkaline tolerance in soybean. Till data only one study has attempted to identify the QTLs for alkaline tolerance in soybean [5], and this study has identified only four QTLs viz., pH sensitivity 1-1, pH sensitivity 1-2, pH sensitivity 1-3 and pH sensitivity 1-4. These QTLs have been detected through traditional mapping approach involving low-throughput markers, which has resulted their low resolution as well as prevent their use in the soybean breeding [9]. Recent breakthroughs in the genome sequencing, genotyping platforms and GWAS approach has provided the avenues for the high-resolution mapping of genetic loci and underlying genes for the particular trait of interest [34]. In soybean GWAS approach has been used to explore the genetic basis of multiple traits [1,35]; however, negligible efforts of GWAS study have been attempted to unravel the genetic architecture of alkaline tolerance in soybean. By keeping the above into consideration, the GWAS approach involving the high-throughput marker system were used to determine the genetic architecture of alkaline tolerance in soybean. We identified a total of 28 SNPs significantly linked with alkaline tolerance-related traits. In agreement to our findings, the earlier studies have documented the similar results for salinity/alkaline tolerance in soybean [5,10,33].

The identification of candidate genes followed by their proper validation of function is the extreme objective of the researchers [36,37]. Previously, efforts to detect the genes for alkaline tolerance in soybean has been negligible. In our study, an in silico analysis detected 17 candidate genes underlying five QTLs viz., of qAT1, qAT4, qAT14, qAT18 and qAT20, and these genes were defined as possible candidates modulating alkaline tolerance in soybean. Based on the available information the gene functions related to salt tolerance, abiotic stress tolerance, protein phosphorylation, heat shock proteins, basic helix-loop-helix (bHLH) transcription factors, S-ribonuclease binding protein 1, GRAS family transcription factor, pentatricopeptide repeat (PPR-like) superfamily protein, C2H2-like zinc finger protein, abiotic stress tolerance, auxin signaling, kinase activity, ethylene signaling and salicylic acid signaling directly or indirectly govern related to alkaline tolerance. Hence, among the 17 identified genes, the functions of six genes, viz., Glyma.01G113400, Glyma.14G084500, Glyma.04G252100, Glyma.20G072600, Glyma.20G072700 and Glyma.20G072900 were found to be related to salt stress (Table 2). Alkali stress is often associated with salt stress, and alkaline stress is caused by excess of alkaline salts (NaHCO₃, Na₂CO₃) [38,39]. The Glyma.04G252300 function is related to auxin signaling, and studies in Arabidopsis showed that under saline-alkaline stress, the increase in the transcript level of AUX1 and auxin biosynthesis-related genes promoted auxin accumulation, and relieved the inhibition of root elongation under saline-alkaline stress [40,41]. The Glyma.04G253000 function is involved in ethylene signaling, and it has been demonstrated that H⁺-ATPase is involved in ethylene-mediated inhibition of rice growth under alkaline stress. Role of the protein phosphorylation and kinase activity in the alkaline stress response has been well documented [41,42], and the function of genes viz., Glyma.18G150300, Glyma.04G252700, Glyma.04G253100 and Glyma.20G072500 has been related to protein phosphorylation and kinase activity. The Glyma.14G083700 function as the heat shock protein (HSP), and the peroxisome-localized small heat shock protein Hsp17.6CII activates catalase by interacting with CAT2 and thus confers alkaline and salt stress tolerance in plants [43]. The Glyma.04G251900 is a GRAS family protein, and the GRAS protein has been reported to regulate salt tolerance in Melilotus albus [8] and kiwifruit [44]. The Glyma.14G083900 functions as E3 ubiquitin ligase, and the plant E3 ligases are involved in abiotic stress tolerance that negatively regulates salt/alkaline stress response [45]. The Glyma.04G252500 function as Pentatricopeptide repeat protein, and this protein has been reported in the Arabidopsis to regulate the salt/alkaline stress [25]. Thus, utilization of the above genes in the soybean breeding requires the verification of their function first either through overexpression or knockout experiments; after successful functional validation these genes can be used in the soybean breeding. Besides, we also identified haplotype alleles varying from two-four underlying the four candidate genes viz., Glyma.04G252300, Glyma.04G253100, Glyma.14G083700 and Glyma.20G072500. These haplotypes alleles govern the different phenotypes of alkaline tolerance ranging from lowest to highest through moderate levels, thus providing an opportunity to improve soybean alkaline tolerance. The use of these haplotypes in soybean breeding will generate alkaline-tolerant cultivars of soybean, and these cultivars will substantially influence the soybean productivity in the soils affected by alkaline stress.

4. Materials and Methods

4.1. Plant Materials and Experimental Design

In the present study, the plant material used for the GWAS analysis are a set of 326 diverse soybean germplasms collected from different agro-ecological regions of China (Table S3). To optimize the alkaline stress used for the evaluation of the GWAS panel, we selected two soybean cultivars viz., DS118 (alkaline-tolerant cultivar) and JY 86 (alkaline-sensitive cultivar) that are subjected to different levels of alkaline stresses [NaHCO₃: Na₂CO₃ in a 5:1 ratio (pH 9.0 ± 0.1)] viz., 50 nmol/L, 100 nmol/L, 150 nmol/L and 200 nmol/L in the seedlings that are two-week old after germination.

The soybean germplasms were cultivated in a controlled conditions using vermiculite for pot culture (8 cm × 8 cm × 8 cm). The growth conditions were set to a photoperiod of 16 light hours and 8 dark hours; the temperature of the light-phase and dark-phase is 28°C and 22°C, respectively, and a relative humidity was maintained at 80%. Uniform, full, and intact soybean seeds were selected for planting. A total of 18 pots were used for each soybean cultivar, with one seed sown per pot. After one week of germination, seedlings with similar growth vigor were selected for alkaline stress treatment (AT) using a solution of Murashige and Skoog (MS) medium supplemented with NaHCO₃: Na₂CO₃ in a 5:1 ratio (pH 9.0 ± 0.1). However, standard MS medium (pH 7.0 ± 0.1) was used as the control (CK). For each cultivar, the three pot plants under both the CK and AT were used. For each treatment, three liters of the corresponding solution were added to a fixed tray in which the pots are put, and the solution was replaced every three days. After seven days, the seedlings of both the CK and AT were used for the collection of the phenotypic parameters. Statistical design of complete randomized design (CRD) with three replications were used for the planting of germplasm.

4.2. Phenotypic Data Analysis

In the present study, we collected the phenotypic information of eight alkaline tolerance-related traits viz., seedling fresh weight (SFW), seedling dry weight (SDW), root dry weight (RDW), root fresh weight (RFW), relative chlorophyll content (CC), root number (RN), root tips number (RTN) and root length (RL) under both CK and AT. The RN, RTN and RL were measured using the Intelligence Root Analysis System (Zhejiang Top Cloud-Agri root analyzer GXY-B, China). The CC was estimated for the terminal leaflet of the first trifoliate leaf of soybean seedlings using the Konica Minolta SPAD-502Plus Portable Chlorophyll Meter (Soil and Plant Analyzer Development, SPAD, Japan). Measurements for CC were taken at the top, middle, and bottom of the leaflet, and the average of these three values was considered the CC of the leaf. For the measurement of fresh weights viz., SFW and RFW, the seedlings at the seventh day of AT were cut at the cotyledonary node to separate the shoot and root tissues. A precision electric balance (Model JA203H, Changzhou Xingyun Electronic Equipment Co., Ltd., China) was used to estimate the RFW and SFW for both CK and AT samples. These samples were then dried in an oven at 65°C for five days to obtain the RDW and SDW using the precision electric balance (Model JA203H, Changzhou Xingyun Electronic Equipment Co., Ltd., China).

The individual environmental data control (CK), alkaline stress (AT) and AT/CK (ACR) were used for the estimation of the CE data, and the “lme4” package was used for the estimation of the CE data following Di et al. [46]. The predicted means (BLUPs) were calculated by following the method of Yu et al. [1]. Descriptive statistical analysis parameters for alkaline tolerance-related traits were estimated according to Yu et al. [3]. The analysis of variance (ANOVA) were estimated following Yu et al. [1].

4.3. Genotyping, Population Structure and Linkage Disequilibrium (LD)

The association mapping panel consisting of 326 accessions were re-sequenced using the Illumina HiSeq platform (Novogene Company, Ltd., China), at an average depth of 25× with a total data volume of 8.57T. The adapter sequences, low-quality reads (where bases with a quality value Q≤20 accounted for over 50% of the read), and reads with >10% N content were removed by using the fastp software [47], and only the clean reads were retained. Alignment of the sequencing data with the reference genome Wm82.a4.v1 were performed using BWA [48]. Variants including SNPs, insertions/deletions (InDels), structural variations (SVs), and copy number variations (CNVs) were identified and marked using the Genome Analysis Toolkit (GATK) [49]. Annotation of the variants were performed by comparing each sample to the reference genome, with non-standard variants removed by using the quality control measures such as QD < 4.0, FS > 50.0, MQ < 40.0, and SOR > 3.0, resulting the 25,442,248 genetic variants. By using the further quality control measures viz., MAF<0.05 and a missing genotype at 0.01 to filtering out low-quality SNPs, we obtained a total of 3,311,166 high-quality SNPs.

To estimate genome-wide linkage disequilibrium (LD) the PopLDdecay software [50] was utilized to calculate the squared correlation coefficient (r²) between SNPs at known genomic positions. This software was also employed to estimate the expected r² under drift equilibrium, and to depict the same results over physical distance (kb). LD decay curves were fitted to scatter plots at the genomic level using smoothing spline regression lines [51], thereby illustrating the relationship between LD and physical distance.

4.4. GWAS Analysis

The GAPIT3 package [52] was utilized to conduct the GWAS. Seven GWAS models viz., CML [53], MLM [54], CMLM [55], SUPER [56], FarmCPU [57], MLMM and BLINK [58] fitted in the GAPIT3 package—were used in the present study. The CLM, MLM and CMLM are the single-locus model, whereas the remaining four models are multilocus models. To correct the population structure, the optimal number of principal component analyses (PCAs) was estimated using PLINK [59]. The visualization of the population structure was performed in the R environment using the “ggplot2” package. Significant SNP associations were determined utilizing the default parameters in GAPIT3 following Yu et al. [1].

4.5. Candidate Genes Identification

For candidate genes identification, we downloaded all the model genes within the physical interval of five major QTLs from SoyBase (https://www.soybase.org/) using the Williams 82 (Wm82.a4.v1) gene model. The annotations of these genes were also downloaded from SoyBase (https://www.soybase.org/). Based on the gene function annotation and literature search, the candidate genes underlying the physical interval of ±71.6 kb of the stable QTLs were selected [3].

4.6. Haplotype Analysis

Haploview 4.2 software were used to calculate the level of LD among the SNPs pairs [60]. The closest SNPs within the ±71.6 kb genomic interval form the haplotype block [4]. The “confidence intervals” algorithm were used to define the haplotype blocks [61]. The haplotype allele analysis and their phenotypic effect analysis was carried out by following the method as described by Yu et al. [1].

5. Conclusions

In this study, an integrated strategy of GWAS, QTL analysis, candidate gene and haplotype analysis was utilized to elucidate the genetic architecture of alkaline tolerance in soybean. Our study identified 28 significant SNPs and five QTLs viz., qAT1, qAT4, qAT14, qAT18 and qAT20 associated with alkaline tolerance. All the five QTLs identified were novel detected for the first time. The 17 genes lying within the physical interval of the five QTLs were regarded as putative candidates based on the in silico analysis. The results of haplotype analysis revealed that haplotype alleles of four genes viz., Glyma.04G252300, Glyma.04G253100, Glyma.14G083700 and Glyma.20G072500 showed significant difference in the regulation of alkaline tolerance. Hence, these four genes are regarded as the most probable candidate genes governing alkaline tolerance. Besides, two haplotypes blocks were identified on the Chr.04 and Chr.14. The haplotype alleles of the haplotypes identified in the present study varies from two to four regulating the alkaline tolerance from lowest to highest through moderate levels. Overall, the outcome of the present study such as significant SNPs, stable QTLs, candidate genes and haplotypes will be an efficient genetic resource for the enhancement of alkaline tolerance in soybean. The identified SNPs and QTLs will be further validated in the multiple genetic backgrounds for their direct utilization in the soybean breeding. Besides, the gene validation tests should be performed for the functional confirmation of the candidate genes before their direct deployment in the soybean improvement. In conclusion, the current study performed the detailed investigation of diverse soybean germplasms, and novel genetic/genomic resources for the improvement of soybean alkaline tolerance were identified.