Genome-Wide Identification and Abiotic Stress Response Analysis of the Isopentenyl Transferase (IPT) Gene Family in Soybean (Glycine max L.)

TL;DR

This study identifies and analyzes soybean IPT genes, revealing their roles in abiotic stress response and providing insights for soybean breeding.

Contribution

The study provides genome-wide identification and functional analysis of the soybean IPT gene family, including stress response and haplotype associations.

Findings

15 GmIPT genes were identified and characterized in soybean, distributed across 12 chromosomes.

GmIPT7, GmIPT10, and GmIPT15 show differential responses to drought, salt, and low-temperature stress.

A putative elite haplotype of GmIPT15 is associated with cold tolerance based on germination index.

Abstract

Isopentenyltransferase (IPT) is the rate-limiting enzyme in cytokinin biosynthesis and plays a critical role in plant acclimation to abiotic stress. To explore soybean IPT genes, we performed genome-wide identification, bioinformatics analysis, and molecular experimental validation to systematically characterize the features and functions of the soybean IPT (GmIPT) gene family. We identified 15 GmIPT genes in the soybean genome, which are unevenly distributed across 12 chromosomes; their evolutionary expansion is primarily driven by whole-genome duplication events. Phylogenetic analysis of soybean IPT proteins with those from Arabidopsis, rice and maize clustered them into four groups, exhibiting lineage-specific functional specialization. GmIPT genes exhibit significant variations in conserved motifs, gene structure, and cis-acting elements; their promoter regions are enriched in light-responsive, abiotic stress-responsive, and hormone-responsive elements, indicating their involvement in complex transcriptional regulatory networks. Tissue expression profiling revealed that GmIPT7 and GmIPT10 are highly expressed in various tissues, whereas GmIPT14 shows specific expression in flowers and the shoot apical meristem. Transcriptomic analysis and qRT-PCR validation demonstrated that GmIPT7, GmIPT10 and GmIPT15 respond differentially to drought, salt and low-temperature stress, with GmIPT15 exhibiting a transient upregulation at 3 h (p < 0.01) followed by a gradual decline to levels close to the pre-treatment control at 6–12 h under low-temperature stress. We further performed haplotype analysis of GmIPT15 and identified a putative elite haplotype (hap1) associated with cold tolerance based on low-temperature germination index assessment. This study provides useful insights for the future functional characterization of plant IPT genes and offers potential genetic resources and molecular markers that may support molecular-assisted breeding for soybean abiotic stress tolerance.

Figures

Full text

1. Introduction

Soybean is a vital source of vegetable oil and protein in global agricultural production, providing rich nutritional supply for humans and livestock. Abiotic stressors like extreme temperatures, salinity, and droughts have very detrimental effects on the soybean’s yield [1,2,3,4,5]. Breeding high-quality and stress resistant soybean varieties with high yield has been a fundamental goal in breeding soybean despite the challenges of population increase, better living conditions, and alterations in the environment.

Cytokinins (CTKs) are essential plant hormones that play critical roles in regulating growth, development, and environmental adaptation, and are involved in cell division, shoot apical meristem maintenance, delaying leaf senescence, nutrient signal transduction, and stress responses [6,7,8,9,10,11]. Plant growth status can be improved through optimization of cytokinin biosynthesis genes, which improve their content and distribution, achieve better performance, enhance regenerative capacity, and increase tolerance to abiotic stress [12,13,14,15,16].

Isopentenyl transferases (IPTs) serve as the rate-limiting enzyme in cytokinin biosynthesis. Their core function is to catalyze the transfer of the isopentenyl group from dimethylallyl diphosphate (DMAPP) to the N^6^ position of adenine nucleotides, generating isopentenyl adenine nucleotides (iPMP, iPDP, iPTP). IPT activity directly determines the efficiency and type of cytokinin synthesis, thereby participating in plant responses to abiotic stress [6,14,16,17,18,19,20]. Therefore, genome-wide identification and functional characterization of the IPT family represent a central approach to dissecting cytokinin regulatory networks and identifying target genes for crop genetic improvement and have become a research hotspot in the field of plant molecular breeding [21,22,23,24,25,26,27,28,29,30].

In Arabidopsis, AtIPT family members regulate CK synthesis through the ATP/ADP pathway and the tRNA pathway; the quadruple mutant of AtIPT1/3/5/7 shows significantly enhanced tolerance to salt and drought stress [20]. In rice, the 10 OsIPT genes are divided into four subfamilies, and their tissue-specific expression and stress response patterns are directly associated with grain yield and drought tolerance [21]. In cucumber, eight CsIPT genes are highly expressed in roots and male flowers and exhibit differential expression under low-light and salt stress, with stress-responsive cis-acting elements in the promoter regions serving as the core regulators [30]. In oilseed rape, 26 BnIPT genes respond to nitrogen deficiency through upregulation in shoots and downregulation in roots, participating in nitrogen transport and remobilization [24]. Ectopic expression of apple MdIPT5b can stabilize redox balance by accumulating proline, thereby significantly improving cold resistance [22]. In wheat, TaIPT2/7/8 are rapidly induced under drought within 0.5–1 h; TaIPT8 enhances drought tolerance by regulating tZ-type CK synthesis, and its drought-inducible overexpression improves tolerance without affecting plant growth, whereas mutants show reduced drought tolerance and accumulation of reactive oxygen species [28]. These studies indicate that the IPT gene family plays a central role in stress adaptation in different plants through gene structure divergence, expression pattern specialization, and crosstalk with hormone signaling, although functional mechanisms exhibit clear species-specific differences. However, detailed genome-wide characterization and functional differentiation of the IPT gene family in soybean, especially regarding natural variation and abiotic stress tolerance, remain largely unexplored. Most existing studies focus on general genome-wide gene identification, while haplotype variation and stress-associated allele mining of GmIPT genes have not been systematically investigated. This study aims to fill this knowledge gap and provide a targeted foundation for understanding the functions of IPT genes in soybean abiotic stress tolerance.

