To identify the GAPDH gene sequence of soybean, systematic BLASTP was conducted against the soybean reference genome database (https://www.soybase.org/) and the Phytozome database (https://phytozome-next.jgi.doe.gov/) using the published Arabidopsis GAPDH as alignment sequence. The screening threshold was set to E-value (<10⁻¹⁰), and the protein length was greater than 200 aa. The candidate GAPDH genes were determined by SMART (http://smart.embl-heidelberg.de) and Pfam (http://pfam.xfam.org/) software with both Gp_dh_N and Gp_dh_C domains. The open reading frame length was obtained from the Phytozome database. The molecular weight and the isoelectric point values were downloaded from the ExPASy (Artimo et al., 2012) software (https://web.expasy.org/protparam/). CELLO 2.5 (Yu et al., 2004) was used to predict subcellular localization.

The protein sequences of GAPDHs from soybean (Glycine max), maize (Zea mays), rice (O. sativa), and Arabidopsis (A. thaliana) were used to construct a phylogenetic tree using the neighbor-joining method and bootstrap test set at 1,000 replications through the MEGA7.0 software (Kumar et al., 2016). The exon/intron structures of GmGAPDHs were demonstrated at the GSDS online server (Hu et al., 2015). The coding and genomic sequences of GmGAPDH were collected from the Phytozome database. The conserved domains of GAPDH were determined by SMART (http://smart.embl-heidelberg.de) (Letunic et al., 2015) and Pfam (http://pfam.xfam.org/) (Finn et al., 2016) software, including Gp_dh_N and Gp_dh_C domains, and the structure of GAPDH proteins was visualized using the IBS 6.0 software (Liu et al., 2015).

To investigate the critical cis-acting elements in the promoter of GmGAPDH genes, the sequence at 2.0 kb upstream of the position of the ATG codon in these genes was obtained from the Phytozome database (https://phytozome-next.jgi.doe.gov/).The plant CARE database was used to predict the cis-acting regulatory elements, including motifs related to plant growth and development, plant hormone responses, and abiotic and biotic stress responses.

The expression patterns of soybean GmGAPDHs at different tissues were obtained using the Phytozome database. The heat maps were generated by cluster analysis with the TBtools software (Chen et al., 2020), and the expression data were log₂-transformed. To explore the expression patterns of GmGAPDHs in seeds at different developmental stages, soybean seeds (DN50) were collected at 10, 20, 30, and 40 days after flowering (DAF), and the total RNA extraction of each sample was performed to analyze the expression patterns of GmGAPDHs under abiotic stresses, including cold, salt, NaHCO₃, and drought. Briefly, soybeans (DN50) were grown in a plant incubator. There were two different plant cultivation methods used: (a) for low temperature treatment, seeds of soybean were sown in soil and vermiculite (v:v/1:1) and (b) for the salinity, NaHCO₃, and drought treatments, soybeans were grown in a hydroponic culture, and the growth condition was 24°C and a 16-h/8-h (day/night) daily photoperiod cycle. Second-trifoliolate-stage seedlings of uniform growth were subjected to cold treatment with a low temperature of 4°C, salt treatment with 150 mM salt, drought treatment with 20% polyethylene glycol (PEG, 6,000 g/M), and alkali treatment with 100 mM NaHCO_3. The soybean leaves were sampled at 0, 6, 12, and 24 h after the treatments. The sample total RNA was extracted using Trizol reagent (Invitrogen). The expressions of GmGAPDH genes in soybean seed samples at 10 DAF were used as a calibrator. The expressions of GmGAPDH genes in soybean samples at 0 h were used as a calibrator. GmACTIN4 (GenBank accession no. AF049106) was used as an internal reference. Quantitative real-time RT-PCR (qRT-PCR) was conducted using the CFX Connect TMreal-time system (BIO-RAD) with SYBR Select Master Mix RT-PCR (SYBR Green, TOYOBO, Osaka, Japan). Three biological replicates with three technical replicates were applied to each sample. The expression levels of GmGAPDHs were calculated using the 2^–ΔΔct method (Hong et al., 2010), and all primers used for the expression analysis were listed in Supplementary Table S1 .

Soybean cultivar DN50 was used for the Agrobacterium rhizogenes strain K599 transformation in soybean hypocotyls. The cDNA of GmGAPDH14 was directly ligated into the vector pCambia3300. The recombinant plasmid and empty pCambia3300 vector (EV) were transferred into Agrobacterium rhizogenes strain K599 and then injected into the hypocotyls following a previous report (Kereszt et al., 2007; Yu et al., 2021). The transgenic plants were identified by PCR amplification, and the non-transgenic hairy roots in the seedlings were removed.

The hairy root soybean plants were grown in Hoagland nutrient solution, and the growth chamber was set at a 16/8-h light–dark daily photoperiodic cycle. The transgenic plants were treated with 0 and 150 mM salt for 3 days (d), respectively.

The leaves of overexpression GmGAPDH14 (OE-GmGAPDH14) and EV seedlings were analyzed to measure physiological indicators. The measurements of superoxide dismutase (SOD) and malonaldehyde (MDA) were conducted according to corresponding assay kit protocols (Cominbio, Suzhou, China). All measurements were obtained with three biological replicates.

A total of 131 soybean germplasms were collected and grown in Harbin for 2 consecutive years (2019–2020). All accessions were treated with 0 and 150 mM salt. The relative swelling rate was obtained according to Zhang’s method (Zhang et al., 2014). According to genomic re-sequencing data, the SNPs in genomic regions including the promoter, 5′UTRs, exon, intron, and 3′UTRs of GmGAPDH14 gene were analyzed in 131 soybean lines using the generalized linear model (GLM) method as conducted with Tassel version 5.0 software (Bradbury et al., 2017).

The statistical significance was evaluated using Student’s t-test as performed with the SPSS 22.0 software. The significance levels were *p < 0.05 and **p < 0.01. The standard deviation (mean ± SD) was calculated with at least three biological replicates.

To identify the GAPDH gene family of soybean, GAPDH genes were selected by BLASTp tool with Arabidopsis as the alignment sequence. A total of 16 GAPDH genes (GmGAPDH1–GmGAPDH16) were retrieved ( Supplementary Table S2 ). The full-length CDS sequences of GmGAPDH1–GmGAPDH16 varied from 786 to 1,362 bp. The isoelectric points of GmGAPDHs ranged from 6.54 to 8.71, and the molecular weight of GmGAPDHs ranged from 28.2 to 48.4 kDa ( Supplementary Table S2 ).

To understand the evolutionary relationship of GmGAPDHs in soybean, the amino acid sequences of GmGAPDHs from A. thaliana (seven), Z. mays (12), and O. sativa (seven) were obtained from the NCBI database (https://www.ncbi.nlm.nih.gov/), and a phylogenetic tree was built. The phylogenetic tree indicated that the GAPDH proteins of soybean were clearly divided into three clusters (I–III) ( Figure 1A ). Cluster I, consisting of GmGAPDH4–6, GmGAPDH8–10, and GmGAPDH14–15, corresponded with cytosol isoforms containing AtGAPC proteins. Cluster II, including GmGAPDH2 and GmGAPDH12, corresponded with plastid isoforms containing AtGAPCp1 and AtGAPCp2 proteins. Cluster III, covering GmGAPDH1/3/7/11/13/16, corresponded with AtGAPA1, AtGAPA2, and AtGAPB proteins ( Figure 1A ).

The structures of soybean GAPDH genes were characterized with the GSDS software. As shown in Figure 1 , cluster II (GmGAPDH2 and GmGAPDH12) had the largest number of exons, including 12 exons. The exon number of cluster I ranged between 10 (GmGAPDH10) and 12 (GmGAPDH14); the remaining GmGAPDHs contained nine exons. The exon number of cluster III ranged from five to nine, and only GmGAPDH3/7 had nine exons; the remaining GmGAPDHs had five exons ( Figure 1B ).

The conserved domains analysis for 16 GmGAPDHs indicated that the GmGAPDHs revealed a multiple-domain protein, including Gp_dh_N (PF00044) and Gp_dh_C (PF02800) domains ( Figure 1B ). The Gp_dh_N domain (INGFGRIGR) and Gp_dh_C (GAAKAV) sequences were identified as highly conserved in the GmGAPDHs ( Supplementary Figures S1 , S2 ). A conserved active site (PS00071: ASCTTNCL) was found in most GmGAPDHs, except for GmGAPDH1 and GmGAPDH13 ( Supplementary Figure S1 ). The similarity of the gene structure and conserved domains of soybean GAPDH genes implies that they have undergone gene duplication during evolution.

To obtain the cis-elements of GmGAPDHs, sequencing of 2,000 bp upstream of all GmGAPDHs gene was performed based on the PlantCARE software. As shown in Figure 2A , a total of 22 cis-elements were found with plant growth and development, phytohormone-responsive, and abiotic and biotic stresses in the upstreams of 16 GmGAPDH genes. The ERE and ARE cis-elements were found with almost all GmGAPDH genes. The GCN4_motif (endosperm expression) elements were discovered in GmGAPDH6-7 and GmGAPDH9-11. Meanwhile, GmGAPDH15 and GmGAPDH14 were analyzed only involving AuxRR-core (auxin-responsive) and GC-motif (anoxic-specific inducibility) elements, respectively. It was noteworthy that there were five GmGAPDHs that harbored low temperature responsiveness (LTR) elements while five GmGAPDHs contained MBS (drought-responsive) elements. In addition, O₂ site (gliadin metabolic regulatory), CCGTCC box (specific activation), and CAT box (meristem expression) were found in the GmGAPDHs gene. TC-rich repeats (defense and stress responsiveness), GC motif (involved in anoxic-specific inducibility), and WUN motif (mechanical injury response) elements were observed in six, one, and 10 GmGAPDH genes, respectively ( Figure 2B ). These results showed that the GAPDH family may play an important role in growth and development and response to environmental stress in soybean.

To further characterize duplicated events within the soybean genome, a synteny analysis of GmGAPDH genes was performed. As shown in Figure 3A , the GmGAPDH genes were scattered on nine of the 20 soybean chromosomes. The nine soybean chromosomes distributed one to three GmGAPDH genes ( Figure 3A ). The replication relationship with the soybean GAPDH genes was analyzed. A total of 10 duplicated gene pairs were identified within the soybean ( Figure 3A ). Meanwhile, a synteny analysis was conducted from G. max and A. thaliana. As shown in Figure 3B , GmGAPDHs had a replication relationship with AtGAPDHs, including nine replication relationship pairs between G. max and A. thaliana (GmGAPDH1/AtGAPC1, GmGAPDH10/AtGAPC2, GmGAPDH11/AtGAPA2, GmGAPDH11/AtGAPA1, GmGAPDH11/AtGAPA2, GAPDH13/AtGAPC1, GAPDH14/AtGAPC2, GAPDH16/AtGAPA1, and GAPDH16/AtGAPC1) ( Figure 3B ). The duplicated genes showed their common genomic origin and maybe functional similarity.

To determine the expression pattern of soybean GAPDH genes in different development phases, we retrieved the high-throughput sequencing data of the Phytozome database and conducted an expression analysis. As demonstrated in Figure 4 , the expression of GmGAPDHs was revealed in diverse tissues. GmGAPDH4, 5, 8, 9, and 14 were found to have a higher expression in different tissues. Meanwhile, the GmGAPDH1, GmGAPDH6, and GmGAPDH15 genes were feebly expressed in nine different tissues. Moreover, GmGAPDH3, 7, 11, and 16 were expressed only in the shoot apical meristem and leaf tissues. GmGAPDH2 and GmGAPDH12 were especially expressed only in nodules ( Figure 4A ). Moreover, the expression of GmGAPDHs during the development of soybean seed is shown in Figure 4B (10 to 40 DAF). The expression levels of GmGAPDH4, GmGAPDH5, GmGAPDH9, and GmGAPDH11 were found to be upregulated in seeds at 20 DAF. The expression of the GmGAPDH8 and GmGAPDH14 genes was significantly upregulated in seeds at 30 DAF.

To confirm the role of GmGAPDH gene responses to various abiotic stresses, the transcription level of 16 GmGAPDH genes under NaHCO₃ (100 mM), PEG (20%), Cold (4°C), and Slat (150 mM) stresses was determined by qRT-PCR. Under simulated alkali stress using NaHCO₃, GmGAPDH 4, 5, 8, 10, 12, and 14 were extraordinarily upregulated (more than six folds) and peaked at 6 and 12 h, respectively ( Figure 5A ). After simulated drought stress using PEG, GmGAPDH14 was upregulated by more than 20 folds in 24 h after the treatment. In comparison, the expression of GmGAPDH16 significantly decreased in 6, 12, and 24 h, respectively ( Figure 5B ). In response to cold treatment, most of GmGAPDHs had upregulated expression, especially GmGAPDH4, which had a significantly higher level of expression at 24 h ( Figure 5C ). Furthermore, the GmGAPDH genes were shown to have a higher expression level at different timepoints of cold stress. The expression levels of GmGAPDHs were different under simulated salt stress with salt ( Figure 5D ). The GmGAPDH14 gene was shown to have the highest expression level at 6, 12, and 24, which was upregulated by more than 150 folds than that at 24 h. Meanwhile, the GmGAPDH4 gene exhibited a higher expression level at 6 and 12 h (20 and 40 folds). In contrast to other stresses, the whole expression of GmGAPDHs under NaHCO₃ stress was found to be relatively low ( Figure 5A ). Remarkably, the expression of GmGAPDH14 was sharply induced under salt stress, indicating that this gene might play a role in salt stress resistance.

To deeply illustrate how GmGAPDH14 genes respond to salt stress, hairy roots with Agrobacterium rhizogenes K599 containing pCambia3300-GmGAPDH14 plasmid or the pCambia3300 empty vector were transformed. Eight soybeans were proven through the PCR method to be positive transgenic. Furthermore, it was found that the effect of overexpression of GmGAPDH14 is such that it can regulate soybean hair roots in response to salt tolerance. Furthermore, 2-week soybean positive lines, involving 0 or 150 mM salt for 3 days, were transferred.

Previous studies showed that plant GAPDHs are involved in functions such as response to oxidative stress (Guo et al., 2012). Hence, the overexpression of GmGAPDH14 may further increase the antioxidant level under salt stress in this study. Therefore, the SOD activity and the MDA content were tested in soybeans positive transgenic at 3 d after 150 mM salt treatment. After 3 d of 150 mM salt stress, the overexpression of GmGAPDH14 was exhibited to enhance the resistance to salt stress than those of the control (transformed by pCambia3300 EV) ( Figures 6A, B ). As shown in Figure 6 , the OE-GmGAPDH14 lines had a higher activity of SOD than that of the EV lines ( Figure 6C ). The above-mentioned results showed that GmGAPDH14 could participate in regulating the ROS level. The MDA content showed that the contents of the OE-GmGAPDH14 lines were significantly lower than that of the EV lines after the salt treatment ( Figure 6D ).

To further confirm the potential effects of GmGAPDH14 gene for salt stress, gene-based association analysis was applied through the GLM method. A total of five SNPs in the GmGAPDH14 gene were identified among 131 lines ( Supplementary Table S3 ). All the five SNPs were significantly associated with salt stress, and they were located in the exon, intron, UTR region, and upstream regions of the GmGAPDH14 gene, respectively ( Table 1 ). Four haplotypes of the GmGAPDH14 gene were defined by the five SNPs. Haplotypes 2, 3, and 4 were composed of the combination of TAG, TTT, and AAT alleles and were beneficial for the salt tolerance of soybean. The carriers of haplotype 1 were composed of a combination of AAG alleles which tended to be sensitive to salt stress ( Figure 7 , Supplementary Figure S3 ). The difference of salt tolerance between the carriers of two haplotypes reached a very significant level.