Skip to main content
Download PDF

The members of zinc finger-homeodomain (ZF-HD) transcription factors are associated with abiotic stresses in soybean: insights from genomics and expression analysis

BMC Plant Biology volume 25, Article number: 56 (2025) Cite this article

Abstract

Background

Zinc finger homeodomain (ZF-HD) belongs to the plant-specific transcription factor (TF) family and is widely involved in plant growth, development and stress responses. Despite their importance, a comprehensive identification and analysis of ZF-HD genes in the soybean (Glycine max) genome and their possible roles under abiotic stress remain unexplored.

Results

In this study, 51 ZF-HD genes were identified in the soybean genome that were unevenly distributed on 17 chromosomes. All GmZF-HD genes contained a conserved ZF-HD_dimer domain and had diverse physicochemical features. Furthermore, the GmZF-HD gene structures exhibited 3 to 10 conserved motifs, and most of them showed intronless gene structures. Phylogenetic analysis categorized them into eight major groups with the highest closeness to dicots including Brassica rapa and Malus domestica. The cis-element analysis recognized plant growth and development (10%), phytohormones (31%) and stress-responsive (59%) elements. Synteny analysis identified 73 segmental and 1 tandem duplicated genes that underwent purifying selection. The collinearity analysis revealed that GmZF-HD genes showed higher homology with dicot species, indicating common ancestors with close evolutionary relationships. A total of 94 gma-miRNAs from 41 diverse miRNA families were identified, targeting 40 GmZF-HD genes, with GmZF-HD6 being most targeted by 7 miRNAs, and gma-miR4993 emerging as the dominant miRNA family. Different TFs including ERF, LBD, BBR-BPC and MYB, etc., were predicted in all 51 GmZF-HD genes upstream regions and visualized in the network. Expression profiling through RNA-Seq showed diverse expressions of GmZF-HD genes in different tissues including seeds, roots, shoots and leaves under diverse conditions. Further, the qRT-PCR analysis demonstrated that all tested GmZF-HD genes were significantly induced in soybean leaves, mainly the GmZF-HD5/6/13/39 and GmZF-HD45 genes were significantly upregulated (2.5 to 8.8 folds) under the tested stress treatments compared to control, highlighting their potential roles in response to stresses in soybean.

Conclusion

Overall, this study reveals comprehensive insights into the ZF-HD genes in soybeans and provides a valuable contribution towards functional studies for soybean improvement under stress conditions.

Peer Review reports

Background

Soybean (Glycine max) is one of the major leguminous crops with significant economic value and a source of edible oil, protein, phytochemicals, and nutrients for both human consumption and animal feed [1]. Despite its importance, soybean cultivation faces significant environmental challenges from various biotic and abiotic stresses including extreme temperatures, drought, nutrient deficiencies, salinity, and heavy metals [2,3,4,5]. Excessive accumulation of metal ions such as cadmium (Cd2+), cobalt (Co2+), manganese (Mn2+), zinc (Zn2+), iron (Fe2+,3+), copper (Cu2+) and nickel (Ni2+) can lead to the production of reactive oxygen species (ROS), which negatively affect plant growth and development [6, 7]. In general, environmental stress adversely affects plant growth, and development and causes a significant threat to crop yields and food security [8,9,10]. Furthermore, under environmental stress, the plant immune system triggers molecular signaling pathways that protect plants from various stress damages. Transcription factors (TFs) are crucial in regulating stress responses by binding to the promoters of target genes and regulating their expression to tolerate multiple stresses [11]. Different TFs including bZIP, MYB, ERF, WRKY, bHLH and NAC have been reported to play crucial roles under various abiotic stresses including drought, cold, salt, hormones and heavy metal resistance [11,12,13,14]. TFs respond to adverse conditions by binding to promoters of target genes that regulate the secondary metabolites and produce stress-induced phytochemicals that play an integral role in plant immunity [15, 16]. Furthermore, TFs confer tolerance to heavy metals by regulating antioxidants and lipids [17]. Due to their important roles, TFs are promising candidates for enhancing crop resilience through genetic engineering.

Among plant TFs, zinc finger homeodomain (ZF-HD) genes are plant-specific transcription factors, that belong to zinc finger proteins and play key roles in plant growth, development, and various stress responses [18,19,20]. ZF-HD genes were first identified in the Flaveria as a potential regulator of photosynthesis gene encoding phosphoenolpyruvate carboxylase (PEP Case) [21]. ZF-HD genes contained two highly conserved domains including the C-terminal homeodomain (HD) domain and the N-terminal C2H2-type zinc finger (ZF) domain [22]. The HD domain consists of 60 amino acids forming three α-helices and is well-known as a DNA binding domain that is involved in plant growth, development and target gene expressions [19, 21, 23]. The ZF domain extensively controls regulatory proteins that promote DNA, RNA, protein, and DNA-RNA interactions, playing an important role in transcriptional and translational regulation [24]. ZF-HD proteins specifically bind to DNA sequences with an ATTA core consensus and form homodimers or heterodimers and play key roles in various biological processes, including stress responses and developmental regulation [21, 25]. The ZF-HD gene family is divided into two subfamilies (ZHD and MIF (mini zinc finger)) based on the conserved motifs and phylogenetic relationships. ZHD contains both (HD and ZF) domains, while MIF only has the ZF domain [20, 26].

Several studies have functionally characterized and highlighted the crucial role of ZF-HD genes including plant growth and development [27,28,29,30], hypocotyl elongation [30] and floral development [25]. Transient heterologous expression of soybean ZF-HD1 and ZF-HD2 genes in Arabidopsis protoplasts highlighted their response to pathogen infection and activation of CaM4 gene expression [31]. Recently, the cassava (Manihot esculenta Crantz) MeZHD7 gene was functionally characterized and found to be potentially involved in conferring resistance to cassava bacterial wilt (CBB) [32]. In addition, ZF-HD genes played significant roles in plant tolerance to various stresses such as drought tolerance in Arabidopsis thaliana [18], Chrysanthemum nankingense [33] and CsHDZ3 in tea (Camellia sinensis) [34]. Further, ZF-HD genes regulated the expression of OsDREB1 under cold stress in rice [35]. In addition, Yong et al., [36] reported that the LlZFHD4 gene was upregulated under cold, salt, and drought stress conditions in lily (Lilium lancifolium). Khatun et al., [37] also identified that ZF-HD genes in tomato (Solanum lycopersicum) were involved in regulating fruit development and enhancing stress tolerance under heat, cold, drought, salt, and abscisic acid (ABA) conditions. Moreover, HvZFHD1 enhanced the drought, salt, heat methyl jasmonate (MJ), salicylic acid (SA), and ethephone stress tolerance in barley (Hordeum vulgare) [38]. Numerous studies on several plants have shown that ZF-HD genes are involved in regulating stress responses at the transcriptional level, and similar functions are expected to be exerted in soybeans.

Recently, genome-wide analyses have identified and characterized the ZF-HD gene family members in several plant species including Arabidopsis (Arabidopsis thaliana) [26], tomato (Solanum lycopersicum) [37], wheat (Triticum aestivum) [39], chilli (Capsicum annuum) [40], barley (Hordeum vulgare) [20], apple (Malus domestica) [41], tea (Camellia sinensis) [42], cotton (Gossypium hirsutum) [43] and cassava (Manihot esculenta) [44], with distinct potential functions. Although ZF-HD genes have been studied in various species, so far, there is no report on the comprehensive identification and characterization of ZF-HD genes in the soybean genome. The whole genome sequencing of soybeans [45] provides strong support to study the comprehensive identification and characterization of ZF-HD genes and their possible functions.

In this study, ZF-HD genes were identified in the soybean genome through a screening process and further characterized through various in silico analyses including physicochemical properties, chromosomal localization, conserved domain, motifs, intron–exon structures, cis-elements, phylogenic, synteny, collinearity, transcription factors and miRNAs prediction, gene ontology and expression analyses. The relative expression profiles of selected candidate GmZF-HD genes in response to stress treatments were also determined by qRT-PCR in soybean leaves. This study provides the basis for further functional research on GmZF-HD genes and would be beneficial in improving the stress tolerance of soybeans by genetic manipulation.

Materials and methods

Identification of ZF-HD genes in the soybean genome