Based on the above, the present study used soybean as the material and combined bioinformatics and molecular biology approaches to conduct the following investigations: (1) Systematic identification of soybean IPT family members and analysis of their gene structure, conserved domains, chromosomal localization, and evolutionary relationships, to establish a basic genomic profile of the soybean IPT gene family. (2) Prediction of cis-acting elements in promoter regions to reveal potential transcriptional regulatory networks, thus facilitating understanding of their stress-related transcriptional regulation. (3) Screening of candidate genes using transcriptome data, combined with qRT-PCR validation of the expression patterns of GmIPT genes under drought, salt and low-temperature abiotic stresses, to identify core stress-responsive genes for subsequent functional research. (4) Analysis of haplotype variation in the core low-temperature-responsive gene GmIPT15 in natural soybean populations and mining of superior alleles associated with cold tolerance, providing valuable genetic targets for molecular breeding of soybean cold tolerance. The results of this study are expected to enrich research on the evolution and function of the plant IPT gene family and to provide important genetic resources, molecular markers, and technical support for improving soybean stress tolerance and breeding high-yielding, stable-yielding new varieties through molecular design breeding.

2. Results

2.1. Identification, Phylogenetic Analysis, and Chromosomal Localization of Soybean IPT Genes

To better and systematically search and find members of the soybean IPT gene family, we performed a genome-wide analysis. The searches were conducted by the use of nine Arabidopsis IPT protein sequences as a reference used in BLASTP searches. The Pfam database and conserved domain databases (CDD and SMART) also confirmed candidate genes to ascertain the correct identification of functional domains, and all of them had the IPPT domain (Table S1). The genome of the soybean identified 15 high-confidence GmIPT genes (Tables S1 and S2).

These GmIPT proteins were characterized in detail, showing considerable variation in their physicochemical properties (Table S3). The length of the amino acids was 206 to 478 residues with a molecular weight of 24.06 to 54.10 kDa. Theoretical isoelectric points (pI) ranged between 5.39 and 9.36 and indicated that there might have been some functional diversity in the various cellular environments. All proteins were found to be stable, with instability indices of 24.63–60.94 and aliphatic indices of 74.94–105.62 that showed thermal stability. A hydrophobicity analysis was performed and revealed that all GmIPT proteins were hydrophilic (GRAVY < 0). The subcellular localization of each GmIPT protein was accurately predicted as follows: GmIPT2, GmIPT3, GmIPT5, GmIPT9, GmIPT11, and GmIPT15 were localized in the chloroplast; GmIPT1 and GmIPT6 were localized in the chloroplast and cytoplasm; GmIPT4 was localized in the chloroplast, cytoplasm, and mitochondrion; GmIPT10 was localized in the chloroplast, cytoplasm, and peroxisome; GmIPT12 and GmIPT14 were localized in the chloroplast and mitochondrion; GmIPT7 was localized in the chloroplast, cytoplasm, mitochondrion, and peroxisome; GmIPT13 was localized in the peroxisome; and GmIPT8 was confined to the nucleus. These physicochemical properties and subcellular localizations suggest the development of functional specialization among GmIPT family members across distinct cellular organelles.

As an additional step to explain the evolutionary connections between IPT proteins, we cited the findings of the work on Arabidopsis IPT genes and carried out a comprehensive phylogenetic examination of the 44 protein sequences of four plant species: Arabidopsis, soybean, rice, and maize. The neighbor-joining tree which was created by multiple sequence alignments indicated clearly that the IPT proteins have four distinct clades (Figure 1A). Clade D was the most numerous branches with the majority of IPT members; clade A was the smallest cluster which included six conserved members. It is interesting to note that the clade B seemed to be Poaceae species specific with all the eight members of the clade having a rice and maize origin. The clade C had most of the conserved members of soybean, amounting up to 11. This phylogenetic distribution suggests that the evolutionary history of IPT genes in both functional conservation and functional divergence has occurred, with both Poaceae-specific clade B and the soybean-enriched clade C probably representing lineage-specific functional specialization.

Chromosomal localization studies showed that the 15 GmIPT genes noted were unevenly distributed in 12 chromosomes (Figure 1B) with most of the genes being distributed at the end of the chromosomes, with GmIPT1, GmIPT7, and GmIPT9 being in proximity to centromeric regions. Chromosomes Chr13 and Chr18 had two genes each, as compared to the other chromosomes, which had only a single gene. This is a non-random arrangement of the genomics, and this could possibly be the evidence of the evolutionary and functional specialization of the gene family.

2.2. Analysis of Gene Structure, Conserved Motifs, and Domains

To elucidate the structural characteristics and functional foundation of the GmIPT gene family, this study synthesized phylogenetic tree analysis, conserved motif occurrence, gene arrangement study, and an evaluation of domains to assess it extensively. The phylogenetic tree (Figure 2A) split the 15 GmIPT genes into three branches, in which the representatives of the various branches were dramatically different in their structural features. An analysis of conserved motifs (Figure 2B) revealed that most of the members had 9–10 motifs, but GmIPT10, with a very specific structure, only had motifs 2, 1, and 5, and GmIPT7 had an enrichment of motifs 8, 6, and 4 at the C-terminus leading to a more complex structure. This distinct motif variation across GmIPT members may imply the functional differentiation of this gene family, yet the specific biochemical functions of most conserved motifs have not been experimentally characterized and require further validation through functional experiments. The abundance of transcript lengths was evident, at the level of gene structure (Figure 2C), where the length of UTR and CDS varied significantly. Domain annotation (Figure 2D) revealed functionally relevant conserved structures. Pfam annotation indicated that the IPT domain was predominantly concentrated at the N-terminus in most members, while the IPPT domain appeared in the central or C-terminal regions of members such as GmIPT7 and GmIPT1. The ZnF_C2H2 domain was present in GmIPT7, GmIPT4, and others, whereas low complexity region and transmembrane domain were detected in GmIPT10, GmIPT14, and similar members. These findings demonstrate that during evolutionary processes, GmIPT family genes exhibit a pattern of structural diversification. The variations in the combinations of conserved motifs, gene structure and functional domains give the structural basis of their performance of specific biological functions.

2.3. Synteny Analysis of IPT Genes

Synteny relationships were compared within the soybean genome and between soybean and Arabidopsis, rice, and maize to thoroughly and systemically explain the evolutionary origins and chromosomal distribution properties of the soybean IPT gene family (Figure 3A). Intra-genomic synteny pattern revealed that most of the GmIPT genes are in 20 chromosomes with the formation of complex and dense collinear relationships among various chromosomes. In specific areas of chromosomes Gm03, Gm08, Gm10, Gm14, and Gm18, there was evident pairwise or even multi-pair collinearity. These findings suggest that the massive growth of the soybean IPT gene family followed extensive segmental duplications or whole-genome duplication, but not single tandem repeats. The wide spread of cross-chromosomal collinearity also indicates that soybean has experienced several events of genome repetition in the long evolution history. The process did not only add tremendously to the number of members of the IPT gene family but also could have facilitated the possible divergence in functionality.

