Genome-Wide Identification and Exogenous Hormone and Stress Response Expression Analysis of the GARP Gene Family in Soybean (Glycine max)

Abstract

The GARP transcription factor family is integral to the regulation of plant growth, development, hormone signaling pathways, circadian rhythms, and responses to both biotic and abiotic stressors. Despite its recognized importance, a comprehensive characterization of the GARP gene family in Glycine max remains unexplored. In this study, we identified 126 GmGARP genes across the 16 chromosomes of G. max and elucidated their diverse physicochemical properties. Phylogenetic classification grouped the GmGARP genes into eight distinct subfamilies, based on conserved motifs and gene structures, suggesting functional and evolutionary conservation within these clusters. The discovery of 56 segmentally duplicated gene pairs highlights gene duplication as a key driver of family expansion. Promoter analysis revealed various cis-regulatory elements, while expression profiling demonstrated that GmGARP genes possess significant tissue specificity. Furthermore, qRT-PCR analysis indicated that GmGARP genes are highly responsive to exogenous hormones, such as ABA, MeJA, and GA, as well as to abiotic stresses, including cold, salt, and drought conditions. Notably, GmGARP120 and GmGARP98 contain specific cis-elements linked to hormone responses, with their interaction verified through yeast two-hybrid (Y2H) and bimolecular fluorescence complementation (BiFC) assays. Additionally, 11,195 potential target genes were predicted, underscoring the regulatory potential of the GmGARP transcription factors. These findings provide significant insights into the GmGARP gene family, laying a strong foundation for future studies on its role in G. max development and adaptive responses to environmental stressors.

1. Introduction

Plants rely on the coordinated expression of numerous genes and transcription factors to regulate essential processes such as tissue growth, nutrient assimilation, metabolic activity, and signal transduction during development [1,2]. The concept of genes and transcription factors, while related, reflects distinct biological elements with unique structural and functional roles. Gene sequences encode proteins or functional RNAs, guiding their synthesis through transcription and translation [3]. In contrast, transcription factors are proteins encoded by these genes that regulate gene expression by binding to specific DNA sequences during transcription [4]. Various transcription factor families, including MYB [5], WRKY [6], AP2/ERF [7], and Basic Helix–Loop–Helix (bHLH) [8], have been extensively studied in plant science. Investigating transcription factor families is crucial in plant science, offering insights into their evolution, diversification, and regulatory roles in plant growth, development, and environmental adaptation [4,6]. Such in-depth studies can enhance our understanding of plant life processes and facilitate the development of crop varieties with superior traits, thereby boosting agricultural productivity [9].

The GARP transcription factor family, in particular, plays a central role in regulating a variety of processes, including growth, development, hormone regulation, circadian rhythms, and responses to biotic and abiotic stresses [10]. This gene family is named after the Golden 2 (G2) protein in Zea mays, the type-B authentic response regulator (ARR-B) protein in A. thaliana, and the phosphate starvation response 1 (PSR1) protein in Chlamydomonas [11]. These transcription factors share a conserved GARP DNA-binding domain, approximately 60 amino acids long, with a helix–turn–helix structure that binds to the promoter regions of target genes [11,12]. Members of this family, such as GLK transcription factors, are primarily involved in regulating chloroplast development and photosynthesis-related gene expression. In A. thaliana, GLK1 and GLK2 exhibit functional redundancy, with double mutants leading to impaired chloroplast development and reduced photosynthetic efficiency [13]. GLK factors promote chloroplast development and function by activating genes involved in chlorophyll synthesis and photosynthesis [14]. ARR-B transcription factors are essential in the cytokinin signaling pathway. Activated by cytokinins, they undergo phosphorylation and regulate cytokinin-responsive genes [15]. PSR1 transcription factors are crucial in response to phosphorus starvation, becoming activated under phosphorus-deficient conditions to regulate genes related to phosphorus uptake and utilization, aiding survival in low-phosphate environments [12]. In A. thaliana, genes such as ARR1, ARR2, ARR10, ARR11, and ARR12 exhibit significant cytokinin regulation functions [16]. Additionally, KAN1 in A. thaliana regulates responses related to abscisic acid, jasmonic acid, brassinosteroids, ethylene, cytokinin, and gibberellin [17]. KAN transcription factors also contribute to leaf and Oryza sativa spikelet development [18,19]. PHR1-Like 1 (PHL1) and PHL2 in A. thaliana, as well as OsPHR2 in Oryza sativa, regulate phosphorus transport [20,21]. In A. thaliana, the glk1 glk2 double knockout mutant shows decreased Light Harvesting Complex II (LHCII) protein accumulation, responsible for chlorophyll a/b binding, while RuBisCO levels remain unaffected [14]. Overexpression of GLKs in various tissues leads to chlorophyll accumulation and chloroplast formation; for instance, overexpression in roots results in root greening, and in fruits, it increases fruit greenness [22,23]. Overexpressing GLK1 enhances A. thaliana resistance to Fusarium graminearum [24]. AtGLKs are involved in jasmonic acid-dependent susceptibility to the biotrophic pathogen Hpa Noco2 and jasmonic acid-independent resistance to the necrotrophic pathogen Botrytis cinerea [25]. GLKs also regulate the antioxidant system in A. thaliana, enhancing resistance to cucumber mosaic virus [26]. In C. sinensis, CsGARP gene expression significantly changes under abiotic stresses like cold, drought, and salt [27]. In Gossypium hirsutum, silencing GhGLK1 leads to plant damage under drought and cold stress [28].