The genomic sequences of soybean (Glycine max) were downloaded (Wm82.a6.v1) from phytozome (https://phytozome-next.jgi.doe.gov). The ZF-HD hidden Markov models (HMMs) domain (PF04770) was downloaded from the Pfam database (http://pfam-legacy.xfam.org/) and used for HMMER via TBTool V 2.119 [46] under default mode. The conserved domain (ZF-HD_dimer) of ZF-HD genes was confirmed by SMART (http://smart.embl-heidelberg.de/) and NCBI Conserved Domain Database (https://www.ncbi.nlm.nih.gov/Structure/bwrpsb/bwrpsb.cgi). The molecular weight (M.W), isoelectric point (pl), instability index (II) and grand average of hydropathicity (GRAVY) were retrieved by the Expasy tool (https://web.expasy.org/protparam/). The prediction of protein subcellular localization was conducted by CELLO V 2.5 (http://cello.life.nctu.edu.tw/) and ZF-HD genes were renamed according to their chromosomal position in soybean.

Phylogenetic analyses

The evolutionary relationships between ZF-HD proteins were obtained by constructing a neighbor-joining phylogenetic tree between different plant species. The protein sequences of Arabidopsis thaliana (AtZF-HD), Oryza sativa (OsZF-HD), Solanum lycopersicum (SlZF-HD) and Malus domestica (MdZF-HD) were downloaded from Phytozome (https://phytozome-next.jgi.doe.gov), while Brassica rapa (BrZF-HD) and Pisum sativum (PsZF-HD) were downloaded from ensemble (https://plants.ensembl.org/index.html). All the ZF-HD protein sequences were aligned and subjected to a neighbor-joining tree with 1,000 bootstrap replicates using MEGA V10.1.8 software (https://www.megasoftware.net/) and visualized with iTOL (https://itol.embl.de/login.cgi). In addition, ZF-HD pairwise protein similarity/identity percentage SIAS (http://imed.med.ucm.es/Tools/sias.html).

Gene structure, conserved motif and cis-regulatory elements analysis

The GmZF-HD gene structure information was retrieved from soybean genomic files and conserved motifs were identified by the online MEME suite (https://meme-suite.org/meme/tools/meme) with a maximum of ten motifs. Moreover, plant growth and development, phytohormones responses, light and stress-responsive related cis-regulatory elements in GmZF-HD upstream 2000 bp regions were predicted by PlantCARE database (https://bioinformatics.psb.ugent.be/webtools/plantcare/html/). The following results were visualized by TBtools.

Synteny and collinearity analysis

The synteny and collinearity analyses of GmZF-HD genes were performed as described in our recent work [47]. Synteny analyses were performed between GmZF-HD and other species of ZF-HD genes including dicotyledonous (A. thaliana; AtZF-HD), (S. lycopersicum; SlZF-HD), (M. domestica; MdZF-HD), (B. rapa; BrZF-HD), (P. sativum; PsZF-HD) and monocotyledonous (O. sativa; OsZF-HD) species respectively. The Ka (non-synonymous) and Ks (synonymous) values for the duplicated GmZF-HD genes were calculated according to Rizwan et al., [47]. Moreover, multiple collinearity analyses of GmZF-HD genes were performed with other species including G. max, A. thaliana, B. rapa, P. vulgaris, M. domestica, V. unguiculata, P. sativum, C. cajan S. tuberosum, O. sativa and S. lycopersicum according to Rizwan et al., [47].

Protein–protein interaction network and 3D structure analysis

The online STRING database (https://string-db.org/cg) was used to perform GmZF-HD protein–protein interactions with a maximum of 10 interactions and protein 3D models were predicted by Phyre2 (http://www.sbg.bio.ic.ac.uk/phyre2/html/) with default parameters.

Prediction of miRNAs, TFs targeting GmZF-HD and GO annotation analysis

Prediction of putative miRNAs, TFs and gene ontology (GO) annotation analyses were performed as described in our recent work [47].

Expression profiling of GmZF-HD genes in various tissues

The gene expression data (FPKM; Fragments per Kilobase of transcript per Million mapped reads) from various soybean tissues such as seed, root, leave and shoot was downloaded from the Phytozome (https://phytozome-next.jgi.doe.gov/info/) and expressed in heat maps by TBtools after transferred to Log2FC.

Plant growth and abiotic-stress treatments of soybean seedlings

The seeds of black soybean were purchased from Shaouguang Xinxiran Horticulture Co. Ltd. and were surface sterilized using 2% NaClO for 10 min and washed subsequently five times with distilled water. The seeds were sown on planting trays in a growth chamber (YK-ZW-300D3, Hefei, China) containing half Hoagland nutrient solution (HNS, pH 6.0) under a temperature of 25 ± 2 °C, 70–75% relative humidity, and 16/8 light/dark conditions. The seedlings with consistent growth were subjected to various abiotic stress treatments after the first trifoliate stage (V1 stage) [48, 49]. Cold and heating stresses were achieved by transferring plants to an identical chamber (growth chamber; YK-ZW-300D3, Hefei, China) set to 4°C and 37°C, respectively [49]. The half HNS was supplemented with 150 mM NaCl for salt stress [50] and 20% PEG6000 (w/v) for drought/osmotic stress [51]; 0.1 mM of abscisic acid (ABA) and salicylic acid (SA) for hormone stress [51, 52]; 0.1 mM CdCl2, 0.1 mM CoCl2, 0.5 mM FeSO4, 1.0 mM MnSO4, and 0.5 mM ZnSO4 for metal ion stress [53]. Seedlings under half HNS were used as a control group. Soybean leaves from twenty individual plants (three biological repeats) of various treatments were collected after 24 h of treatments and immediately frozen in liquid nitrogen and stored at -80°C for subsequent analysis.

RNA extraction and qRT-PCR analysis

The total RNA extraction, cDNA-synthesis, qRT-PCR and relative expression analyses were performed according to the protocol described in our recent work [47]. Three replicates were used for each reaction and the relative gene expression levels were calculated using the 2–ΔΔCT method with GmACTIN as an internal reference gene [54]. All the primers used for qRT-PCR are listed in Table S1. The graphs were generated using GraphPad Prism V 9.0 (https://graphpad.com).

Data analysis

The mean ± SD values of three replicates were used and the least significant difference (LSD) among the mean values was determined by Tukey’s test (p < 0.05 using) in Statistix V 8.1 software (https://statistix.informer.com/8.1/). Heatmaps depicting the gene expression data were visualized through TBTool V 2.119 (https://github.com/CJ-Chen/TBtools-II/releases). Bar plots were calculated and visualized through GraphPad Prism V 9.0 (https://graphpad.com).

Results

Identification, physicochemical features and chromosomal location of GmZF-HD genes

After removing repetitive and non-redundant sequences, a total of 51 ZF-HD genes were obtained in the soybean genome and were renamed according to their order on chromosomes (GmZF-HD1 to GmZF-HD51) (Table 1). All the identified genes were unevenly distributed on 17 soybean chromosomes except for chromosomes Gm03, Gm10 and Gm15. Among chromosomes, chromosome Gm09 had the maximum number of genes (6 genes) followed by chromosomes Gm02, Gm08 and Gm20 (5 genes) respectively. The least only one gene was located on each chromosome including Gm11, Gm12 and Gm14 respectively (Table 1). Physicochemical analyses exhibited a wide range of protein lengths ranging from 73 aa (GmZF-HD31) to 358 aa (GmZF-HD19/47), with molecular weight (M.W) ranging from 8.11 to 39.16 KDa respectively. Similarly, the protein isoelectric points (pI) varied from 6.5 (GmZF-HD9/10/15) to 9.87 (GmZF-HD22) and the mean of hydropathicity (GRAVY) scores ranged from -1.20 (GmZF-HD2) to -0.21 (GmZF-HD43). The majority of 46 GmZF-HD genes consisted of 1 exon and only 5 GmZF-HD genes consisted of 2 exons including GmZF-HD8/9/10/15 and GmZF-HD38 (Table 1). In addition, the predicted subcellular localization results showed the possible localization of GmZF-HD genes is in periplasmic, extracellular, and in periplasmic (Table 1). The detailed results of GmZF-HD genes have been provided in Table 1 and the protein sequences are available in Table S2.

Table 1 Physicochemical features of GmZF-HD genes
Full size table

Phylogenetic analyses of GmZF-HD proteins

The phylogenetic tree categorized GmZF-HD proteins into eight major groups (I-VIII), similar to the tomato [37] (Fig. 1). The ZF-HD protein sequences used for the phylogenetic tree have been provided in Table S2. Group-I contained the maximum of 36 ZF-HD members followed by VIII, VII, VI, V, IV, III and II groups contained 36, 30, 20, 16, 14, 11, and 11 ZF-HD members (Fig. 1). Group VIII contained a maximum of 19 GmZF-HD (37% of the total), while groups IV and V contained at least three GmZF-HD members. Overall, GmZF-HD proteins clustered with other species of ZF-HD proteins that show a close relationship with their ancestral species and may perform similar functions. Conversely, the cluster of group II only consists of 3 AtZF-HD and 8 BrZF-HD proteins, indicating that during the evolution of soybean species, some genes may have been lost. In addition, the results of GmZF-HD proteins pair-wise similarity percentage exhibited that GmZF-HD7/48 and GmZF-HD49 proteins had the highest similarity (93–96%), followed by GmZF-HD12/13 and GmZF-HD17 (77%) respectively. Similarity percentage results are consistent with the phylogenetic analysis, in which the genes with the highest similarity percentage clustered in the same groups such as group I (GmZF-HD12/13 and GmZF-HD17) and group VIII (GmZF-HD7/48 and GmZF-HD49) (Fig. 1, Table S3).

Fig. 1
figure 1

Neighbor-joining (NJ) phylogenetic tree based on the amino acid sequences alignment among Arabidopsis thaliana (AtZF-HD), Glycine max (GmZF-HD), Oryza sativa (OsZF-HD), and Solanum lycopersicum (SlZF-HD), Malus domestica (MdZF-HD), Brassica rapa (BrZF-HD) and Pisum sativum (PsZF-HD) sequences with 1000 bootstraps. All the ZF-HD members were divided into 8 groups and presented in different colors

Full size image

Gene structure, domain and conserved motifs analysis of GmZF-HD genes

Furthermore, the gene structure, conserved domain and conserved motifs of GmZF-HD genes were analyzed (Fig. 2). All the 51 GmZF-HD genes were phylogenetically (NJ) clustered into six groups (I-VI, Fig. 2A) and group I contained the maximum 16 members followed by groups VI (11 members), V (8 members), II-III (6 members) and least 3 members in group IV (Fig. 2A). Further, the MEME online tool was utilized and a total of 10 conserved motifs (motif 1–10) were identified among all GmZF-HD proteins (Fig. 2B). The results exhibited that motifs 1, 2 and 6 were the most common motifs in all proteins and members of group I contained all of the 10 motifs respectively. In terms of genes, GmZF-HD7/49/48/50/41/24/23/25/1/4/8/36/15/16/11/9/10/19 and GmZF-HD44 genes contained the maximum of 10 motifs followed by GmZF-HD3/32/27 and GmZF-HD28 genes (9 motifs), GmZF-HD5/37/38/51/20/39/34/46/2/6 and GmZF-HD47 genes (8 motifs), GmZF-HD18 and GmZF-HD26 genes (7 motifs) GmZF-HD22/33/21 and GmZF-HD14 genes (6 motifs), GmZF-HD30 gene (5 motifs) and GmZF-HD42/29/43/31/35/45/13/40/12 and GmZF-HD17 genes contained the minimum of 3 motifs respectively (Fig. 2B). The conserved motif results are consistent with the GmZF-HD genes phylogenetic tree and clustered in the same group (Fig. 2A-B). GmZF-HD genes motifs logo has been provided in Fig-S1.

Fig. 2
figure 2

The gene structure, domain, and motif analysis of GmZF-HD genes. A Unrooted neighbor-joining phylogenetic tree among GmZF-HD proteins and clustered into VI clades, B Ten conserved motifs were represented via boxes, and different colors represent different motifs, C Conserved domain in GmZF-HD, D GmZF-HD gene structures; yellow color represent the exons (CDS), green color represents the untranslated 5′ and 3′-regions, and gray color indicates the introns

Full size image

Additionally, the CDD database analysis revealed that all of the 51 GmZF-HD proteins contained a highly conserved “ZF-HD_dimer” (PF04770) domain, confirming their classification within the ZF-HD gene family (Fig. 2C). Exon–intron analysis showed that most GmZF-HD genes (43 out of 51) consisted of one CDS and lacked introns. For example, GmZF-HD8/9/10 and GmZF-HD38 genes contained only two CDS, while GmZF-HD27/29 and GmZF-HD42 genes had one CDS with no introns (Fig. 2D). Overall, the consistency between phylogenetic grouping, conserved motifs, domain structures, and gene architecture suggests that GmZF-HD genes share highly conserved amino acid sequences, and genes within the same group may have similar functions.

Cis-regulatory elements in the upstream sequence of GmZF-HD genes

The cis-regulatory elements (CREs) analysis was conducted to understand the possible involvement of GmZF-HD genes in plant growth, development and stress responses. A total of 1750 CREs were identified and were categorized into four main groups (Fig. 3A). The results revealed that the highest CREs were from light-responsive (724, 41%) followed by phytohormone-responsive (547, 31%), stress-responsive (328, 18%) and the least from metabolism (173, 10%) (Fig B). Overall, all the 51 GmZF-HD genes contained light-responsive CREs followed by phytohormone-responsive (45 genes), stress-responsive (44 genes) and only 34 genes contained growth and development-related CREs (Fig. 3 Ba). The main categories of CREs were further described based on their distributions and subfamilies such as light responsive (724 CREs) contained the maximum of Box-2 (278, 39%), GT1-motif (89, 12%), TCT-motif (57, 10%) CERs and the least AT1-motif (8, 1%) (Fig. 3 Bd).

Fig. 3 
figure 3

The cis-regulatory element analysis in GmZF-HD genes (A-B), sum number of GmZF-HD genes involved in four categories of cis-elements (Ba), percentage (%) ratio of the numerous cis-elements from each category is presented in pie charts including plant growth and development (Bb), phytohormones-responsive (Bc), light-responsive (Bd) and stress-responsive (Be)

Full size image

In the phytohormone responsive CREs group (547, 31%), the maximum cis-elements belongs to ABRE (138, 25%), CGTCA/TGACG-motifs (111, 20%), ERE (78, 14%) and the least to TGA-box (1, 0.1%) (Fig. 3 Bc) CREs. In the stress-responsive category (328, 18%), the maximum CREs included ARE (328, 35%), MBS (48, 15%) and the least CARE (7, 2%) (Fig. 3 Be). In the metabolism category (173, 10%), the maximum CREs consisted of CAT-box (59, 34%), O2-site (46, 27%) and the least MSA-like (2, 1%) (Fig. 3 Ba). In addition, GmZF-HD2 and GmZF-HD5 contained the highest 49 and 44 CREs, whereas GmZF-HD30 was found to have the lowest 19 CREs in their promoter regions (Fig. 3 and Table S4). The presence of these cis-elements in the promoter regions of GmZF-HD genes indicated that GmZF-HD genes may play important roles in soybean growth and development and respond to different stresses while requiring further research.

Synteny analyses of GmZF-HD genes

The expansion of GmZF-HD genes was further investigated through gene duplication events. A total of 74 GmZF-HD duplicated genes were identified and all gene pairs underwent segmental duplication events, except for one tandem duplication (GmZF-HD4/GmZF-HD8) that was located on chromosome Gm02. The segmentally duplicated genes were located on different chromosomes including chromosomes Gm01, Gm04, Gm05, Gm06, Gm07, Gm08, Gm09, Gm11, Gm13, Gm17, Gm18 and Gm19. The largest number of duplicated genes were positioned on chromosomes Gm02 and Gm15 (Fig. 4A, Table S5). In addition, the Ka / Ks ratios of all duplicated GmZF-HD genes were found to be less than 1 (0.45) indicating that the GmZF-HD gene pairs had undergone purifying selection during evolution. The divergence rate between segmented and tandem GmZF-HD genes varied from 1.43 MYA (million years ago) (GmZF-HD13/GmZF-HD40) to 84.11 MYA (GmZF-HD6/GmZF-HD46) (Table 2). These results suggest that gene duplication played a crucial role in the development of the GmZF-HD gene family.

Fig. 4
figure 4

Circos illustrations of GmZF-HD gene duplication and multicollinearity analysis. A Gene duplication of ZF-HD genes in G. max (GmZF-HD). B Orthologous of GmZF-HD genes with AtZF-HD genes (A. thaliana) and MdZF-HD genes (M. domestica) and BrZF-HD genes (B. rapa). C Orthologous of GmZF-HD genes with OsZF-HD genes (O. sativa), SlZF-HD genes (S. lycopersicum) and PsZF-HD genes (P. sativum). The background gray lines show all the syntenic blocks among genomes and different colored connected lines show duplication link regions among the ZF-HD genes. The first ring from outside indicates the approximate location of ZF-HD genes and are labeled with a short gray line outside with gene names. The second and third rings from the outside represent the density of genes on the chromosomes. The blue-to-red scale bar indicates the number of SNPs within 1 Mb window size. D-E Multicollinearity analysis of GmZF-HD genes with the genomes of the different plant species including A. thaliana, B. rapa, P. vulgaris, M. domestica, P. sativum, C. cajan, S. tuberosum, O. sativa, S. lycopersicum and V. unguiculata species. The red lines represent the GmZF-HD syntenic genes between other species genomes and the gray lines in the background represent all of the orthologous genes from the above species genomes. The detailed results information can be found in Table S5-6

Full size image
Table 2 Ka/Ks analysis and estimated divergence time for duplicated GmZF-HD genes
Full size table

Additionally, a comprehensive synteny analysis was performed among GmZF-HD genes and ZF-HD genes from other plant species including dicotyledonous AtZF-HD, SlZF-HD, BrZF-HD, PsZF-HD and MdZF-HD and monocotyledonous OsZF-HD. Among all 51 GmZF-HD genes, only 33 GmZF-HD genes showed synteny with ZF-HD genes from other plants and a total of 314 duplicated gene pairs were found between them. In terms of species, GmZF-HD genes showed the highest synteny with B. rapa (BrZF-HD, 108 gene pairs) followed by M. domestica (MdZF-HD, 75 gene pairs), A. thaliana (AtZF-HD, 59 gene pairs), S. lycopersicum (SlZF-HD, 30 gene pairs), P. sativum (PsZF-HD, 25 gene pairs) and the lowest with O. sativa (OsZF-HD, 17 gene pairs). Among GmZF-HD genes, the GmZF-HF28 gene exhibited the highest homology of 15 gene pairs with ZF-HD genes from other species, such as maximum of 6 gene pairs with B. rapa (BrZF-HD3/5/7/16/17/19), 4 gene pairs with A. thaliana, (AtZF-HD8/9/10/11), 4 gene pairs with M. domestica (MdZF-HD7/14/15) and minimum 1 gene pair with S. lycopersicum (SlZF-HD18). While the GmZF-HD31 gene was found to be paired with only one SlZF-HD9 gene (Fig. 4B-C, Table S5).

Moreover, to further understand the orthologous homology of GmZF-HD genes between the genomes of different plant species, we further performed multicollinearity analysis on GmZF-HD genes with other plant species genomes including A. thaliana, B. rapa, P. vulgaris, M. domestica, P. sativum, C. cajan, S. tuberosum, O. sativa, S. lycopersicum and V. unguiculata (Fig. 4D-E, Table S6). The results exhibited that among 51 GmZF-HD genes, only 37 GmZF-HD genes exhibited collinearity with 615 genes from all the above species. Further, the highest degree of collinearity was found between G. max (GmZF-HD) and M. domestica (108 genes) followed by V. unguiculata (103 genes), P. vulgaris (86 genes), B. rapa (74 genes), A. thaliana (64 genes), S. lycopersicum (62 genes), S. tuberosum (48 genes), C. cajan (28 genes), P. sativum (25 genes) and the lowest was in O. sativa (17 genes) genomes (Fig. 4D-E, Table S6). Synteny and orthologous collinearity results exhibited that GmZF-HD genes were highly homologous to dicot species compared with monocots, suggesting that they may belong to a common ancestor and may perform similar functions.

Protein-to-protein network interaction and 3D models of GmZF-HD proteins

The results of the STRING database showed that all of the 51 GmZF-HD proteins showed homology with Arabidopsis ZHD and MIF proteins and interacted with the identified proteins (Fig. 5A, Table S7). It was found that the maximum of 19 proteins including GmZF-HD1/7/8/9/10/11/14/15/16/21/23/24/25/33/36/41/48/49 and GmZF-HD50 were homologous with AtZHD1-2 protein, followed by 7 proteins (GmZHD13/29/31/35/42/43 and GmZHD45) with AtMIF2, 6 proteins (GmZf-HD3/5/22/32/37 and GmZF-HD38) with AtZHD6, 5 proteins (GmZF-HD18/20/26/39 and GmZF-HD51) with AtZHD4 protein, 5 proteins (GmZF-HD2/6/19/34 and GmZF-HD47) with AtZHD11 protein, 3 proteins (GmZF-HD27/28 and GmZF-HD44) with AtZHD10 and 3 proteins (GmZFHD12/17 and GmZF-HD40) showed homology with AtMIF3. Moreover, only one GmZF-HD46 was homologous with AtZHD8, GmZF-HD4 with AtZHD1 and GmZF-HD30 with AtMIF1 respectively. Furthermore, the homologous proteins (AtMF1/2/3, AtZHD1/4/6/8/10/100 and AtZHD1-2) were further interacted with different proteins including ATHB-21, SPEAR2/4, ZHD5/6/7/9, LCR81, NET4B, GASA7, TRFL8, DOF3.4, PRH, IBL1, KNU, TBL31 and MJB20.15, respectively (Fig. 5A, Table S7). The thicker the line between proteins reveals a strong interaction and may perform have functions.

Fig. 5
figure 5

Protein–protein interaction and predicted 3D models of GmZF-HD proteins. A The network was constructed using online STRING software. The proteins are displayed at network nodes with 3D structure of the proteins within nodes and the colors of the lines indicate different data sources. The higher the interaction coefficient, the thicker the network lines between proteins and vice versa. B Predicted 3D models of GmZF-HD proteins by homology modeling using the online Phyre2 server with default mode

Full size image

In addition, we also performed 3D modeling of GmZF-HD proteins using Phyre2 under default parameters. Overall, all 51 GmZF-HD proteins were modeled with 7 unique templates including d1wh7a, d1wh5a, c3ufkA, c1p58C, c6ra0A, c7nf9A and c1s1iY. The results exhibited that the confidence levels of GmZF-HD proteins with the templates were approximately 98%, sequence identities of 27–71% and alignment coverage levels of 21–52% (Fig. 5B, Table S8). A maximum of 38 proteins (GmZF-HD1/2/3/5/10/11/14/15/16/18/19/20/21/23/24/25/26/27/28/31/ 32/33/34/36/37/38/39/41/44/46/47/48/49/50 and GmZF-HD51) were modeled with “d1wh7a” template, followed by 7 proteins (GmZF-HD3/12/13/17/29/30/ and GmZF-HD45) with “c3ufkA” template and 2 proteins (GmZF-HD4 and GmZF-HD22) were modeled with “d1wh5a” template. However, the following proteins were modeled individually including GmZF-HD31 modeled with “c1p58C”, GmZF-HD35 with “c6ra0A” GmZF-HD40 with “c7nf9A”, and GmZF-HD4 was modeled with “c1s1iY” template, respectively (Fig. 5B, Table S8).

Prediction of miRNAs and TFs in GmZF-HD genes

Recently, miRNAs have been demonstrated to play crucial roles in various plant signal transduction and stress responses. Thus, in this study, 94 gma-miRNAs belonging to 41 different miRNA families were identified that targeted 40 GmZF-HD genes, except for the following 11 genes including GmZF-HD13/14/24/25/29/31/34/40/42/46 and GmZF-HD48 genes (Fig. 6, Table S9). Among the targeted GmZF-HD genes, GmZF-HD6 was found to be the highest targeted by 7 gma-miRNAs belonging to 6 different miRNA families (gma-miR169j-3p, gma-miR9742, gma-miR4993, gma-miR1535, gma-miR4397-5p and gma-miR9752), followed by GmZF-HD2 targeted by 6 miRNAs (gma-miR9742, gma-miR5774, gma-miR169j-3p, gma-miR4993 and gma-miR4397-5p) and GmZF-HD5 targeted by 5 miRNAs (gma-miR5041-3p, gma-miR1535, gma-miR9752, gma-miR2108b and gma-miR5376) respectively. Interestingly, only one miRNA from different families has targeted the following genes including GmZF-HD4/7/8/9/10/11/12/17/ 21/23/30/33/35/39/41/43/44/49 and GmZG-HD51. Moreover, gma-miR4993 was found to be the unique miRNA that targeted 10 genes including GmZF-HD1/2/6/8/19/27/28/36/41 and GmZF-HD47 (Fig. 6A-B, Table S9). Consequently, further studies are needed to confirm the biological role of predicted miRNAs and targeted GmZF-HD genes.

Fig. 6
figure 6

The putative miRNA targeted GmZF-HD genes. A Network illustration of predicted miRNA targeting GmZF-HD genes. Green color boxes represent GmZF-HD genes and lavender color boxes represent miRNAs. B The schematic diagram indicates the GmZF-HD genes targeted by miRNAs

Full size image

Furthermore, plant TFs have been reported to be involved in numerous functions, therefore, using the plantregmap database, a total of 10,444 TFs were identified in upstream 1000 bp regions of all 51 GmZF-HD genes that belong to 25 distinct TF families (Fig. 7, Table S10). The results revealed that the TF families with the highest abundance were ERF (6858), LBD (1287), BBR-BPC (532), MYB (457), C2H2 (225), GATA (191), TALE (179), MIKC_MADS (158), AP2 (130), TCP (105), bHLH (65), Dof (52), Nin-like (44), NAC (37), bZIP (34), Trihelix (33) and CAMTA (17). Whereas the least abundant TF families were SBP and C3H with only one TF. In addition, all GmZF-HD genes were targeted by different TF families, for example, the top 10 highly targeted genes were GmZHD8 (882), GmZHD36 (711), GmZHD44 (562), GmZHD27 (557), GmZHD28 (557), GmZHD15 (511), GmZHD16 (511), GmZHD1 (462), GmZHD4 (451) and GmZHD34 (444) (Fig. 8B, Table S10). Whereas, the least targeted genes were GmZHD12 (4), GmZHD33 (4) and GmZHD22 (2) respectively (Fig. 8, Table S10). Additionally, the existence of different types of TFs suggests that they can be considered for further functional studies.

Fig. 7
figure 7

Putative transcription factor in GmZF-HD genes and regulatory network analysis. A Green to dark red color shows (low to high) enrichment, lines represent the co-expression network between genes and TFs and node size represents the degree of interaction and enrichment between genes and TFs. B Top 10 highly enriched and targeted GmZF-HD genes with TFs, green nodes represent the TFs and red colors show the GmZF-HD genes

Full size image
Fig. 8
figure 8

GO annotation analyses of GmZF-HD genes. A Highly enriched GO terms in GmZF-HD genes including molecular (blue color), biological (red–orange color) and cellular (yellow-orange). Further detailed annotation results and numerous significantly enriched terms can be found in Table S11

Full size image

GO annotation analysis of GmZF-HD genes

To further understand the dynamic role of GmZF-HD genes in molecular function (MF), cellular component (CC) and biological process (BP), gene ontology (GO) annotation analysis was performed and identified 66 highly enriched GO terms. Among the identified GO terms, BP was found to have a maximum of 45 terms, followed by MF with 14 terms and CC with at least 7 terms (Fig. 8, Table S11). The highly enriched BP terms were related to floral meristem determinacy (GO:0010582), meristem determinacy (GO:0010022), S-glycoside metabolic process (GO:0016143), glycosinolate metabolic process (GO:0019757, GO:0019760), plant-type ovary development (GO:0035670), post-embryonic development (GO:0009791), carpel development (GO:0048440), reproductive system development (GO:0061458) and reproductive structure development (GO:0048608). The highly enriched MF terms were related to protein dimerization activity (GO:0046983), protein binding (GO:0005515), heterocyclic compound binding (GO:1,901,363), organic cyclic compound binding (GO:0097159) and DNA binding (GO:0003677). Whereas the highly enriched CC terms were related to the nucleus (GO:0005634), intracellular membrane-bounded organelle (GO:0043231), membrane-bounded organelle (GO:0043227), intracellular organelle (GO:0043229), organelle (GO:0043226) and obsolete intracellular part (GO:0044424) (Fig. 8, Table S11). GO annotation analysis indicated that GmZF-HD genes may play functional roles in different biological, cellular, and molecular processes that require further studies.

Expression profiles of GmZF-HD genes in different soybean tissues

The expression patterns of all 51 GmZF-HD genes were analyzed across different soybean organs/tissues including seeds, roots, shoots and leaves under various conditions, based on public data collected from the Phytozome database. The FPKM expression values were converted to Log2FC and were visualized in heatmaps (Fig. 9, Table S12). The results showed that among all 51 GmZF-HD genes, a maximum of 18 GmZF-HD genes (35%) were positively or negatively expressed in seeds, followed by 16 genes (31%) in the shoot, 14 genes (27%) in roots and minimum 8 genes (15%) in leaves under different conditions. In addition, the 33 genes including GmZF-HD1/3/8/11/16/17/18/19/20/21/23/24/27/28/29/30/31/32/33/35/36/37/38/41/42/43/44/46/47/48/49 and GmZF-HD50 were not expressed in all soybean tissues, representing that they may not play roles in these tissue under the specific conditions.

Fig. 9
figure 9

FPKM-based expression profiles of GmZF-HD genes in different tissues including seeds, roots, leaves and shoots under different conditions. FPKM values were transformed by log2FC, and heatmaps were constructed using TBtools software. The side scale bar represents the expression profile from lowest (blue) to highest (red)

Full size image

Conversely, among the expressed GmZF-HD genes, the expression levels also varied in different tissues, such as, the expression levels varied from -8.42 to 4.38 fold in seeds, -6.65 to 4.20 fold in leaves, -5.48 to 3.98 fold in shoots and -2.68 to 10.11 fold in roots. Furthermore, several genes were significantly highly expressed in different tissues under different conditions such as GmZF-HD15 (10.11 fold) and GmZF-HD40 (9.24 fold) in roots under salt stress, GmZF-HD2, GmZF-HD5, GmZF-HD6, GmZF-HD13, GmZF-HD22, GmZF-HD45 (> 6.0 fold) and GmZF-HD51 (> 4.0 fold) in shoots and leaves under drought stress conditions (Fig. 9, Table S12). Overall, it is suggested that genes (GmZF-HD2/5/6/13/15/22/40 and GmZF-HD45) with highly significant expressions can be considered for further genetic engineering to further evaluate their involvement in soybean growth, development and response to stresses.

Expression profiling of GmZF-HD genes in response to different stress treatments

To further investigate the putative functions of GmZF-HD genes in response to different abiotic stresses, we performed qRT-PCR analysis of ten candidate GmZF-HD genes in soybean leaves subjected to cold, heat, salt, drought, phytohormones and heavy metals stress treatments. The candidate GmZF-HD genes were carefully selected based on their significant involvement in cis-elements and FPKM expressions under different conditions. The qRT-PCR expression profiles revealed a different response of individual GmZF-HD genes under specific stress and were mainly consistent with the FPKM results, thereby validating the reliability of existing results (Figs. 10, 11 and 12). Overall, the majority of GmZF-HD genes showed upregulated expression levels under different stress treatments (Figs. 10, 11 and 12). Interestingly, GmZF-HD5 (Figs. 10b, 11b and 12b), GmZF-HD6 (Figs. 10c and 11c), GmZF-HD13 (Fig. 12d), GmZF-HD39 (Figs. 10g and 12g) and GmZF-HD45 (Fig. 12i) genes were significantly upregulated (2.5 to 8.0 fold) under stress treatments, indicating they may play crucial role in enhancing soybean plant tolerance to stress conditions.

Fig. 10
figure 10

qRT-PCR-based relative expressions of GmZF-HD genes in soybean leave tissues under control, cold, heat, salt and drought treatments. The relative gene expression levels were calculated using the 2–ΔΔct. Vertical bars represent means ± SD (n = 3). The asterisks *, **, and *** show significance levels at p ≤ 0.05, p ≤ 0.01, and p ≤ 0.001 and ns means non-significant among treatment and control samples according to the Students t-test

Full size image
Fig. 11
figure 11

qRT-PCR-based relative expressions of GmZF-HD genes in soybean leave tissues under control, ABA and SA treatments. The relative gene expression levels were calculated using the 2–ΔΔct. Vertical bars represent means ± SD (n = 3). The asterisks *, **, and *** show significance levels at p ≤ 0.05, p ≤ 0.01, and p ≤ 0.001 and ns means non-significant among treatment and control samples according to the Students t-test

Full size image
Fig. 12
figure 12

qRT-PCR-based relative expressions of GmZF-HD genes in soybean leave tissues under control, CdCl2, CoCl2, FeSO4, MnSO4, and ZnSO4 treatments. The relative gene expression levels were calculated using the 2–ΔΔct. Vertical bars represent means ± SD (n = 3). The asterisks *, **, and *** show significance levels at p ≤ 0.05, p ≤ 0.01, and p ≤ 0.001 and ns means non-significant among treatment and control samples according to the Students t-test

Full size image

In addition, under cold stress, the expression levels of GmZF-HD5/6/34/ and GmZF-HD39 genes were upregulated (1.2 to 2.4 fold), while GmZF-HD2/13/15/40/45 and GmZF-HD51 genes were downregulated compared to control (Fig. 10). Under heat stress, the expression levels of GmZF-HD2/5/6/13/15/34/39 and GmZF-HD40 were upregulated (1.2 to 2.7 fold), whereas GmZF-HD45/51 expression levels were downregulated compared to control (Fig. 10). Under salt stress, GmZF-HD2/6 and GmZF-HD39 genes were upregulated (1.2 to 2.1 fold), however, GmZF-HD15/34/45 and GmZF-HD51 were downregulated. Remarkably, GmZF-HD5/13 and GmZF-HD40 genes showed no change in expressions compared with the control. Under PEG-induced drought stress, the expression level of GmZF-HD2/515/34 and GmZF-HD51 were upregulated (1.2 to 3.9 fold), while GmZF-HD6/13/40/45 and GmZF-HD51 were downregulated compared to control (Fig. 10).

Under ABA treatment, Gm2/5/34 and GmZF-HD39 genes were upregulated (1.2 to 1.6 fold), while GmZF-HD6/13/15/40/45 and GmZF-HD51 were downregulated compared to control (Fig. 11). Under SA treatment, the expression levels of all the tested genes were upregulated (1.2 to 5.5 fold) compared to the control (Fig. 11). Furthermore, the expression levels of GmZF-HD genes also varied, when treated with different metal ions. For example, under CdCl2 treatment, the expression level of GmZF-HD5/40 and GmZF-HD45 were upregulated (1.2 to 6.9 fold), while GmZF-HD2/6/13/15/34/39 and GmZF-HD51 were downregulated compared to control (Fig. 12). Under CoCl2 treatment, GmZF-HD2/5/6 and GmZF-HD39 genes were upregulated (1.2 to 1.5 fold), while GmZF-HD13/15/34/40/45 and GmZF-HD51 were downregulated compared to control (Fig. 12). Under MnSO4 treatment, GmZF-HD2/5/13/15/34/39 and GmZF-HD45 showed upregulation (1.2 to 3.5 fold), while GmZF-HD6/40 and GmZF-HD51 showed downregulation compared to control (Fig. 12). Under ZnSO4 treatment, the expression levels of all tested genes were upregulated except for GmZF-HD40 (Fig. 12). Interestingly, all tested GmZF-HD genes were upregulated under FeSO4 treatment than control (Fig. 12). Taken together, all the GmZF-HD genes were significantly induced under cold, heat, drought, ABA, SA, CdCl2, CoCl2, MnSO4, ZnSO4 and FeSO4 stress treatments, highlighting their potential roles in response to stresses in soybean. Genes that were significantly up-regulated under stress treatments compared with the control, such as GmZF-HD5/6/13/39 and GmZF-HD45, could be considered for functional studies in soybean improvement under stress conditions.

Discussion

Soybean is an important economic crop source for edible oil, protein, and food nutrition, however, its production faces different environmental challenges from biotic and abiotic stresses [55,56,57,58,59,60,61,62]. Zinc finger homeodomain (ZF-HDs) genes play crucial roles in plant growth, development and various stress responses [18,19,20]. ZF-HD genes have been identified in different plant species including Arabidopsis [26], tomato [37], wheat [39], chilli [40], barley [20], apple [41], tea [42] and cotton [43]. Despite its significant importance, a comprehensive analysis of ZF-HD genes in the soybean genome remains unexplored.

In this, 51 GmZF-HD genes were identified in the soybean genome that was unevenly distributed on 11 chromosomes (Table 1) and were inconsistent and somewhat higher than that of Arabidopsis (17) [26], 22 in tomato [37], 37 in wheat [39], 11 in chili [40], 8 in barley [20], 19 in apple [41], 37 in tea [42] and 37 in cotton [43]. Variations in the numbers of the same family members among different plant species may be due to genome size or WGD that cause an increase or decrease in their numbers [63]. In addition, the physiochemical features of GmZF-HD showed variations in molecular weight (M.W), isoelectric points (pI), GRAVY and exons (Table 1). GmZF-HD proteins exhibited a conserved “ZF-HD_dimer” domain and are consistent with those observed in other plant species [20, 37, 41, 43], indicating that GmZF-HD proteins are conservative as like other plant species.

In addition, a phylogenetic tree between soybean ZF-HD proteins and other plant ZF-HD proteins was constructed (Fig. 1), and divided into eight major groups (I-VIII), which aligns with previous reports [37, 43] (Fig. 1). It has been stated that genes that are highly conserved among species and exist in the same evolutionary branch of the phylogenetic tree belong to similar domains and may perform similar functions [64]. For example, it was found that GmZF-HD2/6/19/27/28/34/44/46/ and GmZF-HF47 clustered with AtZF-HD5 and AtZF-HD11 (Fig. 1, VII). AtZFHD5 and AtZFHD11 have been reported to be associated with drought stress tolerance [18, 65]. Moreover, GmZF-HD1/4/7/8/9/10/11/15/16/23/24/25/30/36/43/48/49 and GmZF-HD50 proteins were clustered into the same group with OsZF-HD1 and OsZF-HD2 proteins (Fig. 1, I). OsZF-HD1 and OsZF-H2 genes have been reported to play roles in plant morphogenesis and leaf development [27]. The MEME results exhibited that there were 10 conserved motifs and motifs 1, 2, and 6 were the most common motifs in all GmZF-HD proteins (Fig. 2). Our results are consistent with previous reports on ZF-HD motifs including apple and tomato [19, 41]. CDD results showed that all 51 GmZF-HD proteins contained a highly conserved “ZF-HD_dimer” (PF04770) domain (Fig. 2C), and gene structure analysis revealed that the majority of GmZF-HD genes (43 out of 51) were intronless (Fig. 2D). These results are consistent with previous reports on intronless ZF-HD gene structures [22, 26, 37], indicating that GmZF-HD proteins are structurally conservative.

Additionally, we performed cis-elements analysis in the promoter region of GmZF-HD genes and found different cis-elements related to phytohormone, development and stress responsiveness (Fig. 3, Table S4). Moreover, several cis-elements were dominant including Box 2, GT1-motif, TCT-motif, ABRE, CGTCA/TGACG-motifs, ERE and TGA-box, ARE, MBS and CARE; CREs, CAT-box and O2-site. Feng et al., [66] reported that ABRE cis-element responded to ABA and was involved in abiotic stress tolerance. CGTCA-motif exhibited stress tolerance [67]. Moreover, under abiotic stresses the upregulation of various rice genes was found due to the cis-elements in the promoter region [67]. In addition, gene duplication plays an important role in the expansion and divergence of gene families and it can be whole-genome duplication and tandem or segmental duplication [68].

In this study, 73 GmZF-HD gene pairs were developed by segmental duplication and only one gene pair (GmZF-HD4/GmZF-HD8) was developed by tandem duplication (Fig. 4A, Table S5), indicating that segmental duplication played a crucial role in the evolution of GmZF-HD gene family. The soybean GmZF-HD gene duplication pattern aligns with observations in other plants ZF-HD genes [19, 41, 69, 70]. The Ka/Ks ratios for duplicated GmZF-HD genes were less than 1 (0.45), indicating that they underwent strong purifying selections during the evaluation process and the rate of divergence was 1.43 to 84.11 MYA (Table 2). Similar trends of gene duplication and strong purification have been reported in the ZF-HD genes of cotton [43], Chinese cabbage [22], and tomato [37]. Comprehensive synteny and collinearity analysis of GmZF-HD genes revealed that soybean ZF-HD genes showed the highest homology with dicots including Chinese cabbage, apple, Arabidopsis, tomato, cowpea, pigeon pea and pea than that of monocots including rice (Fig. 4, Table S6), suggesting that they may have developed from same ancestors with similar functions. Synteny and collinearity results of GmZF-HD genes are consistent with previous soybean gene families, which also reported higher homology with the dicots species than monocots [71, 72].

Protein–protein network interactions play a key role in predicting the structures and functions of target proteins [73, 74]. The STRING database results exhibited that GmZF-HD proteins showed orthology with identified and characterized Arabidopsis ZF-HD and MFI proteins and interacted with different proteins including ATHB, SPEAR, ZHD, LCR81, NET4B, GASA, TRFL, DOF, PRH, IBL, KNU, TBL and MJB2 proteins, respectively (Fig. 5A, Table S7). It has been reported that ATHB was found to be involved in plant growth, development and stress tolerance in Arabidopsis [75] and drought tolerance in maize [76]. Moreover, ZF-HD proteins have been reported to respond to various abiotic stresses including heat, drought, salt, and cold [77, 78]. These results provide an important understanding of GmZF-HD proteins towards the possible regulation of plant development and stress tolerance. Furthermore, Faraji et al., [79] reported that protein structures are closely related to gene functions, therefore, we performed 3D structure analysis of all GmZF-HD proteins and found almost alike protein 3D structures, indicating conservation and consistency with the motif and amino acids (Fig. 5B, Table S8). Gouet et al., [80] demonstrated that proteins with highly conserved domains and structures belong to the same family and may perform similar functions.

MicroRNAs (miRNAs) are short noncoding RNAs (18–24 nucleotides) that are involved in multiple metabolic processes including plant growth, development, and stress responses [81,82,83]. In this study, we identified 94 miRNAs belonging to 41 different miRNA families and targeted 40 GmZF-HD genes (Fig. 6A-B, Table S9). Our study found that a large number of miRNAs, including gma-miR169, gma-miR9742, gma-miR4993, gma-miR1535, gma-miR4397, and gma-miR9752, that targeted most of the GmZF-HD genes, suggesting that they may play critical roles on these genes functions. Yu et al., [84] reported that overexpression of gma-miR169 enhanced drought tolerance in Arabidopsis, gma-miR9742 was involved in antiviral immunity [85], gma-miR4993 played a role in soybean isoflavonoid biosynthetic pathway [86], gma-miR1535 was involved in plant growth and development [87], gma-miR4397 involved in carpel and petal development [88] and gma-miR9752 induced nodule development [24]. Recently Zhou et al., [34] reported that miR166 targeted the CsHDZ3 gene and enhanced drought tolerance in tea. These findings indicated that the predicted miRNAs in soybean ZF-HD genes promoter regions may have certain significant importance that could play crucial roles in multiple processes including stresses, however further studies are needed.

Plant TFs facilitate gene expressions and play crucial roles in numerous biological processes including metabolism, signal transduction, and stress responses [89, 90]. In this study, we identified various TFs including ERF, LBD, BBR-BPC, MYB, bHLH, Dof, NAC and bZIP (Fig. 7, Table S10). Previous research reported that ERF enhanced osmotic stress tolerance in Arabidopsis [91] and Phytophthora sojae resistance in soybeans [92]. Yang et al., [93] investigated and found the role of LBD TF in abiotic stress. BBR-BPC TF was found to be involved in plant development, stress and immunity [94]. MYB TF was involved in isoflavonoid biosynthesis [95] and drought tolerance [96]. bHLH TF enhanced the salt tolerance [97]. Dof TF enhances drought and salt tolerance [98]. NAC TF enhanced drought [99] and salt [100] tolerance. bZIP TF was involved in biotic and abiotic stresses [101]. Furthermore, the results of GO annotation analysis identified various molecular, cellular and biological-related terms (Fig. 8, Table S11), and are consistent with previous studies on ZF-HD genes [39].

The expression profiling of all 51 GmZF-HD genes was carried out in different soybean organs/tissues using the publically available expression dataset. The results exhibited that several GmZF-HD genes showed distinct expression levels in specific tissues under different conditions (Fig. 9, Table S12). In the present study, several genes showed high expressions in different tissues, including GmZF-HD15 and GmZF-HD40 in roots under salt stress, GmZF-HD2/5/6/13/22 and GmZF-HD45 in shoots and GmZF-HD51 in leaves under drought stress conditions (Fig. 9, Table S12), indicating their significant roles under stress conditions. In previous studies, ZF-HD genes were also found to have obvious spatiotemporal expression differences under different abiotic stresses [39, 40, 69]. Khatun et al., [37] also reported that ZF-HD genes were differentially expressed in different tissues under different conditions and play a vital role in plant growth and development. Besides, ZF-HD genes have been extensively recognized for their key roles in various stresses, such as NtZF-HD21 was found to enhance drought stress tolerance in tobacco [102], AtZF-HD10 in Arabidopsis plays a synergistic role in regulating hormone signaling in response to stress [30]. SlZF-HD genes have been reported to be involved in fruit development and stress resistance in tomatoes [37]. It has been reported that overexpression of the CqZF-HD14 gene promoted enhanced drought tolerance in quinoa (Chenopodium quinoa) [103].

In addition, to gain insights into the identified GmZF-HD genes under various abiotic stress conditions, we performed qRT-PCR and results exhibited that most of the GmZF-HD genes responded to various stress treatments positively or negatively, for example, GmZF-HD5/6/13/39 and GmZF-HD45 genes were significantly upregulated in soybean leaves under different stress treatments compared with the control. Interestingly, GmZF-HD2/5 and GmZF-HD39 genes were upregulated under all treatments compared with control (Figs. 10, 11 and 12). Further, the GmZF-HD6/13/15/40/45 and GmZF-HD51 genes were upregulated in SA treatment, whereas, downregulated in ABA treatments. GmZF-HD15 and GmZF-HD34 genes were upregulated under heat, SA, MnSO4, ZnSO4 and FeSO4 treatments, while downregulated under CdCl2 and CoCl2 treatments (Figs. 10, 11 and 12).

Our findings are consistent with previous reports, such as Jing et al., [104] also found that the maize ZmZF-HD2 gene was upregulated under salt and drought treatments, while downregulated under ABA treatments. Additionally, the ZmZF-HD3 gene was upregulated under salt and ABA treatments, while downregulated under drought treatment. Zhou et al., [42] also found a similar expression trend in tea genes such as the CsZF-HD15 gene was upregulated under SA treatment, while downregulated under cold stress treatments. In addition, the expression patterns of GmZF-HD genes under heavy metal treatment were also consistent with the report of El-Sappah et al., [53] on soybean MTP genes, for example, GmaMTP2.2 was upregulated in leaves under CoCl2 and FeSO4 treatments, while downregulated under CdCl2, MnSO4 and ZnSO4 treatments. Overall, these findings highlighted the potential roles of GmZF-HD genes in response to stresses in soybeans and laid the foundation for further functional investigations.

Conclusion

In this study, for the first time, 51 GmZF-HD genes were comprehensively identified and characterized in the soybean genome. In addition, their physiochemical features, chromosomal location, phylogenetic evolution, gene structure, conserved domain and motifs, cis-elements, synteny and collinearity, protein-to-protein interaction, 3D modeling, putative miRNAs and TFs prediction and GO annotations analyses were performed to gain insight into the GmZF-HD genes. FPKM-based expression profiling of all GmZF-HD genes in diverse soybean tissues under diverse conditions was evaluated. qRT-PCR analysis highlighted the potential roles of GmZF-HD genes especially GmZF-HD5/6/13/39 and GmZF-HD45 genes in response to different stress treatments in soybean leaves. This comprehensive study provided valuable insights at the genomic level and highlighted the potential roles of GmZF-HD genes under different abiotic stresses and laying the foundation for further functional investigation of GmZF-HD genes for soybean improvement under adverse conditions.

Data availability

The authors declare that all the data and plant materials will be available without restrictions. All the plant materials and the raw data used in the manuscript will be available on request to Mingfu Wang (mfwang@szu.edu.cn).

References

  1. Goettel W, Zhang H, Li Y, Qiao Z, Jiang H, Hou D, Song Q, Pantalone VR, Song B-H, Yu D. POWR1 is a domestication gene pleiotropically regulating seed quality and yield in soybean. Nat Commun. 2022;13(1):3051.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  2. Gavili E, Moosavi AA, Haghighi AAK. Does biochar mitigate the adverse effects of drought on the agronomic traits and yield components of soybean? Ind Crops Prod. 2019;128:445–54.

    Article  CAS  Google Scholar 

  3. Shaffique S, Hussain S, Kang S-M, Imran M, Kwon E-H, Khan MA, Lee I-J. Recent progress on the microbial mitigation of heavy metal stress in soybean: overview and implications. Front Plant Sci. 2023;14:1188856.

    Article  PubMed  PubMed Central  Google Scholar 

  4. Deshmukh R, Sonah H, Patil G, Chen W, Prince S, Mutava R, Vuong T, Valliyodan B, Nguyen HT. Integrating omic approaches for abiotic stress tolerance in soybean. Front Plant Sci. 2014;5:244.

    Article  PubMed  PubMed Central  Google Scholar 

  5. Raza A, Bashir S, Khare T, Karikari B, Copeland RG, Jamla M, Abbas S, Charagh S, Nayak SN, Djalovic I. Temperature-smart plants: A new horizon with omics-driven plant breeding. Physiol Plant. 2024;176(1):e14188.

    Article  CAS  Google Scholar 

  6. El-Sappah AH, Rather SA. Genomics approaches to study abiotic stress tolerance in plants. In: Plant abiotic stress physiology. Apple Academic Press; 2022:25–46. ISBN9781003180579.

  7. Liu J, Gao Y, Tang Y, Wang D, Chen X, Yao Y, Guo Y. Genome-wide identification, comprehensive gene feature, evolution, and expression analysis of plant metal tolerance proteins in tobacco under heavy metal toxicity. Front Genet. 2019;10:345.

    Article  PubMed  PubMed Central  Google Scholar 

  8. Raza A, Razzaq A, Mehmood SS, Zou X, Zhang X, Lv Y, Xu J. Impact of climate change on crops adaptation and strategies to tackle its outcome: A review. Plants. 2019;8(2):34.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  9. Bukhari SAH, Peerzada AM, Javed MH, Dawood M, Hussain N, Ahmad S. Growth and development dynamics in agronomic crops under environmental stress. Agronomic Crops: Volume 1: Production Technologies. 2019:83–114.

  10. Raza A, Bhardwaj S, Rahman MA, García-Caparrós P, Copeland RG, Charagh S, Rivero RM, Gopalakrishnan S, Corpas FJ, Siddique KH. Fighting to thrive via plant growth regulators: Green chemical strategies for drought stress tolerance. Physiol Plant. 2024;176(6):e14605.

    Article  PubMed  CAS  Google Scholar 

  11. Khan S-A, Li M-Z, Wang S-M, Yin H-J. Revisiting the role of plant transcription factors in the battle against abiotic stress. Int J Mol Sci. 2018;19(6):1634.

    Article  PubMed  PubMed Central  Google Scholar 

  12. Huang C, Zhou J, Jie Y, Xing H, Zhong Y, Yu W, She W, Ma Y, Liu Z, Zhang Y. A ramie bZIP transcription factor BnbZIP2 is involved in drought, salt, and heavy metal stress response. DNA cell biology. 2016;35(12):776–86.

    Article  PubMed  CAS  Google Scholar 

  13. Chen H, Lai Z, Shi J, Xiao Y, Chen Z, Xu X. Roles of Arabidopsis WRKY18, WRKY40 and WRKY60 transcription factors in plant responses to abscisic acid and abiotic stress. BMC Plant Biol. 2010;10:1–15.

    Article  CAS  Google Scholar 

  14. Ma Z, Hu L. WRKY Transcription Factor Responses and Tolerance to Abiotic Stresses in Plants. Int J Mol Sci. 2024;25(13):6845.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  15. Tang W, Charles TM, Newton RJ. Overexpression of the pepper transcription factor CaPF1 in transgenic Virginia pine (Pinus virginiana Mill.) confers multiple stress tolerance and enhances organ growth. Plant Mol Biol. 2005;59:603–17.

    Article  PubMed  CAS  Google Scholar 

  16. Meraj TA, Fu J, Raza MA, Zhu C, Shen Q, Xu D, Wang Q. Transcriptional factors regulate plant stress responses through mediating secondary metabolism. Genes. 2020;11(4):346.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  17. Wang C, Qiao F, Wang M, Wang Y, Xu Y, Qi X. PvERF104 confers cadmium tolerance in Arabidopsis: Evidence for metal-responsive element-binding transcription factors. Environ Exp Bot. 2023;206: 105167.

    Article  CAS  Google Scholar 

  18. Tran LSP, Nakashima K, Sakuma Y, Osakabe Y, Qin F, Simpson SD, Maruyama K, Fujita Y, Shinozaki K, Yamaguchi-Shinozaki K. Co-expression of the stress-inducible zinc finger homeodomain ZFHD1 and NAC transcription factors enhances expression of the ERD1 gene in Arabidopsis. Plant J. 2007;49(1):46–63.

    Article  PubMed  CAS  Google Scholar 

  19. Hu J, Gao Y, Zhao T, Li J, Yao M, Xu X. Genome-wide identification and expression pattern analysis of zinc-finger homeodomain transcription factors in tomato under abiotic stress. J Am Soc Hortic Sci. 2018;143(1):14–22.

    Article  CAS  Google Scholar 

  20. Liu M-D, Liu H, Liu W-Y, Ni S-F, Wang Z-Y, Geng Z-H, Zhu K-Y, Wang Y-F, Zhao Y-H. Systematic Analysis of Zinc Finger-Homeodomain Transcription Factors (ZF-HDs) in Barley (Hordeum vulgare L.). Genes. 2024;15(5):578.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  21. Windhövel A, Hein I, Dabrowa R, Stockhaus J. Characterization of a novel class of plant homeodomain proteins that bind to the C4 phosphoenolpyruvate carboxylase gene of Flaveria trinervia. Plant Mol Biol. 2001;45:201–14.

    Article  PubMed  Google Scholar 

  22. Wang W, Wu P, Li Y, Hou X. Genome-wide analysis and expression patterns of ZF-HD transcription factors under different developmental tissues and abiotic stresses in Chinese cabbage. Mol Genet Genomics. 2016;291:1451–64.

    Article  PubMed  CAS  Google Scholar 

  23. Ariel FD, Manavella PA, Dezar CA, Chan RL. The true story of the HD-Zip family. Trends Plant Sci. 2007;12(9):419–26.

    Article  PubMed  CAS  Google Scholar 

  24. Piya S, Lopes-Caitar VS, Kim WS, Pantalone V, Krishnan HB, Hewezi T. JFiMB: Hypermethylation of miRNA genes during nodule development. Front Mol Biosci. 2021;8:616623.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  25. Tan QK-G, Irish VF. The Arabidopsis zinc finger-homeodomain genes encode proteins with unique biochemical properties that are coordinately expressed during floral development. Plant physiology. 2006;140(3):1095–108.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  26. Hu W, DePamphilis CW, Ma H. Phylogenetic analysis of the plant-specific zinc finger-homeobox and mini zinc finger gene families. J Integr Plant Biol. 2008;50(8):1031–45.

    Article  PubMed  CAS  Google Scholar 

  27. Xu Y, Wang Y, Long Q, Huang J, Wang Y, Zhou K, Zheng M, Sun J, Chen H, Chen SJP. Overexpression of OsZHD1, a zinc finger homeodomain class homeobox transcription factor, induces abaxially curled and drooping leaf in rice. Planta. 2014;239:803–16.

    Article  PubMed  CAS  Google Scholar 

  28. Hong S-Y, Kim O-K, Kim S-G, Yang M-S, Park C-M. Nuclear import and DNA binding of the ZHD5 transcription factor is modulated by a competitive peptide inhibitor in Arabidopsis. J Biol Chem. 2011;286(2):1659–68.

    Article  PubMed  CAS  Google Scholar 

  29. Bueso E, Muñoz-Bertomeu J, Campos F, Brunaud V, Martínez L, Sayas E, Ballester P, Yenush L, Serrano R. ARABIDOPSIS THALIANA HOMEOBOX25 uncovers a role for gibberellins in seed longevity. Plant Physiol. 2014;164(2):999–1010.

    Article  PubMed  CAS  Google Scholar 

  30. Perrella G, Davidson ML, O’Donnell L, Nastase A-M, Herzyk P, Breton G, Pruneda-Paz JL, Kay SA, Chory J, Kaiserli E. ZINC-FINGER interactions mediate transcriptional regulation of hypocotyl growth in Arabidopsis. Proc Natl Acad Sci. 2018;115(19):E4503–11.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  31. Park HC, Kim ML, Lee SM, Bahk JD, Yun D-J, Lim CO, Hong JC, Lee SY, Cho MJ, Chung WS. Pathogen-induced binding of the soybean zinc finger homeodomain proteins GmZF-HD1 and GmZF-HD2 to two repeats of ATTA homeodomain binding site in the calmodulin isoform 4 (GmCaM4) promoter. Nucleic Acids Res. 2007;35(11):3612–23.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  32. Li J, Li M, Shen T, Guo Q, Zhang R, Chen Y, Zhang Y, Luo K. Molecular and characterization of cassava zinc finger-homeodomain (ZF-HD) transcription factors reveals their role in disease resistance. Int J Biol Macromol. 2024;279:134846.

    Article  PubMed  CAS  Google Scholar 

  33. Chen C-y, Yang X-m, Chen S-y, Chen B, Yue L. Expression Analysis of the ZF-HD Gene Family in Chrysanthemum nankingense Under Drought and ABA Treatment. Biotechnol Bull. 2023;39(11):270.

    Google Scholar 

  34. Zhou C, Yang N, Tian C, Wen S, Zhang C, Zheng A, Hu X, Fang J, Zhang Z, Lai Z. The miR166 targets CsHDZ3 genes to negatively regulate drought tolerance in tea plant (Camellia sinensis). Int J Biol Macromol. 2024;264:130735.

    Article  PubMed  CAS  Google Scholar 

  35. Figueiredo DD, Barros PM, Cordeiro AM, Serra TS, Lourenço T, Chander S, Oliveira MM, Saibo NJ. Seven zinc-finger transcription factors are novel regulators of the stress responsive gene OsDREB1B. J Exp Bot. 2012;63(10):3643–56.

    Article  PubMed  CAS  Google Scholar 

  36. Yong Y, Zhang Y, Lyu YJP. Functional characterization of Lilium lancifolium cold-responsive Zinc Finger Homeodomain (ZFHD) gene in abscisic acid and osmotic stress tolerance. PeerJ. 2021;9:e11508.

    Article  PubMed  PubMed Central  Google Scholar 

  37. Khatun K, Nath UK, Robin AHK, Park J-I, Lee D-J, Kim M-B, Kim CK, Lim K-B, Nou IS, Chung M-Y. Genome-wide analysis and expression profiling of zinc finger homeodomain (ZHD) family genes reveal likely roles in organ development and stress responses in tomato. BMC Genomics. 2017;18:1–16.

    Article  Google Scholar 

  38. Saeid A-R, Khaldoun A-H. Novel zinc finger-homeodomain gene from barley (HvZFHD1) is differentially regulated during spike development and under hormonal treatments and abiotic stresses. Notulae Botanicae Horti Agrobotanici Cluj-Napoca. 2017;45(1):89–96.

    Article  Google Scholar 

  39. Liu H, Yang Y, Zhang LJP. Zinc finger-homeodomain transcriptional factors (ZF-HDs) in wheat (Triticum aestivum L.): identification, evolution, expression analysis and response to abiotic stresses. Plants. 2021;10(3):593.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  40. Islam MAU, Nupur JA, Shafiq M, Ali Q, Sami A, Shahid MA. In silico and computational analysis of zinc finger motif-associated homeodomain (ZF-HD) family genes in chilli (Capsicum annuum L). BMC Genomics. 2023;24(1):603.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  41. Zheng X-b, Wu Y, Wang H, Song S-w, Bai T-h, Jiao J, Song C-h, Pang H-g, Wang M-m. Genome-wide investigation of the zinc finger-homeodomain family genes reveals potential roles in apple fruit ripening. Front Genet. 2022;12:783482.

    Article  PubMed  PubMed Central  Google Scholar 

  42. Zhou C, Zhu C, Xie S, Weng J, Lin Y, Lai Z, Guo Y. Genome-wide analysis of zinc finger motif-associated homeodomain (ZF-HD) family genes and their expression profiles under abiotic stresses and phytohormones stimuli in tea plants (Camellia sinensis). Sci Hortic. 2021;281:109976.

    Article  CAS  Google Scholar 

  43. Abdullah M, Cheng X, Cao Y, Su X, Manzoor MA, Gao J, Cai Y, Lin Y. Zinc finger-homeodomain transcriptional factors (ZHDs) in upland cotton (Gossypium hirsutum): Genome-wide identification and expression analysis in fiber development. Front Genet. 2018;9:357.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  44. Li J, Li M, Shen T, Guo Q, Zhang R, Chen Y, Zhang Y, Luo K. Molecular characterization of cassava zinc finger-homeodomain (ZF-HD) transcription factors reveals their role in disease resistance. Int J Biol Macromol. 2024;279:134846.

    Article  PubMed  CAS  Google Scholar 

  45. Schmutz J, Cannon SB, Schlueter J, Ma J, Mitros T, Nelson W, Hyten DL, Song Q, Thelen JJ, Cheng JJ. Genome sequence of the palaeopolyploid soybean. Nature. 2010;463(7278):178–83.

    Article  PubMed  CAS  Google Scholar 

  46. Chen C, Chen H, Zhang Y, Thomas HR, Frank MH, He Y, Xia R. TBtools: an integrative toolkit developed for interactive analyses of big biological data. Mol Plant. 2020;13(8):1194–202.

    Article  PubMed  CAS  Google Scholar 

  47. Rizwan HM, He J, Nawaz M, Cheng K-W, Wang M. Comprehensive in silico characterization of soybean (Glycine max L.) isoflavone reductase genes and their expressions in response to spermidine and ultrasonication. Plant Stress. 2024;11:100392.

    Article  CAS  Google Scholar 

  48. Fehr WR, Caviness CE. Stages of soybean development. Iowa: Iowa State University, Ames; 1977. p. 11.

    Google Scholar 

  49. Gong Y, Wang D, Xie H, Zhao Z, Chen Y, Zhang D, Jiao Y, Shi M, Lv P, Sha Q. Genome-wide identification and expression analysis of the KCS gene family in soybean (Glycine max) reveal their potential roles in response to abiotic stress. Front Plant Sci. 2023;14:1291731.

    Article  PubMed  PubMed Central  Google Scholar 

  50. Huai D, Xue X, Li Y, Wang P, Li J, Yan L, Chen Y, Wang X, Liu N, Kang Y. Genome-wide identification of peanut KCS genes reveals that AhKCS1 and AhKCS28 are involved in regulating VLCFA contents in seeds. Front Plant Sci. 2020;11:406.

    Article  PubMed  PubMed Central  Google Scholar 

  51. Wang S, Guo J, Zhang Y, Guo Y, Ji W. Genomics: Genome-wide characterization and expression analysis of TOPP-type protein phosphatases in soybean (Glycine max L.) reveal the role of GmTOPP13 in drought tolerance. Genes Genomics. 2021;43:783–96.

    Article  PubMed  CAS  Google Scholar 

  52. Zhang B, Liu Z, Zhou R, Cheng P, Li H, Wang Z, Liu Y, Li M, Zhao Z, Hu Z. Genome-wide analysis of soybean DnaJA-family genes and functional characterization of GmDnaJA6 responses to saline and alkaline stress. Crop J. 2023;11(4):1230–41.

    Article  CAS  Google Scholar 

  53. El-Sappah AH, Abbas M, Rather SA, Wani SH, Soaud N, Noor Z, Qiulan H, Eldomiaty AS, Mir RR, Li J. Genome-wide identification and expression analysis of metal tolerance protein (MTP) gene family in soybean (Glycine max) under heavy metal stress. Mol Biol Rep. 2023;50(4):2975–90.

    Article  PubMed  CAS  Google Scholar 

  54. Schmittgen TD, Livak K. Analyzing real-time PCR data by the comparative CT method. Nat Protoc. 2008;3(6):1101–8.

    Article  PubMed  CAS  Google Scholar 

  55. Dilawari R, Kaur N, Priyadarshi N, Prakash I, Patra A, Mehta S, Singh B, Jain P, Islam MA. Soybean: A key player for global food security. In: Soybean Improvement: Physiological, Molecular and Genetic Perspectives. Cham: Springer International Publishing; 2022: 1–46.

  56. Ren C, Wang H, Zhou Z, Jia J, Zhang Q, Liang C, Li W, Zhang Y, Yu G. Genome-wide identification of the B3 gene family in soybean and the response to melatonin under cold stress. Front Plant Sci. 2023;13:1091907.

    Article  PubMed  PubMed Central  Google Scholar 

  57. Mazarei M, Routray P, Piya S, Stewart CN Jr, Hewezi T. Overexpression of soybean GmNAC19 and GmGRAB1 enhances root growth and water-deficit stress tolerance in soybean. Front Plant Sci. 2023;14:1186292.

    Article  PubMed  PubMed Central  Google Scholar 

  58. Wang L, Zhang J, Li H, Zhang G, Hu D, Zhang D, Xu X, Yang Y, Huang Z. Genome-Wide Identification of the Phytocyanin Gene Family and Its Potential Function in Salt Stress in Soybean (Glycine max (L.) Merr.). Agronomy. 2023;13(10):2484.

    Article  CAS  Google Scholar 

  59. Sabagh AE, Hossain A, Islam MS, Iqbal MA, Fahad S, Ratnasekera D, Llanes A. Consequences and mitigation strategies of heat stress for sustainability of soybean (Glycine max L. Merr.) production under the changing climate. Plant Stress Physiol. 2020;1 (7):119-40.

  60. Peña Calzada K, Olivera Viciedo D, Habermann E, Calero Hurtado A, Lupino Gratão P, De Mello PR, Lata-Tenesaca LF, Martinez CA, Ajila Celi GE, Rodríguez JCJA. Exogenous application of amino acids mitigates the deleterious effects of salt stress on soybean plants. Agronomy. 2022;12(9):2014.

    Article  Google Scholar 

  61. Hussain A, Shah M, Hamayun M, Iqbal A, Qadir M, Alataway A, Dewidar AZ, Elansary HO, Lee I-J. Phytohormones producing rhizobacteria alleviate heavy metals stress in soybean through multilayered response. Microbiol Res. 2023;266:127237.

    Article  PubMed  Google Scholar 

  62. Bisht A, Saini DK, Kaur B, Batra R, Kaur S, Kaur I, Jindal S, Malik P, Sandhu PK, Kaur A. Multi-omics assisted breeding for biotic stress resistance in soybean. Mol Biol Rep. 2023;50(4):3787–814.

    Article  PubMed  CAS  Google Scholar 

  63. Neale DB, Martinez-Garcia PJ, De La Torre AR, Montanari S, Wei X-X. Novel insights into tree biology and genome evolution as revealed through genomics. Annu Rev Plant Biol. 2017;68:457–83.

    Article  PubMed  CAS  Google Scholar 

  64. Segata N, Huttenhower C. Toward an efficient method of identifying core genes for evolutionary and functional microbial phylogenies. PLoS ONE. 2011;6(9):e24704.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  65. Wang L, Hua D, He J, Duan Y, Chen Z, Hong X, Gong Z. Auxin Response Factor2 (ARF2) and its regulated homeodomain gene HB33 mediate abscisic acid response in Arabidopsis. PLoS Genet. 2011;7(7):e1002172.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  66. Feng R-J, Ren M-Y, Lu L-F, Peng M, Guan X, Zhou D-B, Zhang M-Y, Qi D-F, Li K, Tang W. Involvement of abscisic acid-responsive element-binding factors in cassava (Manihot esculenta) dehydration stress response. Sci Rep. 2019;9(1):12661.

    Article  PubMed  PubMed Central  Google Scholar 

  67. Shariatipour N, Heidari B. Meta-analysis of expression of the stress tolerance associated genes and uncover their-regulatory elements in rice (L.).  Open Bioinform J. 2020;13(1):13-39.

  68. Qiao X, Yin H, Li L, Wang R, Wu J, Wu J, Zhang S. Different modes of gene duplication show divergent evolutionary patterns and contribute differently to the expansion of gene families involved in important fruit traits in pear (Pyrus bretschneideri). Front Plant Sci. 2018;9:161.

    Article  PubMed  PubMed Central  Google Scholar 

  69. Wang X, Sun W, Ma Z, Zheng T, Huang L, Wu Q, Tang Z, Bu T, Li C, Chen H. Genome-wide investigation of the ZF-HD gene family in Tartary buckwheat (Fagopyrum tataricum). BMC Plant Biol. 2019;19(1):248.

    Article  PubMed  PubMed Central  Google Scholar 

  70. Li J, Li M, Shen T, Guo Q, Zhang R, Chen Y, Zhang Y, Luo K. Molecular characterization of cassava zinc finger-homeodomain (ZF-HD) transcription factors reveals their role in disease resistance. Int J Biol Macromol. 2024;279:134846.

    Article  PubMed  CAS  Google Scholar 

  71. Shan Q, Zhou B, Wang Y, Hao F, Zhu L, Liu Y, Wang N, Wang F, Li X, Dong Y. Genome-wide identification and comprehensive analysis of the Ftsh gene family in soybean (Glycine max). Int J Mol Sci. 2023;24(23):16996.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  72. Wang L, Ding X, Gao Y, Yang S. Genome-wide identification and characterization of GRAS genes in soybean (Glycine max). BMC Plant Biol. 2020;20:1–21.

    Article  Google Scholar 

  73. Rao VS, Srinivas K, Sujini G, Kumar GS. Protein-protein interaction detection: methods and analysis. Int J Proteomics. 2014;2014(1):147648.

    PubMed  PubMed Central  Google Scholar 

  74. Park J, Son A, Kim HJSR. A protein–protein interaction analysis tool for targeted cross-linking mass spectrometry. Sci Rep. 2023;13(1):22103.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  75. Carabelli M, Sessa G, Baima S, Morelli G, Ruberti I. The Arabidopsis Athb-2 and-4 genes are strongly induced by far-red-rich light. Plant J. 1993;4(3):469–79.

    Article  PubMed  CAS  Google Scholar 

  76. Jiao P, Jiang Z, Wei X, Liu S, Qu J, Guan S, Ma Y. Overexpression of the homeobox-leucine zipper protein ATHB-6 improves the drought tolerance of maize (Zea mays L.). Plant Sci. 2022;316:111159.

    Article  PubMed  CAS  Google Scholar 

  77. Shalmani A, Muhammad I, Sharif R, Zhao C, Ullah U, Zhang D, Jing X-Q, Amin B, Jia P, Mobeen Tahir MJEB. Zinc finger-homeodomain genes: Evolution, functional differentiation, and expression profiling under flowering-related treatments and abiotic stresses in plants. Evol Bioinforma. 2019;15:1176934319867930.

    Article  Google Scholar 

  78. Zhang Y, Gao P, Yuan JS. Plant protein-protein interaction network and interactome. Curr Genomics. 2010;11(1):40–6.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  79. Faraji S, Rasouli SH, Kazemitabar SK. Genome-wide exploration of C2H2 zinc finger family in durum wheat (Triticum turgidum ssp. Durum): insights into the roles in biological processes especially stress response. Biometals. 2018;31(6):1019–42.

    Article  PubMed  CAS  Google Scholar 

  80. Gouet P, Robert X, Courcelle E. ESPript/ENDscript: extracting and rendering sequence and 3D information from atomic structures of proteins. Nucleic Acids Res. 2003;31(13):3320–3.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  81. Kumar R. biotechnology: Role of microRNAs in biotic and abiotic stress responses in crop plants. Appl Biochem Biotechnol. 2014;174:93–115.

    Article  PubMed  CAS  Google Scholar 

  82. Zhang B. MicroRNA: a new target for improving plant tolerance to abiotic stress. J Exp Bot. 2015;66(7):1749–61.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  83. Sun X, Lin L, Sui N. Regulation mechanism of microRNA in plant response to abiotic stress and breeding. Mol Biol Rep. 2019;46:1447–57.

    Article  PubMed  CAS  Google Scholar 

  84. Yu Y, Ni Z, Wang Y, Wan H, Hu Z, Jiang Q, Sun X, Zhang H. Overexpression of soybean miR169c confers increased drought stress sensitivity in transgenic Arabidopsis thaliana. Plant Sci. 2019;285:68–78.

    Article  PubMed  CAS  Google Scholar 

  85. Ramesh SV, Gupta GK, Husain SM. Soybean (Glycine max) microRNAs display proclivity to repress begomovirus genomes. Curr Sci. 2016;110 (3):424-8.

  86. Gupta OP, Dahuja A, Sachdev A, Kumari S, Jain PK, Vinutha T, Praveen S. Conserved miRNAs modulate the expression of potential transcription factors of isoflavonoid biosynthetic pathway in soybean seeds. Mol Biol Rep. 2019;46:3713–30.

    Article  PubMed  CAS  Google Scholar 

  87. Xu S, Liu N, Mao W, Hu Q, Wang G, Gong Y. Identification of chilling-responsive microRNAs and their targets in vegetable soybean (Glycine max L.). Sci Rep. 2016;6(1):26619.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  88. Kulcheski FR, Molina LG, da Fonseca GC, de Morais GL, de Oliveira LF, Margis RJG. Novel and conserved microRNAs in soybean floral whorls. Genes. 2016;575(2):213–23.

    CAS  Google Scholar 

  89. Huang Y, An J, Sircar S, Bergis C, Lopes CD, He X, Da Costa B, Tan F-Q, Bazin J, Antunez-Sanchez J. HSFA1a modulates plant heat stress responses and alters the 3D chromatin organization of enhancer-promoter interactions. Nat Commun. 2023;14(1):469.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  90. Yuan HY, Kagale S, Ferrie AM. Multifaceted roles of transcription factors during plant embryogenesis. Front Plant Sci. 2024;14:1322728.

    Article  PubMed  PubMed Central  Google Scholar 

  91. Zhao M-J, Yin L-J, Liu Y, Ma J, Zheng J-C, Lan J-H, Fu J-D, Chen M, Xu Z-S, Ma Y-Z. The ABA-induced soybean ERF transcription factor gene GmERF75 plays a role in enhancing osmotic stress tolerance in Arabidopsis and soybean. BMC Plant Biol. 2019;19:1–14.

    Article  Google Scholar 

  92. Dong L, Cheng Y, Wu J, Cheng Q, Li W, Fan S, Jiang L, Xu Z, Kong F, Zhang D. Overexpression of GmERF5, a new member of the soybean EAR motif-containing ERF transcription factor, enhances resistance to Phytophthora sojae in soybean. J Exp Bot. 2015;66(9):2635–47.

    Article  PubMed  CAS  Google Scholar 

  93. Yang H, Shi G, Du H, Wang H, Zhang Z, Hu D, Wang J, Huang F, Yu D. Genome-Wide Analysis of Soybean LATERAL ORGAN BOUNDARIES Domain-Containing Genes: A Functional Investigation of GmLBD12. Plant Genome. 2017;10(1):plantgenome 2016.2007.0058.

    Article  Google Scholar 

  94. Sahu A, Singh R, Verma PKJP. Plant BBR/BPC transcription factors: unlocking multilayered regulation in development, stress and immunity. Planta. 2023;258(2):31.

    Article  PubMed  CAS  Google Scholar 

  95. Yi J, Derynck MR, Li X, Telmer P, Marsolais F, Dhaubhadel S. A single-repeat MYB transcription factor, GmMYB176, regulates CHS8 gene expression and affects isoflavonoid biosynthesis in soybean. Plant J. 2010;62(6):1019–34.

    PubMed  CAS  Google Scholar 

  96. Wang N, Zhang W, Qin M, Li S, Qiao M, Liu Z, Xiang F. Physiology C: Drought tolerance conferred in soybean (Glycine max. L) by GmMYB84, a novel R2R3-MYB transcription factor. Plant Cell Physiol. 2017;58(10):1764–76.

    Article  PubMed  CAS  Google Scholar 

  97. Liu X, Pi B, Du Z, Yang T, Gu M, Sun S, Yu BJE, Botany E. The transcription factor GmbHLH3 confers Cl−/salt tolerance to soybean by upregulating GmCLC1 expression for maintenance of anion homeostasis. Environ Exp Bot. 2022;194:104755.

    Article  CAS  Google Scholar 

  98. Wei J-T, Zhao S-P, Zhang H-Y, Jin L-G, Yu T-F, Zheng L, Ma J, Chen J, Zhou Y-B, Chen M. GmDof41 regulated by the DREB1-type protein improves drought and salt tolerance by regulating the DREB2-type protein in soybean. Int J Biol Macromol. 2023;230:123255.

    Article  PubMed  CAS  Google Scholar 

  99. Yang C, Huang Y, Lv P, Antwi-Boasiako A, Begum N, Zhao T, Zhao J. NAC transcription factor GmNAC12 improved drought stress tolerance in soybean. Int J Mol Sci. 2022;23(19):12029.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  100. Li M, Chen R, Jiang Q, Sun X, Zhang H, Hu Z. GmNAC06, a NAC domain transcription factor enhances salt stress tolerance in soybean. Plant Mol Biol. 2021;105:333–45.

    Article  PubMed  CAS  Google Scholar 

  101. Chai M, Fan R, Huang Y, Jiang X, Wai MH, Yang Q, Su H, Liu K, Ma S, Chen Z. GmbZIP152, a soybean bZIP transcription factor, confers multiple biotic and abiotic stress responses in plant. Int J Mol Sci. 2022;23(18):10935.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  102. Sun J, Xie M, Li X, Li Z, Wang Q, Ding A, Wang W, Sun YJA. Systematic investigations of the ZF-HD gene family in tobacco reveal their multiple roles in abiotic stresses. Agronomy. 2021;11(3):406.

    Article  CAS  Google Scholar 

  103. Sun W, Wei J, Wu G, Xu H, Chen Y, Yao M, Zhan J, Yan J, Wu N, Chen HJPS. CqZF-HD14 enhances drought tolerance in quinoa seedlings through interaction with CqHIPP34 and CqNAC79. Plant Sci. 2022;323: 111406.

    Article  PubMed  CAS  Google Scholar 

  104. Jing X, Li C, Luo C, Yao C, Zhang J, Zhu T, Wang J, Liu C. Identification and characterization of ZF-HD genes in response to abscisic acid and abiotic stresses in maize. Phyton Int J Exp Bot. 2022;92(3):707–23.

    Google Scholar 

Download references

Acknowledgements

The authors extend their gratitude to the College of Civil and Transportation Engineering, College of Chemistry and Environmental Engineering, Institute for Advanced Study, and Mingfu Wang, Shenzhen University, China for providing the experiment environment and funding support.

Funding

This work was supported by grants from the Shenzhen Science and Technology Program (ZDSYS20220117155800001).

Author information

Authors and Affiliations

  1. College of Civil and Transportation Engineering, Shenzhen University, Shenzhen, 518060, China

    Hafiz Muhammad Rizwan & Mingfu Wang

  2. Shenzhen Key Laboratory of Food Nutrition and Health, College of Chemistry and Environmental Engineering, Shenzhen University, Shenzhen, 518060, China

    Hafiz Muhammad Rizwan, Jiayi He & Mingfu Wang

  3. Institute for Advanced Study, Shenzhen University, Shenzhen, 518060, China

    Muhammad Nawaz & Keyu Lu

Authors
  1. Hafiz Muhammad Rizwan
    View author publications

    Search author on:PubMed Google Scholar

  2. Jiayi He
    View author publications

    Search author on:PubMed Google Scholar

  3. Muhammad Nawaz
    View author publications

    Search author on:PubMed Google Scholar

  4. Keyu Lu
    View author publications

    Search author on:PubMed Google Scholar

  5. Mingfu Wang
    View author publications

    Search author on:PubMed Google Scholar

Contributions

Hafiz Muhammad Rizwan and Mingfu Wang conceived the research plan. Mingfu Wang supervised the entire experiment. Hafiz Muhammad Rizwan organized and conducted the whole experiment. Jiayi He, Muhammad Nawaz and Keyu Lu provided technical support during the experiment. Hafiz Muhammad Rizwan analyzed the data and wrote the first draft of the manuscript. Mingfu Wang reviewed and edited the final version of the manuscript. Mingfu Wang contributed to the funding and correspondence. All authors have read and agreed to the published version of the manuscript.

Corresponding author

Correspondence to Mingfu Wang.

Ethics declarations

Ethics approval and consent to participate

All experimental research on plants, including the collection of plant materials has been used by following the relevant institutional, national, and international guidelines and legislation.

Competing interests

The authors declare no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

12870_2024_6028_MOESM1_ESM.xlsx

Additional file 1: Table S1. Primers used in this study. Table S2. Protein sequences were used in this study. Table S3. GmZF-HD protein percent similarity index. Table S4. Cis-regulatory element analysis. Table S5. Synteny analysis detailed information. Table S6. Multi-collinearity analysis detailed information. Table S7. Protein-to-protein interaction analysis detailed information. Table S8. Protein 3D structure information. Table S9. Putative predicted miRNA information. Table S10. Predicted transcription factors in upstream regions of GmZF-HD genes. Table S11. Gene ontology analysis. Table S12. Gene expression analysis information.

Additional file 2: Fig S1. Motifs logo of all the 51 GmZF-HD proteins.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Rizwan, H.M., He, J., Nawaz, M. et al. The members of zinc finger-homeodomain (ZF-HD) transcription factors are associated with abiotic stresses in soybean: insights from genomics and expression analysis. BMC Plant Biol 25, 56 (2025). https://doi.org/10.1186/s12870-024-06028-x

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12870-024-06028-x

Keywords

BMC Plant Biology

ISSN: 1471-2229

Contact us