Synteny relationships of IPT genes among Arabidopsis, maize, and rice were also examined (Figure 3B). The highest number of collinear genes was observed with Arabidopsis (15 genes), followed by rice (1 gene), while no collinearity was detected with maize. This suggests higher evolutionary conservation of IPT genes in dicotyledonous plants. However, compared with monocotyledonous plants, IPT genes retained homologous relationships in only a few chromosomal segments, indicating that significant structural rearrangements and functional divergence may have occurred after the divergence of monocots and dicots.

2.4. Cis-Acting Elements in the Promoter Regions of Soybean IPT Genes

Cis-acting elements play an important role in the transcriptional regulation of gene expression. To investigate the potential regulatory mechanisms of GmIPT genes in depth, we analyzed 2 kb sequences upstream of GmIPT genes using the Plant CARE tool and successfully identified 15 types of cis-acting elements. These putative elements include core promoters, enhancers, and silencers that may be involved in biotic and abiotic stress responses (Figure 4A; Table S4).

The findings demonstrated that the promoter regions of the GmIPT gene family contain several cis-regulatory elements such as promoters, enhancers and silencers among others with notable variation in the location and number of these elements among the different genes. Not only does this indicate the complexity of the expression regulation in this gene family, but it also suggests the nature of functional divergence.

Among them, light-responsive elements predominate in all GmIPT promoters, with nearly all showing high-density distribution, suggesting that IPT genes may potentially be widely involved in light signal-mediated regulation of growth and development processes. In addition, many putative defense and stress response elements, drought-responsive elements, low-temperature-responsive elements, and anaerobic/hypoxia-inducible elements were identified in most GmIPT promoters. The enrichment of these stress-related cis-elements in GmIPT promoters suggests that these genes may potentially contribute to the response to diverse abiotic stresses in soybean. Notably, there were significant differences in the types and numbers of cis elements among different genes. Within the GmIPT gene family, the promoters of GmIPT1, GmIPT6, GmIPT10, and GmIPT11 contained relatively large numbers of putative stress-related elements, whereas the overall number of elements in the GmIPT14 promoter was relatively low, suggesting that its transcriptional regulation may be more conserved.

Further statistical analysis of hormone-related cis-acting elements (Figure 4B) showed that GmIPT promoters generally contain multiple putative plant hormone response elements, such as methyl jasmonate (MeJA), abscisic acid (ABA), salicylic acid (SA), auxin (AUXIN), and gibberellin (GA) response elements. Among these, MeJA and ABA response elements were present in greater numbers in most GmIPT promoters and were particularly prominent in GmIPT2, GmIPT9, GmIPT10, GmIPT11, and GmIPT13, implying that these genes may potentially participate in hormone-mediated regulation of stress responses and defense reactions. In contrast, AUXIN and GA response elements were more limited in distribution, appearing only in certain genes, reflecting functional divergence of IPT family members in different hormone signaling pathways. Combined with the distribution patterns of environmental stress-related elements (Figure 4C, Table S4), it is inferred that the expression of soybean IPT genes is predicted to be not only broadly regulated by light signals but also co-regulated by multiple hormone and stress signals, thereby potentially exerting multi-level regulatory roles in cytokinin biosynthesis, plant growth and development, and abiotic stress adaptation.

2.5. Expression Analysis of GmIPT Genes in Different Tissues

To clarify the expression characteristics of the soybean GmIPT gene family across different tissues, transcriptome data were used in this study to systematically analyze the expression levels of GmIPT1-GmIPT15 genes in multiple soybean tissues, including roots, root hairs, nodules, stems, leaves, flowers, shoot apical meristems, pods, and seeds (Figure 5A,B), in order to compare and validate the expression features of the GmIPT gene family. Broadly, the expression pattern and difference across tissues indicated that various GmIPT genes differentiate differently, and diverse patterns of expression imply that the genes could play various biological roles during the growth and development of soybeans. Other genes had rather widespread patterns of expression; GmIPT7 and GmIPT10 were highly expressed in most vegetative and reproductive organs, with significant further upregulation in stems, leaves, flowers, and the apical meristem of the shoots. suggesting they may be associated with key cytokinin production pathways and aboveground growth control. Oppositely, GmIPT14 presented clear tissue-specific high-level expression in flowers and apical meristems of the shoot, whereas the expression levels were minimal in other tissues, and thus GmIPT14 may be related to reproductive growth or maintenance of meristem activity. Also, GmIPT9 was not highly expressed in roots and related tissues, while GmIPT3 was highly expressed in root hairs and nodules, and GmIPT12 was highly expressed in roots and nodules. In contrast, GmIPT1 and GmIPT5 showed some level of expression in nodules but not in roots or root hairs, whereas GmIPT6 was relatively expressed in roots or root hairs, which may indicate a potential association with development of the underground organs or related physiological events. Clustering analysis also found that the similar expression pattern of GmIPT genes could be able to cluster into a subclade (Figure 5B), indicating that this gene family has functional divergence and coordinated regulation at the level of tissue expression.

These findings provide important clues for subsequent investigations into the roles of GmIPT genes in soybean development and yield formation.

2.6. Expression Pattern Analysis of the Soybean IPT Gene Family Under Abiotic Stress

Transcriptome analysis revealed distinct global expression profiles of the entire GmIPT gene family in response to drought, salt, and cold stress, with clear differences in the overall responsiveness and dynamic patterns among the three abiotic stress treatments (Figure 6A–C). Collectively, the GmIPT family exhibited selective and divergent transcriptional responses to abiotic stresses: most members maintained low basal expression across all three stress treatments, while a subset of genes showed significant induction.

Exposure to 8% PEG8000-induced drought stress for 0 h, 6 h, and 12 h (Figure 6A) led to distinct transcript level changes in GmIPT genes: the expression of GmIPT1 and GmIPT10 was reduced at 6 h and restored at 12 h, while GmIPT15 displayed a continuous decline from 0 h to 12 h. These findings suggest that GmIPT1/GmIPT10 and GmIPT15 exhibit divergent expression patterns in response to drought stress. Notably, GmIPT10 maintained relatively high expression across all time points (0 h, 6 h, and 12 h), implying its potential involvement in the drought response. In contrast, the expression of most other GmIPT family members remained low and unchanged throughout the 0–12 h treatment.