Soybean (Glycine max) is a vital oil crop with considerable economic importance, yet it is highly vulnerable to abiotic stress, which negatively affects both yield and quality. The GARP transcription factor family, characterized by its diversity and functional specificity, plays an essential role in plant growth, development, and environmental adaptation. By regulating various gene expression pathways, these transcription factors enable plants to maintain normal growth and development under a wide range of physiological and environmental conditions [29]. Glyma.05G033000, Glyma.19G132300 [30], Glyma.15G145200 [31], and Glyma.03G130000 [30] demonstrate significant expression changes under drought conditions, while Glyma.10G281000 [32] and Glyma.11G155100 [33] are highly responsive to salinity. Additionally, Glyma.19G047800 [34] has been strongly linked to cold tolerance. Despite these insights, the functions and stress response mechanisms of most GARP family genes remain largely unexplored, warranting further investigation to clarify their roles in soybean stress adaptation.

In this study, we aim to systematically identify and characterize members of the Glycine max GARP gene family through comprehensive analyses of their phylogeny, conserved motifs, gene structure, gene duplication events, cis-regulatory elements, and expression profiles. Additionally, we will examine the tissue-specific expression patterns of the GARP family and investigate their responses to various abiotic and hormonal stresses. Protein interactions within the GARP family will be explored using yeast two-hybrid (Y2H) and bimolecular fluorescence complementation (BiFC) assays. Together, these analyses will establish a foundation for further investigations into the functional roles of GmGARP genes in Glycine max growth, development, and stress adaptation, providing valuable insights into the mechanisms by which these genes contribute to enhancing Glycine max resilience to environmental stresses.

2. Materials and Methods

2.1. Identification of GARP Family Members in Glycine max Genome

Glycine max genomic data were obtained from the Ensembl Plants database [35]. To identify members of the GARP gene family, we utilized 56 A. thaliana GARP genes and their corresponding protein sequences, which were downloaded from the UniProt database. These protein sequences were then aligned with the Glycine max genome transcript and protein sequences using the BLAST tool [36], and we aligned the Glycine max genome transcript protein sequences with the A. thaliana GARP protein sequences. After eliminating redundant genes, we identified 126 GmGARP genes, which were retained based on an E-value threshold of ≤e⁻⁵. The domain structures of these genes were analyzed using the InterPro database [37], confirming that all genes belonged to the GARP transcription factor family. The physicochemical properties of the 126 GARP proteins were further analyzed using the ExPASy tool [38], and their subcellular localizations were predicted using the Plant-mPLoc tool [39].

2.2. Phylogeny, Motifs, and Gene Structures

To explore the evolutionary relationships and classify the GmGARP transcription factors, we constructed phylogenetic trees using MEGA X (v.11). Protein sequences were aligned with ClustalW, and the Neighbor-Joining (NJ) method was used for tree construction with 1000 bootstrap replications. The Poisson correction model and pairwise deletion were applied to refine the analysis [40]. The classification of GmGARP family members followed the established A. thaliana GARP classification system, and gene groupings were visualized using the iTOL online tool [41]. Motif identification was conducted using the MEME database, specifying 10 motifs with default parameters [42]. Gene structure analysis, including exon and intron counts, was performed using GFF3 data from the Glycine max genome, and the results were visualized using TBtools [43].

2.3. Gene Locations, Gene Duplications, and Ka/Ks

The chromosomal positions of the GmGARP genes were determined using Glycine max GFF3 data, and a chromosomal location map was generated using TBtools. Gene duplication events were identified with MCScanX, distinguishing between segmental and tandem duplication events, and a corresponding map of segmental duplication gene pairs was created using the Circos tool [44,45]. A. thaliana and Oryza sativa were selected for homology comparisons, and homologous GARP gene pairs between Glycine max, A. thaliana, and Oryza sativa were identified and mapped. Ka (non-synonymous substitution rate) and Ks (synonymous substitution rate) values for gene duplication events were calculated using the KaKs_Calculator 3.0 tool [46].

2.4. Cis-Acting Elements Annotation

The promoter regions (2000 bp upstream) corresponding to the GmGARP genes were extracted using a Perl script, and the cis-acting elements were annotated with the PlantCARE online tool [47]. The cis-elements related to hormone regulation, stress response, and growth/development were mapped to visualize their distributions within the promoter regions of GmGARP genes [48].

2.5. Expression Pattern