For 0.9% NaCl-induced salt stress with treatment durations of 0 h, 1 h, 2 h, 4 h, 24 h, and 48 h (Figure 6B), the GmIPT gene family showed dynamic expression variations. GmIPT7 expression was upregulated at 1 h, dropped back to pretreatment levels at 2 h, and then increased continuously from 2 h to 48 h, peaking in the middle-to-late stages of the treatment. GmIPT10 was strongly induced during the early and middle phases of salt stress. By comparison, GmIPT1, GmIPT5, and GmIPT15 only had high expression at 0 h, with their transcript levels gradually decreasing as the treatment proceeded. All other GmIPT members kept low expression levels during the entire salt stress treatment.

Under 4 °C cold stress with treatment durations of 0 h, 1 h, and 24 h (Figure 6C), GmIPT1, GmIPT7, GmIPT10, and GmIPT15 responded rapidly at 1 h after treatment, among which GmIPT7, GmIPT10, and GmIPT15 showed more significant upregulation. At 24 h, the transcript levels of GmIPT1, GmIPT7, and GmIPT15 declined, whereas GmIPT10 maintained relatively high expression, suggesting distinct expression patterns. In addition, GmIPT5 showed no obvious change at 1 h and began to be upregulated at 24 h. Most other members of the GmIPT gene family remained at low expression levels throughout the treatment.

Based on the above expression profiles under the three abiotic stresses, cold stress, GmIPT7, GmIPT10, and GmIPT15 particularly exhibited obvious transcriptional responses to multiple stresses, suggesting their potential roles in plant stress responses.

In order to further clarify the transcriptional regulatory mechanisms of GmIPT7, GmIPT10 and GmIPT15 in soybean in response to major abiotic stresses, namely drought, salt, and low temperature, we determined the relative expression levels of the targeted genes using qRT-PCR at different time points (0 h as the control) following PEG-simulated drought (10% PEG6000), salt stress (100 mM NaCl), and low-temperature stress (10 °C) (Figure 7A–C). These findings revealed that these genes had stress-specific and gene-specific expression patterns. Under 10% PEG6000-simulated drought treatment (Figure 7A), GmIPT7 expression was upregulated and reached a significant level at 6 h post-treatment and declined hereafter. GmIPT10 showed a sharp increase in expression at 3 h post-treatment (p < 0.01) and continued to decrease with the extension of treatment time. In contrast, GmIPT15 exhibited sustained and significant downregulation starting from 3 h post-treatment (p < 0.01), and this inhibitory trend persisted until 48 h. Under 100 mM NaCl salt stress (Figure 7B), the expression of GmIPT7 showed no significant changes at all time points prior to 48 h, while its expression at 48 h was significantly higher than that of the 0 h control (p < 0.01); GmIPT10 was significantly downregulated at 3 h post-treatment (p < 0.01), with slight recovery at 6 h, and fell back to a low expression level at 48 h; GmIPT15 was significantly downregulated starting from 3 h post-treatment (p < 0.01), with a slight recovery only at 48 h but still remaining significantly lower than that of the 0 h control (p < 0.01). It should be noted that the expression patterns of GmIPT10 and GmIPT15 under salt stress showed some differences between the transcriptome profile (Figure 6B) and the qRT-PCR validation (Figure 7B). These variations are mainly attributed to differences in salt concentration, treatment time points, and detection methods used in the two experiments. At low temperatures (10 °C) (Figure 7C), The relative expression of GmIPT7 was significantly reduced (p < 0.01) at 3 h, 6 h and 12 h, maintaining a relatively low level throughout the stress period; GmIPT10 exhibited a pattern of rapid upregulation followed by sustained high expression: it was significantly upregulated at 3 h post-treatment (p < 0.01), reached its expression peak at 6 h, and maintained an significantly high expression level at 12 h (p < 0.01); GmIPT15 exhibited a unique transient upregulation pattern among the three genes: it was significantly upregulated at 3 h post-treatment (p < 0.01), then declined to a level close to the 0 h control during 6–12 h.

These findings confirmed that GmIPT7, GmIPT10, and GmIPT15 show divergent regulatory expression patterns when soybeans respond to drought, salt and low-temperature stresses, which implies that these genes mediate cytokinin metabolism to play divergent roles in soybean adaptation to abiotic stress.

2.7. Haplotype Analysis of GmIPT15 Associated with Low-Temperature Germination Index

Based on the expression characteristic of GmIPT15 that its expression level is transiently upregulated and gradually restored to the pre-treatment level under low-temperature stress, we speculate that this gene plays an important role in the low-temperature stress response of soybean. To further elucidate the genetic variation characteristics of this gene in natural soybean accessions and verify the differences in low-temperature stress resistance among its different alleles, we performed haplotype analysis on GmIPT15.

Sequence analysis of GmIPT15 in 1113 soybean accessions revealed one SNP locus and one InDel variation in the promoter region of this gene, which were located 279 bp and 366 bp upstream of the ATG start codon, respectively. Previous comparison of promoter elements showed that neither of these two variations occurred in the cis-acting element regions (Figure 4C and Figure 8A). In addition, one SNP locus was identified at position 507 in the exon region of GmIPT15, which caused a substitution of aspartic acid (Asp) with glutamic acid (Glu) at the 169th amino acid of its encoded protein, and this mutation site was in the functional IPPT domain (Figure 2D and Figure 8A). These three variation loci together formed five haplotypes in natural accessions, among which hap1 and hap5 were the dominant haplotypes, accounting for 98.38% of the total accessions combined.

To compare the differences in low-temperature resistance between these two dominant haplotypes, 54 accessions each carrying hap1 and hap5 were selected from natural soybean accessions, and their germination indices were determined after 15 days of low-temperature treatment at 8 °C. The results show that the germination index of hap1-type accessions was 0.56 ± 0.22, while that of hap5-type accessions was 0.45 ± 0.18. The germination index of hap1-type accessions was significantly higher than that of hap5-type accessions (p < 0.01), suggesting that accessions carrying hap1 may have stronger low temperature tolerance during germination (Figure 8B).

The five haplotypes exhibited significant genetic differentiation among distinct evolutionary populations (wild soybean, landrace, improved cultivar) (Figure 8C). Hap5 was the dominant haplotype in wild soybean, accounting for 68.75%; Hap1 was a rare haplotype in wild soybean, with a proportion of only 12.50%; while hap2-hap4 were wild soybean-specific haplotypes, each accounting for 2.50%, and were not detected in landraces or improved cultivars. It is inferred that these three haplotypes might have been eliminated by artificial selection during the domestication bottleneck. In landraces, the frequencies of hap1 and hap5 were almost equivalent, at 49.25% and 50.75%, respectively and combined to form the dominant haplotypes (Figure 8C). These results suggest that hap1, a haplotype putatively associated with enhanced low-temperature tolerance, may have undergone a considerable frequency increase during domestication from wild soybean to landraces, which could reflect positive selection for its genotype that is highly compatible with core soybean domestication traits. In the subsequent breeding stage, hap1 was retained at a relatively high proportion but with a slight frequency decline, which might imply partial loss of its selection advantage in modern breeding.

To investigate the geographical distribution of distinct haplotypes across the three major soybean cultivation regions in China (Huang-Huai-Hai, Northeast, and Southern Planting Regions) (Figure 8D), we statistically analyzed the geographic origins of soybean landraces and cultivars. For the geographic distribution patterns, hap1 showed the highest frequency in the Huang-Huai-Hai Planting Region (71.58%), followed by the Northeast Planting Region (56.00%) and the lowest in the Southern Planting Region (39.60%). In contrast, hap5 exhibited the highest frequency in the Southern Planting Region (60.40%), the lowest in the Huang-Huai-Hai Planting Region (28.42%), and a moderate level in the Northeast Planting Region (44.00%). Although this distribution pattern did not perfectly match the gradient of low-temperature stress intensity across the three regions (Northeast > Huang-Huai-Hai > South), it clearly indicated a strong correlation with the ecological adaptability of soybean. Specifically, hap1 was significantly more frequent in the two northern low-temperature regions (Huang-Huai-Hai and Northeast Planting Regions) than in the Southern Planting Region, where low-temperature stress is notably milder. This finding further confirms the selective effect of regional ecological adaptability on haplotype distribution. We hypothesize that the differences in low-temperature stress characteristics (e.g., stress duration and occurrence period) between the Huang-Huai-Hai and Northeast Planting Regions may underpin the stronger selective advantage of hap1 in the Huang-Huai-Hai Planting Region.

Combining the high-frequency distribution of hap1 in low-temperature planting regions with its significantly higher expression level compared to hap5, hap1 is suggested to represent a cold-adaptive haplotype. This haplotype may enhance soybean cold tolerance through sequence variation or differences in regulatory patterns, relying on its higher expression level, and has thus been continuously enriched under artificial selection in northern low-temperature planting regions. In contrast, hap5 shows relatively weaker cold-tolerance potential and declined in proportion under selection pressure in northern regions, but it is speculated to be potentially associated with core agronomic traits such as leaf shape development and photosynthetic efficiency, although no direct phenotypic data were obtained in this study to support this inference, thereby maintaining a relatively high retention frequency in improved cultivars. As the core soybean production area, the Northeast planting region requires breeding objectives to balance cold tolerance and yield traits; this demand ultimately led to the balanced distribution of hap1 and hap5, reflecting a compromise optimization strategy under artificial selection for dual breeding goals.

3. Discussion

Cytokinins serve as a core hormone regulating plant growth, development, and adaptation to abiotic stress, with their biosynthesis efficiency directly determined by genes of the isopentenyltransferase (IPT) family. The present study systematically characterizes the evolutionary features, expression patterns, and functional potential of the soybean IPT gene family using genome-wide identification, bioinformatics analysis, and molecular experimental validation. These findings provide new insights into cytokinin-mediated regulatory networks for abiotic stress tolerance in soybean and offer key targets and marker resources for molecular breeding.

3.1. Evolutionary and Structural Divergence of the GmIPT Gene Family: Foundation of Functional Specialization

The structural variation and growth of gene families are the most important forces in adaptive evolution of plants. In this study, we found 15 GmIPT genes in the genome of the soybean (Glycine max), which is also higher than those found in Arabidopsis thaliana (9 genes) [31], Oryza sativa (10 genes) [21,31], Zea mays (10 genes) [23], and most of the characterized Rosaceae species [22]. The growth is mainly due to whole-genome duplication (WGD) which is in line with the pattern of expansion of gene families that have been observed in soybean, a paleopolyploid plant [32]. It is also an indication of the adaptive evolution of IPT gene family in soybean to intricate environmental pressures.

Phylogenetic analysis grouped the IPT proteins of soybean, Arabidopsis, rice and maize into four clusters, with Clade B uniquely Poaceae and Clade C uniquely enriched in the soybean (11 members). Such clustering pattern by lineage also suggests that IPT genes have experienced a functional divergence over the long evolutionary periods. The same analyses have been done in IPT gene family analysis in Medicago truncatula [33], implying that the evolutionary divergence of the IPT genes is strictly connected to the species-specific environmental adaptation. Analysis of conserved motifs and gene structure indicated that there were major differences among the members of GmIPT. An example is that GmIPT10 has only 3 conserved motifs, GmIPT7 is enriched in specific motifs (8, 6, and 4), but it contains a Zn_F C2H2 domain, which is a structural domain that has been well documented to mediate specificity in the binding of DNA and protein–protein interaction in abiotic stress signaling [34]. Such structural variations explain functional divergence, which has a molecular foundation.

The analysis of promoter cis-acting elements demonstrated that all GmIPT promoters are enriched with light-responsive elements and most of them have abiotic stress-promotable elements (e.g., drought and low temperature) and hormone promotable elements (e.g., MeJA and ABA). This trend is similar to the findings of IPT gene family studies in many other plants, such as cucumber [30] and oilseed rape [24], where IPT genes are speculated to be involved in the coordinated regulation of light signals, hormone signals, and stress signals. It is noteworthy that GmIPT1 and GmIPT6 have more stress-related elements, which agrees with their possible high ability to respond to stress unlike GmIPT14 which has less elements indicating more conservation in transcriptional regulation. It is probable that these differences in cis-acting factors play a major role in the differences in the expression trends of GmIPT genes during stress.

3.2. Expression Patterns: Functional Differentiation in Adaptation to Soybean Development and Stress Tolerance