The expression levels (FPKM) of GmGARP family members were analyzed across nine Glycine max tissues (cotyledon, flower, leaf, pod, seed pod, root, seed, shoot meristem, and stem), were obtained from the Glycine max multi-omics database [49] (https://yanglab.hzau.edu.cn/SoyMD/#/) (accessed on 13 April 2024). Genes with transcripts per kilobase million (TPM) values greater than or equal to 0.8 were considered to exhibit significant expression. These values were analyzed, and heatmaps were generated using the Hiplot Pro online tool (https://hiplot.com.cn/home/index.html) (accessed on 20 April 2024).

2.6. Protein Interaction Network and Target Gene Prediction

Protein interactions among GmGARP family members were predicted using the A. thaliana STRING database (v.12) (https://string-db.org/) (accessed on 25 April 2024), and the predicted interactions were visualized using Cytoscape [50,51]. Binding sites for the GARP11 transcription factor were retrieved from the JASPAR Plantae database [52]. Subsequently, promoter sequences (2000 bp upstream) of all genes in the Glycine max genome were extracted, and genes binding to the GARP11 transcription factor were identified using the MEME suit (v.5.5.7) (https://meme-suite.org/meme/) [42]. Target gene domains were predicted using the PFAM database, and GO (Gene Ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes) enrichment analyses were conducted on the target genes using OmicShare Tools (https://www.omicshare.com/tools) (accessed on 14 May 2024).

2.7. Plant Materials and Stress Treatments

Glycine max cultivars were obtained from Heinong 84 (Heilongjiang Academy of Agricultural Sciences, Jiamusi, China) and cultivated at 25 °C under a 16 h light/8 h dark photoperiod. Stable seedling materials were subjected to stress treatments at the subsequent trifoliate leaves stage. Leaf samples were collected at intervals (0, 1, 2, 4, 8, 16, and 24 h) following low temperature (4 °C), NaCl (120 mM), PEG6000 (20%), and hormone (gibberellin, methyl jasmonate, abscisic acid) treatments.

2.8. RNA Extraction and Real-Time Fluorescence Quantitative PCR

Total RNA from Glycine max leaves was extracted using TRIzol reagent (Vazyme, Nanjing, China) and reverse transcribed to cDNA using a HiScript II 1st Strand cDNA Synthesis Kit (+gDNA wiper) (Vazyme). Real-time fluorescence quantitative PCR (RT-qPCR) was performed using a CFX384 Real-Time PCR detection system (Bio-Rad, Hercules, CA, USA) and ChamQ Universal SYBR^® qPCR Master Mix (Vazyme) for real-time fluorescence quantitative PCR (RT-qPCR). The relative expression levels of GmGARP were measured under exogenous hormone and abiotic stress treatments. Primer sequences are shown in Supplementary Table S1. Three biological replicates were used for each experiment, and results were calculated using the 2^−ΔΔCT method.

2.9. Yeast Two-Hybrid Assay

The coding region sequences of GmGARP120 and GmGARP98 were cloned, ligated into pGBKT7 and pGADT7 vectors, respectively, and transfected into the yeast strain Y2HGold (Anyu Bio, Shanghai, China). pGBKT7-Lam/pGADT7-T- and pGBKT7-53/pGADT7-T-co-transformed yeast cells were used as negative and positive control, respectively. Transfected cells were grown on SD/-Leu/-Trp/-His/-Ade/x-α-gal medium (ELITE-MEDIA, Shanghai, China) for interaction testing.

2.10. BiFC Verification

The coding sequence of GmGARP120, excluding the termination codon, was amplified and inserted into the pUC-nEYFP vector to construct the GmGARP120 fusion protein. Similarly, the coding sequence of GmGARP98 was PCR-amplified, digested, and ligated into the pUC-cEYFP vector to create the GmGARP98 fusion. These constructs were transiently expressed in Nicotiana benthamiana leaf epidermal cells. YFP fluorescence was then observed using confocal laser scanning microscopy, with excitation and emission wavelengths set to 488 nm and 507 nm, respectively, to confirm protein localization.

3. Results

3.1. Analysis of Protein Characteristics of GmGARP Family Members

We identified 126 GARP genes in Glycine max, which we designated as GmGARP1 to GmGARP126 based on their chromosomal positions. The physicochemical properties of the GmGARP proteins showed significant variability, with amino acid lengths ranging from 179 residues (GmGARP72 and GmGARP89) to 1226 residues (GmGARP10) (Supplementary Table S2). Correspondingly, molecular weights varied from 20,420.59 Da (GmGARP72) to 136,976.97 Da (GmGARP10), while theoretical isoelectric points (pI) ranged from 4.99 to 9.45, indicating that 76 of the family members had a pI greater than 7, while 49 had a pI below 7. Further analysis of protein stability showed that the instability indices ranged from 36.66 (GmGARP94) to 76.75 (GmGARP15), suggesting varying degrees of stability among these proteins. The Aliphatic Index, reflecting the relative volume of aliphatic side chains, ranged from 54.05 (GmGARP58) to 102.14 (GmGARP97). The Grand Average of Hydropathicity (GRAVY) values for all GmGARP family members were negative, indicating their hydrophilic nature, and subcellular localization predictions confirmed their nuclear localization, implying a role in transcriptional regulation.

3.2. Phylogeny and Classification of GmGARP Family Members

We constructed a NJ phylogenetic tree of GARP family members from Glycine max and A. thaliana, dividing the Glycine max GmGARP genes into eight groups (I, IIa~IIf, and III), consistent with the classification of A. thaliana GARP members (Figure 1). Each group was annotated in alignment with the known functions of A. thaliana GARP members. Specifically, Group I corresponds to ARR (A. thaliana Response Regulators), while Group III aligns with PHL/PHR (Phosphate Starvation Response). The subgroups IIa to IIf include APRR (A. thaliana Pseudo-Response Regulator), RVE (REVEILLE), MYBC1/LUX, HRS1/HHO (Histidine-containing Phosphotransfer Protein), KAN (KANADI), and GLK (Golden2-Like) transcription factors, reflecting the functional diversity found in A. thaliana (Figure 1). Group I contained 21 members (16.67%), subgroups IIa–IIf together included 47 members (37.30%), and Group III, the largest, comprised 58 members (46.03%). This classification reveals the broad range of regulatory roles GmGARP genes may play in Glycine max, emphasizing their importance in diverse physiological and developmental processes.

3.3. Analysis of Protein Motifs and Gene Structures of GmGARP Family Members

Ten distinct motifs were identified among the GmGARP family members (Supplementary Figure S1). Motifs 1 and 2 were classified as golden-2-like transcription factors, while Motifs 3, 4, 6, and 7 were associated with ARR domains. A NJ phylogenetic tree (Figure 2A,B) was constructed to analyze the motif composition across different groups. In Group I, GmGARP51, GmGARP80, and GmGARP110 predominantly featured Motifs 1 and 2, while the remaining 18 members displayed a more diverse composition, including Motifs 6, 7, 4, 9, 3, 2, and 1. Groups IIa to IIf consistently contained Motifs 1 and 2, with additional motifs varying across subgroups. For example, Group IIa members also had Motifs 7, 9, and 3, while Group IId exclusively contained Motif 10. In Group III, 31 proteins exhibited Motifs 8 and 5 alongside Motifs 1 and 2, while the other 27 members contained Motifs 6, 7, 4, 9, and 3 but lacked Motifs 1 and 2.

The gene structure analysis of GmGARP family members revealed considerable diversity, with exon numbers ranging from 1 to 13 and introns from 0 to 12. Within each phylogenetic group, exon and intron counts were relatively consistent (Figure 2C). In Group I, most members had either 5 exons/4 introns or 6 exons/5 introns, while one member each had 9 exons/8 introns, 4 exons/3 introns, and 8 exons/7 introns. In Group IIa, exon counts ranged from 7 to 11, with corresponding introns from 6 to 10. All three members of Group IIb had 6 exons and 5 introns. Group IIc showed more variation, with members having 1, 2, or 5 exons, and 0, 1, or 4 introns, respectively. In Group IId, members displayed 4 to 6 exons and 3 to 5 introns. Group IIe included seven members with 6 exons/5 introns and one with 3 exons/2 introns. In Group IIf, nine members had 6 exons/5 introns, while five members had 5 exons/4 introns. Group III showed the most variation, with exons ranging from 5 to 13 and introns from 4 to 12, though the majority contained 5 or 6 exons. This structural diversity highlights the functional versatility and adaptability of the GARP family in the Glycine max genome.

3.4. Localization, Replication, and Ka/Ks Values of GmGARP Family Members

The distribution of GmGARP family members across Glycine max chromosomes was systematically analyzed. These genes are located on chromosomes 1, 2, 3, 4, 5, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, and 20, with the respective counts being 8, 12, 10, 8, 7, 7, 7, 8, 9, 5, 9, 4, 9, 6, 11, and 6 genes (Supplementary Table S3, Figure 3). Notably, all GmGARP family members are situated in regions characterized by high gene density on these chromosomes. This distribution pattern indicates that there is no significant correlation between the chromosome length and the number of GmGARP family members, suggesting that gene density rather than chromosomal length might be a more critical factor influencing the localization of these genes.

The GmGARP family shows extensive genetic duplication, with 57 segmental duplication pairs and one tandem duplication pair distributed across chromosomes 1–5 and 10–20 (Figure 4A). The chromosomal locations of these pairs are detailed in Figure 4B. Glycine max shares 44 GARP homologous gene pairs with both A. thaliana and Oryza sativa, as well as 22 homologous gene pairs with Zea mays. These findings suggest that GARP genes in monocotyledonous and dicotyledonous plants may have diverged during the course of evolution (Figure 4C–E). The Ka/Ks ratio analysis, which evaluates the selection pressure on these genes, revealed that all segmental and tandem duplication gene pairs possessed Ka/Ks values of less than 1. This indicates that these genes are undergoing purifying selection, which preserves the gene function by eliminating deleterious mutations (Supplementary Table S4). However, the existence of numerous Ks values for homozygous gene pairs between Glycine max and A. thaliana, Glycine max and Oryza sativa, and Glycine max and Zea mays that could not be calculated suggests that GARP family genes may exhibit a high degree of variability in amino acid composition throughout the evolution of these species.

3.5. Analysis of Cis-Acting Elements of GmGARP Family Members

A total of 114 cis-acting elements were identified in the promoter regions of GmGARP family members, covering a wide range of regulatory functions such as light response, hormone regulation, stress response, and growth and development (Supplementary Table S5). We mapped the positions of key cis-elements in the promoter regions (Figure 5A). Hormone-related elements included the TGACG-motif, CGTCA-motif, TATC-box, AuxRR-core, ABRE, P-box, GARE-motif, TGA-element, and TGA-box. Stress response elements included MBS, LTR, TC-rich repeats, GC-motif, WUN-motif, and ARE, while growth-related elements encompassed MBSI, A-box, GCN4_motif, AACA_motif, MSA-like, HD-Zip 1, circadian, RY-element, CAT-box, and O2-site. Notably, the cis-elements showed limited conservation across different GmGARP groups. Furthermore, we quantified three categories of cis-elements: hormone regulation, stress response, and growth and development. Rare elements like HD-Zip 1, MBSI, MSA-like, and others occurred fewer than 10 times, while elements such as ARE, ABRE, CGTCA-motif, and TGACG-motif were more prevalent, with occurrences of 108, 94, 73, and 73, respectively.

3.6. Analysis of Expression Patterns of GmGARP Family Members

The expression levels (FPKM) of GmGARP family members were analyzed across nine Glycine max tissues: cotyledon, flower, leaf, pod seed, pod, root, seed, shoot meristem, and stem. Genes with FPKM values below 1 were excluded, leaving 109 GmGARP members for further analysis and heat map generation (Supplementary Table S6, Figure 6). The results revealed distinct tissue-specific expression patterns. In the cotyledon, 16 genes, including GmGARP2, GmGARP10, GmGARP15, and GmGARP122, exhibited significant expression. In the flower, 28 genes, including GmGARP7, GmGARP8, GmGARP12, GmGARP19, and GmGARP126, showed significant expression. The leaf displayed significant expression for 26 genes, including GmGARP2, GmGARP6, GmGARP10, and GmGARP120. In the pod seed, four genes, GmGARP15, GmGARP38, GmGARP107, and GmGARP112, were significantly expressed. In the pod, 21 genes, including GmGARP1, GmGARP8, GmGARP9, GmGARP12, GmGARP16, and GmGARP120, showed notable expression levels. The root had 36 genes, including GmGARP1, GmGARP3, GmGARP5, and GmGARP125, with significant expression. GmGARP126 was distinctly expressed in the seed. The shoot meristem exhibited significant expression for 27 genes, including GmGARP3, GmGARP8, GmGARP14, and GmGARP126. In the stem, 29 genes, including GmGARP4, GmGARP9, GmGARP13, GmGARP14, and GmGARP125, were significantly expressed. In conclusion, the GmGARP family members significant tissue-specific expression characteristics.

3.7. Protein Interaction Network of GmGARP Family Members

Potential interactions among GmGARP family members were predicted using the A. thaliana database via the STRING tool. The analysis identified 23 GmGARP genes involved in three distinct protein interaction networks (Supplementary Table S7, Figure 7). The first network comprises 13 members: GmGARP94, GmGARP80, GmGARP9, GmGARP110, GmGARP95, GmGARP111, GmGARP89, GmGARP114, GmGARP115, GmGARP107, GmGARP36, and GmGARP94, with GmGARP94 acting as the core gene due to its central role. The second network includes six members: GmGARP2, GmGARP15, GmGARP54, GmGARP125, GmGARP105, and GmGARP101, with GmGARP54 serving as the core gene. The third network comprises four members: GmGARP70, GmGARP98, GmGARP120, and GmGARP60, with both GmGARP120 and GmGARP60 identified as core genes. Multiple potential interactions were found within these networks, such as GmGARP95 interacting with GmGARP89, GmGARP111, and GmGARP36.

3.8. Target Gene Analysis

A total of 11,195 target genes were predicted in the Glycine max genome base on the A. thaliana ARR11 binding site, which included five matching sequences: AAGATACG, TAGATACG, AAGATTCG, AAGATATG, and TAGATTCG (Figure 8A, Supplementary Table S8). All target genes were annotated with GO and KEGG (Figure 8B,C, Supplementary Tables S9 and S10). Of these, 4577 genes were assigned to GO terms, with the most enriched in biological process enriched in four GO terms: cellular process (GO:0009987), single-organism process (GO:0044699), metabolic process (GO:0008152), and response to stimulus (GO:0050896). In terms of molecular function, most target genes were enriched in three GO terms: catalytic activity (GO:0003824), binding (GO:0005488), and nucleic acid binding transcription factor activity (GO:0001071). For the cellular component, most target genes were enriched in five GO terms: cell (GO:0005623), cell part (GO:0044464), organelle (GO:0043226), membrane (GO:0016020), and organelle part (GO:0044422). KEGG analysis annotated 2176 genes, with most enriched in the metabolism-related pathways, including glycolysis/gluconeogenesis (ko00010), the citrate cycle (TCA cycle) (ko00020), pentose and glucuronate interconversions (ko00040), and starch and sucrose metabolism (ko00500). These findings suggest that GmGARP members play a key role in regulating complex processes in Glycine max.

3.9. Expression Pattern of GmGARP Family Under Hormone Stress

To investigate the expression patterns of GmGARP genes under hormone stress, qRT-PCR was performed on 25 genes selected for their predicted cis-acting elements related to GA, ABA, and MeJA. Glycine max seedlings were treated with GA (50 μM), MeJA (200 μM), and ABA (10 μM) (Figure 9). The analysis revealed distinct and overlapping responses to these hormones. Notably, GmGARP114, GmGARP95, and GmGARP107 showed strong expression under all three treatments. Additionally, GmGARP98, GmGARP112, GmGARP115, GmGARP117, GmGARP118, GmGARP119, GmGARP121, and GmGARP122 responded significantly to two of the hormones. Eleven GmGARP genes were up-regulated following GA treatment, with GmGARP107 showing sustained expression. Peak expression for GmGARP121, GmGARP114, GmGARP118, GmGARP119, GmGARP95, GmGARP98, and GmGARP99 occurred between 2 and 16 h. In contrast, GmGARP124, GmGARP126, and others were down-regulated or unchanged after GA treatment but showed significant up-regulation under ABA, peaking between 2 and 8 h. Similarly, GmGARP122, GmGARP124, GmGARP112, GmGARP115, and GmGARP97 were not induced by GA but responded to MeJA. These findings highlight the important regulatory role of GmGARP genes in plant growth and stress responses mediated by GA, ABA, and MeJA.

3.10. Expression Pattern of GmGARP Family Under Abiotic Stress Conditions

To examine the role of the GARP family under abiotic stress, GmGARP gene expression was measured following exposure to cold (4 °C), salt (120 mM NaCl), and polyethylene glycol (PEG6000, 20%) treatments. The expression profiles under the salt and PEG6000 treatments were similar (Figure 10), with genes such as GmGARP115, GmGARP85, GmGARP87, GmGARP89, GmGARP88, GmGARP80, GmGARP81, GmGARP72, GmGARP73, and GmGARP75 significantly up-regulated. In contrast, these genes showed different or unchanged expression patterns under cold stress. After 2 h of cold exposure, most GmGARP genes were up-regulated, except for GmGARP85, GmGARP87, GmGARP88, GmGARP80, GmGARP16, GmGARP73, and GmGARP75. Peak expression for GmGARP118, GmGARP106, and others occurred between 4 and 16 h. Under PEG6000 treatment, 17 GmGARP genes were up-regulated, with peak expression after 2 h. GmGARP87, GmGARP75, and GmGARP4 continued to rise, while 14 others initially increased before decreasing. During salt stress, 26 genes showed variable patterns, with only GmGARP86 and GmGARP73 continuously up-regulated, peaking between 2 and 8 h. These results underscore the role of GmGARP genes in responding to different abiotic stresses.

3.11. Interactions Between GmGARP120 and GmGARP98

Our present findings suggest that co-regulation of GmGARP120 and GmGARP98 expression was specifically observed at 8 h after MeJA treatment. Therefore, GmGARP120 and GmGARP98 may be associated with the MeJA signaling pathway. To verify the reliability of the protein interactions analyzed for the GmGARP family genes in this study, we first expressed GmGARP120 and GmGARP98 in yeast. As shown in Figure 11a, the results of Y2H assays showed that GmGARP120-BD and GmGARP98-AD exhibited blue on SD/-Leu/-Trp/-His/-Ade/x-α-gal, confirming the interaction specificity between GmGARP120 and GmGARP98. Co-expression of GmGARP120 fused to the amino-terminal half of YFP (nYFP) and GmGARP98 fused to the carboxyterminal half (cYFP) of yellow florescent protein led to visible fluorescence in the nucleus of co-transformed Nicotiana tabacum (Figure 11b). These results further confirmed that there were specific interactions between GmGARP120 and GmGARP98.

4. Discussion

GARP transcription factors are crucial regulatory genes in plants, controlling growth, development, signal transduction, and stress responses [29]. While extensive research has been conducted on GARP family members in model plants such as A. thaliana (56 members) [12], C. sinensis (69) [27], and Brassica napus (146) [53], no comprehensive research has been conducted on the GARP gene family in Glycine max. In this study, we identified 126 GmGARP members through a genome-wide analysis of Glycine max, revealing substantial variability in their amino acid composition, molecular weight, and isoelectric point. Similar to findings in B. napus and C. sinensis, all GmGARP proteins were hydrophilic and primarily localized in the nucleus, which suggests their involvement in key regulatory pathways, particularly in transcriptional regulation. Phylogenetic analysis allowed us to categorize the GmGARP genes into eight distinct groups: Group I, Group II (subdivided into IIa, IIb, IIc, IId, IIe, and IIf), and Group III, which closely correspond to previously identified GARP families in A. thaliana and Spirodela polyrhiza [40]. This classification provides insights into the potential functional conservation of GmGARP genes across plant species. For example, the PHR1/PHL1 subgroup is known to be involved in phosphate signaling, while the HRS1/HHO subgroup is implicated in nitrate signaling [27]. Given that nutrient uptake and signaling pathways are critical for plant growth and development, understanding the roles of these GmGARP genes could offer valuable strategies for improving nutrient use efficiency in Glycine max, especially under variable environmental conditions. Furthermore, the functional diversity observed within the GmGARP family, inferred from the variable molecular characteristics of the proteins, points to their likely involvement in a broad range of physiological processes. This diversity mirrors the complexity of regulatory mechanisms needed to adapt to environmental stressors. The nuclear localization of GmGARP proteins aligns with their putative roles as transcriptional regulators, potentially interacting with other nuclear factors to modulate responses to external stimuli. For instance, genes in the KAN subfamily, which regulate leaf polarity and meristem identity, could play a significant role in optimizing Glycine max plant architecture, directly impacting yield potential.

Conserved protein motifs and gene structures not only elucidate the functional and structural characteristics of proteins, but also provide insights into the evolutionary and taxonomic relationships among family members [54,55]. In our study, we annotated 10 conserved motifs within the GmGARP members. Specifically, Motifs 1 and 2 were classified as golden-2-like transcription factors, while Motifs 3, 4, 6, and 7 corresponded to ARR domains. Our analysis demonstrated that the motif types, their distributions, and the exon–intron structures within GmGARP proteins are highly conserved among group members. However, some members displayed significant divergence in motif composition and exon/intron counts, suggesting varying degrees of evolutionary differentiation and functional diversification within the group. These findings are consistent with previous research on GARP family members across different plant species [27,40,53]. Variations in motifs may relate to the adaptive responses of these proteins to environmental pressures, emphasizing the importance of conserved motifs in plant evolution. Moreover, differences in motif distribution and exon numbers may reflect adaptive evolutionary processes in specific ecological contexts.

GmGARP family members are distributed across all 16 Glycine max chromosomes, with no significant correlation between gene count and chromosome length. Gene duplication is vital for preserving genetic integrity and stability, influencing processes such as cell division, tissue growth, genetic adaptation, and gene expression regulation [56]. In C.sinensis, 40 out of 69 CsGARP genes (57.97%) are classified as segmental duplications [27]. Similarly, S. polyrhiza has been shown to contain four pairs of segmental duplications and two pairs of tandem duplications [40]. In our study, we identified 57 pairs of segmental duplications and one pair of tandem duplications within the GmGARP family, with Ka/Ks values for these duplicated gene pairs all below 1. This indicates that GmGARP members have undergone significant purifying selection during evolution. Analyzing homologous gene pairs within the same family across different species provides insights into species differentiation and the evolutionary dynamics of gene members [57]. For instance, B.napus and A. thaliana share 152 pairs of GARP homologous genes [53], while C. sinensis and A. thaliana share 53 pairs [27]. Glycine max has 44 pairs of GARP homologous genes with both A. thaliana and Oryza sativa, and 22 pairs of GARP homologous genes with Zea mays, and the Ks values of most homologous gene pairs could not be computed to obtain the results, which suggests that the gene sequences of the GARP family members have a high degree of variability among different species. This suggests a high level of conservation in the protein sequences of GARP family members across different species. These findings highlight the evolutionary stability of GmGARP genes, underscoring their potential functional importance across diverse plant lineages. Understanding these evolutionary relationships can enhance our knowledge of the roles that GARP genes play in plant development and stress responses.

Cis-acting elements are crucial components located within the regulatory regions of genes that interact with transcription factors to modulate the transcriptional activity of adjacent genes. Analyzing these elements can provide insights into the mechanisms governing gene expression under various physiological and environmental conditions, including photosynthesis and stress responses [58]. In our study, we identified a diverse array of cis-acting elements within the GmGARP family. In addition to numerous light-responsive elements, we detected elements associated with hormone regulation, stress response, and growth and development. Notably, the number of cis-elements for six specific types—ARE, ABRE, CGTCA, TGACG, TCA-element, and TC-rich repeats—exceeded 50. This suggests that GmGARP members likely play multifaceted roles in regulating processes vital for Glycine max growth and development, contributing to the plant’s overall health [59]. Previous studies have shown that GARP family members in various species, including tomato [60], A. thaliana, B.napus, and C. sinensis, exhibit significant expression across different tissues such as leaves, flowers, and roots. Consistent with these findings, our research shows that GmGARP family members are prominently expressed in nine Glycine max tissues: cotyledon, flower, leaf, pod seed, pod, root, seed, shoot meristem, and stem. Notably, several genes, including GmGARP1, GmGARP3, and GmGARP84, exhibit high expression levels across multiple tissues, suggesting that these GmGARP members may have critical regulatory functions in diverse physiological contexts. Notably, the cis-acting elements ABRE, ARE, and TCA-element are nearly ubiquitous among the GmGARP gene members with significant expression in the shoot meristem. These cis-elements are respectively associated with anaerobic induction, ABA response, and salicylic acid (SA) response during plant growth and development. Studies have shown that anaerobic peptides and proteins synthesized under anaerobic conditions may regulate cell division and differentiation processes in plant meristem tissues and are involved in various stress responses [61]. The role of ABA in plant meristems primarily involves the regulation of cell division, elongation, and differentiation, while also contributing to stress response and adaptation, as well as stomatal opening and transpiration control [62]. Similarly, SA has been demonstrated to promote cell division and differentiation in plant meristems and to regulate overall plant growth and development [63]. These findings suggest that ABRE, ARE, and TCA-element cis-acting elements play a critical role in the regulation of shoot meristem activity in Glycine max. This comprehensive expression profile underscores the potential importance of GmGARP genes in coordinating various developmental and stress response pathways in Glycine max.

GARP transcription factors play a pivotal role in plant hormone responses and abiotic stress. For example, crossing the abi3, abi4, and abi5 mutants with the hrs1-1 mutant restored the germination phenotype of the hrs1-1 mutant, indicating that HRS1 likely functions upstream in the ABA signaling pathway during seed germination [64]. Furthermore, the glk1 glk2 double mutant displayed ABA-hypersensitivity during seed germination and seedling development, along with enhanced osmotic stress resistance [65]. In Medicago truncatula, MtHHO3 has been associated with responses to salt stress and ABA signaling [11]. Moreover, previous studies have shown that membrane ion leakage was significantly reduced in HRS1 mutants and OsNIGT1 mutants, respectively, and that HRS1 mutants and OsNIGT1 mutants continued to grow and fruit after salt stress [66]. This suggests that NIGT1/HRS1 TFs are associated with salt tolerance. In addition, GhGLK1 enhances cold and drought tolerance in Gossypium hirsutum [28]. AQUILO, a grapevine transcription factor related to GARP/MYB, enhances cold tolerance and promotes the accumulation of oligosaccharides [67]. In our study, we observed a significant up-regulation of GmGARP gene expression following treatment with GA and ABA after 2 h, and after 4 and 8 h of MeJA treatment (Figure 9). These results highlight the pivotal role of GmGARP genes in regulating plant growth and stress responses mediated by GA, ABA, and MeJA. The differential expression patterns indicate the complexity of hormonal regulation and suggest that these genes may contribute to stress adaptation mechanisms in Glycine max. Additionally, the promoter regions of GARP genes contain multiple cis-acting elements related to abiotic stresses, including MBS, LTR, TC-rich repeats, GC-motif, WUN-motif, and ARE. While research on the relationship between GARP transcription factors and abiotic stresses is limited, the grape transcription factor VaAQ, homologous to HRS1, HHO2, and HHO3, has been shown to enhance low-temperature tolerance in transgenic A. thaliana and Vitis amurensis. VaAQ regulates the expression of CBFs or DREBs, which are key genes in low-temperature signaling through the ABA pathway, thereby improving tolerance to low-temperature stress [67]. Our findings also indicate that GmGARP gene expression patterns in Glycine max are similar under salt and drought stress but differ under cold stress (Figure 10). For instance, the role of GmGARP86 in responding to drought stress has also been established [31]. This underscores the significant role of GmGARP genes in adapting to various abiotic stresses, revealing distinct regulatory mechanisms.

Transcription factors often regulate diverse biological processes through protein–protein interaction. In our study, we predicted that 23 GmGARP members would engage in potential protein interactions, organized into three distinct network modules. To validate these predictions, we analyzed the interaction potential of specific protein pairs and confirmed molecular interactions between GmGARP120 and GmGARP98 using both yeast two-hybrid and BiFC assays (Figure 11). Our findings indicate that the co-regulation of GmGARP120 and GmGARP98 expression is specifically observed 8 h post-MeJA treatment, while their expression levels do not align consistently under ABA or GA treatments. Therefore, the GmGARP120–GmGARP98 interaction may be associated with the MeJA signaling pathway, but additional molecular experiments are required to confirm this hypothesis. This finding sheds light on the biological functions of GmGARP in Glycine max. While the precise roles of these GmGARP remain to be elucidated, our results lay a foundation for future investigations into their involvement in auxin signal transduction and abiotic stress regulatory mechanisms in Glycine max. Additionally, transcription factors regulate gene expression by binding to specific cis-regulatory sequences in the promoters of target genes. However, research on GARP transcription factor target genes in plants is limited. In this study, we identified 11,195 potential target genes within the Glycine max genome, utilizing the transcription factor binding map of A. thaliana ARR11 as a reference. We also predicted the associated GO functions, KEGG pathways, and domains of these target genes, offering valuable insights for further exploration of GARP gene interactions in Glycine max. This comprehensive approach not only enhances our understanding of GmGARP’s roles in gene regulation, but also provides a framework for investigating their contributions to Glycine max development and stress responses.

5. Conclusions

In this study, we identified 126 GmGARP genes from the Glycine max genome and conducted a comprehensive bioinformatics analysis of the GmGARP family members. The phylogenetic analysis classified the GmGARP family members into eight groups, with highly conserved structures and protein motifs within the same group. Gene duplication events appear to be a major factor in the expansion of the GmGARP family. Expression pattern analysis revealed significant tissue-specific expression characteristics of GmGARP members, which also showed strong responses to exogenous hormones and abiotic stresses. Specifically, proteins GmGARP120 and GmGARP98 were confirmed to interact through bioinformatic predictions, yeast two-hybrid systems, and BiFC assays. These findings provide valuable insights into the characteristics and functions of the GmGARP family members.