The biological functions are closely related to tissue-specific expression patterns of genes. Tissue expression patterns showed that GmIPT7 and GmIPT10 are expressed at significant levels in various tissues, such as stems, leaves, flowers, and apical meristems in the shoot-like tissues, as is the case with rice OsIPT2 and OsIPT9 [31]. This indicates a possibility of GmIPT7 and GmIPT10 being central members of cytokinin production in soybean, which are involved in the control of aboveground growth, meristem development, and reproductive development. Cytokinins have a central role in cell division and senescence inhibition [8,10,11], and the wide distribution of GmIPT7 and GmIPT10 provides a stable cytokinin concentration at different stages of soybean development.

Analysis of abiotic stress response is critical towards determining functional candidate genes. The transcriptome and qRT-PCR data showed that GmIPT7, GmIPT10 and GmIPT15 react differentially to drought, salt and low temperature stress, indicating the divergent regulatory mechanisms of the GmIPT gene family to different stresses. When under drought conditions, GmIPT7 is upregulated at 6h and then downregulated a similar response to the initial rapid response to drought of wheat TaIPT2/7/8 [28]. This implies that GmIPT7 might play a role in the early signal transduction in drought stress by regulating the production of cytokinins, which will in turn activate the downstream stress signals. Cytokinins have the potential to boost plant drought resistance through the augmentation of antioxidant ability and cell turgor [35,36,37] and the prompt induction of GmIPT7 could allow soybean to adjust to drought more rapidly. GmIPT10 is continuously expressed in response to low-temperature stress, and GmIPT15 is uniquely characterized by transient upregulation (upregulation at 3 h, followed by a gradual decline and restoration to basal levels at 12 h), a behavior that is normally indicative of rapid cold acclimation induction and preservation of cellular homeostasis. These expression patterns that are characterized by transient activation followed by gradual recovery to basal levels have been observed in genes of the Arabidopsis C-repeat binding factor (CBF) pathway, which function as molecular switches in the early stages of cold stress [38], though no experimental validation for this similarity was performed in the present study.

The stress response properties of the GmIPT genes have a close relationship with the cis-acting elements within promoters. As an illustration, GmIPT15 has promoter elements that are responsive to low temperature, which is in line with its high sensitivity to low temperature. In addition, cytokinin and abscisic acid (ABA) tend to display an antagonistic relationship in stress responses [37], and the presence of ABA-binding motifs in GmIPTs promoters suggests a potential involvement in crosstalk between cytokinin and ABA signaling pathways to control soybean tolerance of stress; however, this inference remains speculative and lacks direct experimental validation in the present study. These results validate the retained role of IPT genes in plant stress responses in addition to the species-specific response characteristics of soybean GmIPT genes.

3.3. Haplotype Analysis of GmIPT15: Elite Allele for Cold Tolerance and Breeding Value

Although previous studies have identified excellent cold-tolerant alleles related to soybean cold tolerance through haplotype analysis [39,40], existing research still has limitations, including limited coverage of germplasm resources, uneven coverage of cold-tolerant candidate gene families, and insufficient analysis of the domestication and evolution, geographical adaptability, and functional differentiation mechanisms of haplotypes. This study focused on GmIPT15 which contains three polymorphic sites: two are located in the promoter region (279 bp and 366 bp upstream of the start codon ATG, respectively), and one is a non-synonymous SNP in the exon region, leading to the substitution of aspartic acid (Asp) with glutamic acid (Glu) at position 169 of the encoded protein, and this site is located within the IPPT domain.

Systematic haplotype typing of the GmIPT15 gene was performed on 1113 soybean germplasms, and a total of five haplotypes were identified, among which hap1 and hap5 were the dominant haplotypes. An 8 °C low-temperature germination experiment (with germination index as the evaluation index) confirmed that hap1 is a cold-tolerant superior haplotype. The germination index of germplasms carrying hap1 (0.56 ± 0.22) was significantly higher than that of germplasms carrying hap5 (0.45 ± 0.18) (p < 0.01), providing a new molecular target for soybean cold-tolerant molecular breeding.

Analysis from the perspective of domestication and evolution showed that hap1 is a rare haplotype in wild soybeans (accounting for 12.50%), while its frequency surged to 49.25% in landraces, indicating that the GmIPT15 gene has undergone strong artificial selection during soybean domestication. As an excellent cold-tolerant haplotype, the retention and enrichment of hap1 are important genetic foundations for soybeans to adapt to the low-temperature environment in the Northeast Planting Region and Huang-Huai-Hai Planting Region, highlighting its core role in soybean cold-tolerance domestication. Geographical distribution characteristics showed that the distribution of hap1 presented obvious regional differences: it had a low frequency in the Southern Planting Region where low-temperature stress is mild but was highly enriched in the two major low-temperature soybean producing areas, the Huang-Huai-Hai Planting Region (71.58%) and Northeast Planting Region (56.00%). This pattern confirms that the cold-tolerant trait of hap1 is jointly driven by natural selection and artificial selection, and its distribution is highly compatible with the regional low-temperature stress environment. The frequency of hap1 in the Huang-Huai-Hai Planting Region is higher than that in the Northeast Planting Region, which is speculated to be closely related to the differences in low-temperature stress characteristics between the two major producing areas: the low-temperature duration in the Northeast Planting Region is long, but the frequency of low-temperature occurrence during the germination period is relatively stable; the spring low-temperature in the Huang-Huai-Hai Planting Region is volatile and sudden, and is often accompanied by complex adversities such as drought, which imposes higher requirements on soybean cold tolerance during the germination period. Therefore, hap1 is more likely to be selected and retained in the Huang-Huai-Hai Planting Region, forming a high-frequency distribution.

In summary, this study has achieved multiple extensions and improvements on the basis of previous studies: compared with the 324 soybean germplasms used by Zheng et al. [39], this study expanded the number of germplasms to 1113, covering a wider range of evolutionary populations (wild soybeans, landraces) and geographical populations (Southern Planting Region, Huang-Huai-Hai Planting Region, Northeast Planting Region), which improves the reliability and universality of the identification of the excellent haplotype hap1; it fills the gap in the research on cold-tolerant haplotypes of the IPT gene family, systematically analyzes the polymorphic sites and haplotype characteristics of the GmIPT15 gene, especially the role of the non-synonymous SNP within the IPPT domain; through systematic analysis of domestication, evolution and geographical distribution, it deepens the understanding of the adaptation rules of cold-tolerant haplotypes, and reasonably infers the distribution differences between the Huang-Huai-Hai Planting Region and Northeast Planting Region, providing more targeted theoretical support and molecular targets for soybean cold-tolerant molecular marker-assisted selection (MAS) breeding.

3.4. Limitations of the Study

Several limitations of the present study should be acknowledged. First, the functional roles of GmIPT genes in soybean abiotic stress responses were inferred primarily from expression pattern analyses and promoter cis-element predictions, without direct experimental validation using transgenic or mutant approaches. Although the differential expression of GmIPT7, GmIPT10, and GmIPT15 under drought, salt, and cold stress suggests their potential involvement in stress tolerance, the lack of gene knockout, overexpression, or complementation experiments means that the causal relationships between these genes and stress tolerance traits remain to be confirmed. Second, the association between the hap1 haplotype and cold tolerance was established based on germination index data and haplotype frequency analysis in different planting regions, without supporting evidence from physiological indicators (e.g., electrolyte leakage, proline content) or molecular marker validation in diverse genetic backgrounds. Similarly, the speculation that hap5 might be linked to core agronomic traits was not supported by direct phenotypic data in this study. Third, our study was conducted under controlled laboratory conditions, which may not fully reflect the complex and fluctuating environmental stresses encountered in field conditions. Further field trials are needed to validate the stress tolerance phenotypes associated with GmIPT genes and haplotypes in real agricultural settings. These limitations indicate that the conclusions drawn in this study are preliminary and based on correlative evidence, and further functional verification and field validation are required to fully elucidate the roles of GmIPT genes in soybean adaptation to abiotic stresses.

Overall, this study elucidates the evolutionary history and functional differentiation of the soybean GmIPT gene family, determines GmIPT15 and its elite haplotype hap1 as the basic candidates to growth under cold stress, and offers valuable genetic materials and technical assistance to abiotic stress-resistant breeding in soybean. The results will contribute to the understanding of IPT gene roles in plants and provide another example of a case study to understand how a gene family evolves and adapts to stresses in legumes.

4. Materials and Methods

4.1. Identification of IPT Genes in the Soybean Genome and Prediction of Protein Physicochemical Properties

To identify members of the IPT gene family in the soybean genome, the hidden Markov model file for IPT (Pfam number: PF01715) was downloaded from Pfam [41] and used for preliminary screening with the Advanced Hmmer search plugin in TBtools (v2.056) [42]. To balance sensitivity and specificity, a strict E-value threshold of ≤1 × 10^−5^ was applied, which is a widely accepted standard in gene family identification studies to minimize false positives while retaining bona fide family members. At the same time, protein sequences of Arabidopsis IPT gene family members were downloaded from the TAIR [43] website (https://www.arabidopsis.org/) as query sequences, and gene family members were identified using the BLAST plugin in TBtools. The union of results from both methods was obtained. To remove redundant sequences and ensure the uniqueness of identified genes, the resulting sequences were filtered to retain only non-redundant entries, with duplicate accessions and identical sequences excluded from the final dataset. Sequences with incomplete domains were screened and removed using Batch CD-Search [44] and the SMART [45] database. All physicochemical property parameters of the identified soybean IPT proteins, including amino acid number, molecular weight, isoelectric point, instability index, aliphatic index, and grand average of hydropathicity, were analyzed in detail using the Protein Parameter Calc plugin in TBtools; this process was consistent with the method for protein physicochemical property analysis using TBtools software as demonstrated by researchers in video tutorials. Subcellular localization prediction was performed using Cell-PLoc 2.0 [46]. All gene names and sequences are listed in Supplementary Table S2.

4.2. Analysis of Domains, Motifs, and Gene Structure

The MEME [47] online tool was used to analyze conserved motifs of IPT, with the motif number set to 10 and the remaining parameters using default settings. Subsequently, the Gene Structure View plugin in TBtools was employed to synchronously visualize gene structure, motifs, conserved domains, and the phylogenetic tree.

4.3. Construction of Phylogenetic Tree and Chromosomal Localization

Using default parameters in MEGA12 [48] software, IPT protein sequences from soybean, Arabidopsis, rice, and maize were aligned. The aligned sequences were then used to construct a phylogenetic tree by the neighbor-joining (NJ) method in MEGA12 software (Bootstrap value of 1000 replicates), and the tree was beautified using iTOL (v4) [49]. By analyzing the GFF file downloaded from the Glycine max Wm82.a4.v1 [50] genome database in Phytozome 14 (https://phytozome-next.jgi.doe.gov/), the physical positions of soybean IPT genes on chromosomes and chromosome lengths were determined. The IPT genes were mapped to 12 soybean chromosomes and visualized using the Gene Location Visualize from GTF/GFF plugin in TBtools software.

4.4. Promoter Analysis, Gene Duplication, and Synteny Analysis

Using the GXF state and Table Row Manipulate plugins in TBtools software, sequences 2000 bp upstream of the transcription start site of candidate soybean IPT family genes were extracted. The PlantCARE tool [51] was used to predict cis-acting elements in the promoter regions of soybean IPT gene family members based on these 2000 bp sequences.

With the aid of TBtools software, the evolutionary relationships of IPT genes within the genome and synteny among multiple genomes were investigated. The Advanced Circos plugin in TBtools software was used to construct a circular diagram of duplication events in the soybean genome. For inter-genomic synteny analysis, the One Step MCScanX plugin in TBtools was employed to analyze BLASTp search results between soybean and Arabidopsis, rice, and maize (using default parameters).

4.5. Expression Analysis Using RNA-Seq Data

Gene expression levels (FPKM, fragments per kilobase of exon model per million mapped reads) in nine soybean tissues were retrieved as pre-normalized expression values from the SoyBase database [52], and a heatmap of IPT gene expression across the nine different tissues was constructed using the HeatMap plugin in TBtools. Based on quality-controlled and normalized transcriptome data from the SoyOD database [53], the expression patterns of GmIPT genes under drought, salt stress, and low-temperature conditions were compiled.

4.6. qRT-PCR Validation Analysis of Abiotic Stress Expression Patterns of Soybean GmIPT Genes

Seeds of the soybean cultivar Williams 82 were selected, surface-disinfected with 75% ethanol for 30 s, and sown in a vermiculite-perlite (2:1) mixed substrate. Plants were grown under 25 °C, 16 h/8 h photoperiod, and 60–70% relative humidity until the V3 stage. Seedlings were then subjected to PEG 6000-simulated drought, 100 mM NaCl salt stress, and 10 °C low-temperature treatment, respectively. The second trifoliate leaves were collected at 0 h (control), 3 h, 6 h, 24 h, and 48 h (for drought and salt stress) or at 0 h, 3 h, 6 h, and 12 h (for low-temperature stress), flash-frozen in liquid nitrogen, and stored at −80 °C for later use.

Total RNA from leaves was extracted using TRIzol reagent, treated with DNase I to remove genomic DNA contamination, and assessed for RNA purity and integrity. First-strand cDNA was synthesized using 1 μg of purified RNA as template and diluted for use as qRT-PCR template. Gene-specific primers for GmIPT7, GmIPT10, and GmIPT15 were designed using the NCBI Primer-BLAST (https://www.ncbi.nlm.nih.gov/tools/primer-blast/, accessed on 11 January 2026) tool (Table S5), with GmACTIN6 (LOC100792119) as the reference gene. Primer specificity was verified by agarose gel electrophoresis and melting curve analysis. For qRT-PCR primer efficiency validation, standard curves were constructed using a 5× gradient dilution series of mixed cDNA samples. PCR amplification efficiency (E) was calculated according to the equation: E = (10^(−1/slope)^ − 1) × 100%. All primer pairs exhibited satisfactory amplification efficiencies ranging from 90% to 110%, and the correlation coefficients (R^2^) of all standard curves were greater than 0.99, indicating high reliability and specificity of the qRT-PCR results. qRT-PCR was performed using SYBR^®^ Premix Ex Taq™ II reagent (Takara Bio Inc., Kusatsu, Shiga, Japan) on a CFX96 Touch system, with the reaction system and program following the kit instructions.

Relative gene expression levels were calculated using the 2^−ΔΔCt^ method. One-way analysis of variance was performed using GraphPad Prism (v10.4.1) software, and Dunnett’s test was used to compare expression differences at each time point with the 0 h control group. There were three biological and three technical replicas. The results are reported in terms of mean ± standard error (Mean ± SE).

4.7. Haplotype Analysis of GmIPT15 in Soybean Germplasms

A total of 1113 soybean germplasms with different origins were used as test materials in this study, and the haplotype analysis of GmIPT15 was conducted on the basis of whole-genome resequencing data of all test germplasms. SNPs in the 1500 bp upstream regulatory region and coding region of GmIPT15 were first extracted from the whole-genome resequencing VCF file using BCFtools (v1.15) [54]. The extracted SNPs were then annotated using SnpEff (v5.1) [55] to identify nonsynonymous, nonsense, and silent mutations, and only nonsynonymous SNPs that could cause amino acid sequence changes were retained for subsequent analysis, while nonsense and silent mutations were excluded. Low-frequency variants with an allele frequency of less than 5% were removed using VCFtools (v0.1.16) [56] to reduce the interference of low-frequency variants on haplotype typing results and ensure the reliability and accuracy of subsequent analyses. Based on the retained high-confidence functional SNP loci after screening, genotyping of the target gene loci was carried out for the 1113 soybean germplasms, and haplotypes were constructed by concatenating the alleles of retained SNPs for each germplasm, with germplasms sharing identical genotypes at the target gene loci classified into the same haplotype group. Finally, combined with the type classification of the test soybean germplasms, the distribution quantity and proportion of different haplotypes in each type of soybean were counted using custom R scripts (R v4.2.1) [57] and Haploview (v4.2) [58] to clarify the haplotype composition characteristics and proportional differences among various types of soybeans.

4.8. Determination of the Germination Index for Distinct Haplotypes Under Low-Temperature Conditions

A total of 108 representative soybean (Glycine max (L.) Merr.) germplasms were selected based on the haplotype analysis of the target gene. Prior to the experiment, plump, disease- and pest-free seeds with uniform size were selected to ensure germination consistency. The seeds were surface disinfected by soaking in 1% sodium hypochlorite solution for 20 min, then rinsed 5–6 times with distilled water to remove residual disinfectant, and the surface moisture was blotted dry on sterile filter paper afterwards. A completely randomized block design was adopted with three biological replicates for each germplasm. For each replicate, 20 pretreated seeds were evenly placed in 9 cm Petri dishes lined with two layers of moist sterile sponges, and 12 mL of distilled water was added to keep the sponges moist. All Petri dishes were incubated in a constant-temperature incubator at 8 °C in the dark, and 2 mL of distilled water was supplemented every 2 days during incubation to maintain a moist environment for seed germination. Seed germination was defined as radicle protrusion through the seed coat by ≥2 mm, and germinated seed numbers per replicate were recorded daily for 15 consecutive days post-incubation, with mean values across replicates calculated daily. The germination index (GI) was calculated using the formula $[eqn]$ where $[eqn]$ denotes the mean number of germinated seeds on day t, and $[eqn]$ is the number of germination days. This index comprehensively reflects the germination speed and germination rate of soybean seeds, with a higher value indicating stronger low-temperature germination ability. The Shapiro–Wilk test [59,60] was used to verify the normality of the germination index data for the distinct haplotypes, and the results indicated that the data did not conform to a normal distribution. Thus, the Mann–Whitney U test (Wilcoxon rank-sum test) [59,60] was applied to determine the significant differences in germination index between the two haplotypes at the 0.05 probability level, and the results were presented as mean ± standard error (Mean ± SE).

5. Conclusions

In this study, we identified 15 GmIPT genes at the genome-wide level in soybean (Glycine max) and systematically characterized their chromosomal distribution, evolutionary expansion features, phylogenetic relationships, conserved motifs, gene structures, and cis-acting elements, while also investigating in depth their response patterns to abiotic stresses. Our results reveal that GmIPT15 exhibits a transient response to cold stress, thereby representing a potential candidate gene involved in cold stress adaptation in soybean. Subsequent haplotype analysis of GmIPT15 led to the identification of Hap1 as a haplotype associated with cold tolerance at the germination stage, with its domestication characteristics and geographical distribution patterns further elucidated. Collectively, this study advances our understanding of IPT genes’ characteristics in soybean and provides a potential genetic resource and preliminary reference for molecular marker-assisted breeding of cold-tolerant soybean varieties, pending further functional validation and field trials.