Soybean material was grown under conditions as previously described (Sun et al., 2019). Full and consistent growing Williams 82 seeds were sterilized using 70% alcohol for 1 minute and then rinsed with distilled water. Sterilized soybean seeds were germinated on moist filter paper. Four well germinated seeds were then selected and sown in pots containing vermiculite. All seedlings were grown under a photoperiod of 16 h/8 hours (light/dark) and a temperature of 25 °C/20 °C (light/dark) and were regularly watered with Hoagland’s liquid medium. For 2-week-old seedlings, 200 mL of 200 mM NaCl solution and 10% PEG6000 solution were added to simulate salt and drought stress treatments, respectively. Based on previous studies (Chen et al., 2014; Lv et al., 2022; Zhang et al., 2025), leaf tissue was collected at 0, 6, 12, and 24 hours post-treatment. Three biological replicates were prepared at each time point, with each replicate comprising at least three seedling pots. The 0 hour sample served as the untreated control group. All samples were immediately placed in liquid nitrogen upon collection and stored at -80 °C for subsequent RT-qPCR analysis of gene expression patterns.

For hormone treatments refer to previous descriptions (Xiong et al., 2024). In summary, this study treated 2-week-old soybean seedlings with five hormones at a concentration of 100 μmol/L (abscisic acid: ABA, indole-3-acetic acid: IAA, gibberellin: GA, methyl jasmonate: MeJA, and salicylic acid: SA), with samples collected at 0, 6, 12, and 24 hours post-treatment. Samples collected at 0 hour served as untreated controls. Changes in gene expression patterns following hormone treatment were analyzed at 6, 12, and 24 hours. Sampled leaves were rapidly placed in liquid nitrogen and stored at - 80 °C for subsequent RT-qPCR. There were at least three pots of seedlings per treatment time point and at least three biological replicates per treatment time point.

Total RNA from soybean leaves was extracted using the SteadyPure Plant RNA Extraction Kit (AG21019, Accurate Biology, China) according to the manufacturer’s instructions. First-strand cDNA was synthesized using UEIris RT mix with DNase (All-in-One) (R2020S, UElandy, China). To analyze the expression patterns of all genes, RT-qPCR analysis was performed in a CFX ConnectTM real-time system (BIO-RAD, USA) using Universal SYBR Green Super Mix (S2024S, UElandy, China). Data were normalized using Actin11 as an internal control and relative gene expression was calculated using the 2^-ΔΔCT method. The experiment was conducted using three independent biological replicates, with three technical replicates for each sample. Gene-specific primers for RT-qPCR are listed in Supplementary Table 1.

Isolated histones were extracted using EpiQuik Total Histone Extraction Kit (Epigentek, Farmingdale, NY, USA, cod. OP-0006) according to the instructions. Briefly, leaves of soybean seedlings treated with salt and drought stress for 12 hours were ground into powder, resuspended with 1 mL of diluted prelysis buffer (1×), and incubated with ice bath agitation for 10 minutes to remove the plasma membrane. after centrifugation at 4°C, 9391× g for 1 minutes, the precipitate was resuspended in 50 μL of lysis buffer and incubated for 30 minutes on ice. The lysate was centrifuged at 4°C, 13523× g for 5 minutes and the supernatant was transferred to a new tube. To prevent protein aggregates from interfering with histone extraction, 0.3 volume of DTT-free equilibration buffer was immediately added to each sample for western blotting analysis.

Western blotting was performed as previously described (Liu et al., 2023). Proteins were separated on a 12% sodium dodecyl sulfate‐polyacrylamide gel and transferred onto a polyvinylidene difluoride membrane (Millipore, Germany). The following primary antibodies were used: anti‐H3 (ab1791; Abcam, UK), anti‐H3ac (06‐599; Millipore, Corp, USA), anti-H3K9ac (ab10812; Abcam, UK), anti-H3K18ac (ab1191; Abcam, UK), anti‐H4ac (ab177790; Abcam, UK), anti-H4K5ac (ab51997; Abcam, UK), anti-H4K8ac (ab45166; Abcam, UK). Based on previous research findings, the above acetylated antibodies are expressed in plants in response to salt stress or drought stress (Liu et al., 2023; Xue et al., 2024; Yung et al., 2022). Histone H3 was used as a loading control. Signals were detected using the SuperSignal West Pico Plus chemiluminescent substrate (Thermo Scientific; Vazyme, China).

To identify members of the soybean GmHAT and GmHDAC families, protein sequences of the A.thaliana and O.sativa HATs and HDACs genes were downloaded from the TAIR database (https://www.arabidopsis.org) and the Rice Genome Annotation Project (https://rice.uga.edu/), respectively. The protein sequences of HATs and HDACs from Arabidopsis and rice, respectively, were used as input sequences to search the entire soybean genome (Wm82.a6.v1) through the BLASTP (https://phytozome-next.jgi.doe.gov/blast-search) program in the Phytozome 13 database using the default parameters. Next, based on the previous research methodology (Du et al., 2022; Guo et al., 2024), the following entries were downloaded from the Pfam database: (GNAT: PF00583, CBP: PF08214, TAFII250: PF09247, MYST: PF01853; RPD3/HDA1: PF00850, SIR2: PF02146, HD2: PF17800). We employed the same strategy to search for conserved domains in putative GmHATs and GmHDACs, ensuring they belong to the HAT and HDAC families. Finally, verification was performed using the CD-Search database (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) and the SMART tool (https://smart.embl.de/smart/change_mode.cgi).

Chromosomal positions and amino acid numbers of the genes were obtained using gff annotation files obtained from Phytozome 13. The basic physicochemical properties of GmHATs and GmHDACs, including theoretical isoelectric point (pI), molecular weight and instability coefficient, were analyzed using ExPASy (https://www.expasy.org/). The online software Plant-mPLoc (http://www.csbio.sjtu.edu.cn/bioinf/plant-multi/) was used to predict the subcellular localization of GmHATs and GmHDACs.

To classify the GmHATs and GmHDAC, protein sequence comparison was performed using mafft (Katoh et al., 2005) based on the pre-downloaded protein sequences of O.sativa, A.thaliana and G.max. The evolutionary tree is generated using fasttree (Price et al., 2009) and finally visualized using the R package ggtree (Yu et al., 2017). Gene structure mapping based on gff annotation files using the R package transPlotR. The conserved motifs of the genes were analyzed using the online database MEME (http://meme-suite.org/tools/meme) with the maximum conserved motif search value set to 20 and other parameters defaulted. The motif distribution maps were generated by TBtools software using an XML file containing motif pattern information obtained from MEME software. Conserved structural domains of genes were visualized using TBtools (Chen et al., 2020a) based on CD-Search database (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) results.

The GmHATs and GmHDACs were localized to chromosomes using TBtools based on G.max genome and annotation file data. Tandem and fragment duplication events of GmHATs and GmHDACs were analyzed using Multiple Covariance Scanning Kit X (MCScanX) (Zhu et al., 2014), using default parameters. The multiple synteny plot tool in Circos tool and TBtools was used to visualize the covariance between members of the G.max GmHATs and GmHDAC families and other species (dicots: A. thaliana and C. arietinum; monocots: Z. mays and O. sativa).

Sequences 2000 bp upstream of the transcription start site were extracted as promoter sequences. Sequences were submitted to the online database PlantCARE (Lescot et al., 2002) (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) to analyze potential cis-elements in the promoter region. The analysis results were visualized using TBtools. Protein interaction networks for GmHATs and GmHDACs were predicted using the STRING (https://cn.string-db.org/) database. The sequences of all proteins were used as a query, with G.max selected for species and other parameters defaulted. Subsequently, protein interaction networks were exported and visualized using Cytoscape (version 3.9.1).

Download raw RNA-seq data from NCBI for different tissues of Williams 82 varieties that have been publicly reported (stem: PRJDB7219; leaf, root: PRJNA407016; seedling: PRJNA525277; flower, pod: PRJNA236472; Seed: PRJNA677883; Cotyledon: PRJNA262564). The data were aligned to the G.max genome (Wm82.a6.v1) using hisat2 software, and the transcripts per kilobase million (TPM) values of the genes were further calculated using stringtie and visualized using TBtools.

We used A.thaliana and O.sativa HAT and HDAC protein sequences on the basis of the Phytozome 13 BLAST program as input sequences to search in the G.max (Wm82.a6.v1) genome. After identification of the conserved structural domains of HAT and HDAC, a total of 12 GmHATs and 28 GmHDACs were identified. We named the subfamilies according to their position on the chromosome, and analyzed their protein molecular weights, theoretical isoelectric points (pI), subcellular localization, and other physicochemical properties (Supplementary Table 2).

Of the 12 GmHATs, the amino acid numbers ranged from 434 aa (GmHAM1) to 1890 aa (GmHAF2), and the molecular weights ranged from 50.10 kDa (GmHAM2) to 214.56 kDa (GmHAF2). The smallest pI was 5.22 (GmHAG3) and the largest was 8.45 (GmHAC1). Interestingly, 12 GmHAT proteins were classified as unstable, with instability coefficients all exceeding 40, necessitating careful handling in vitro. Concerning 28 GmHDAC proteins, the amino acid numbers ranged from 114 aa (GmHDA10) to 656 aa (GmHDA5) and the molecular weights ranged from 13.01 kDa (GmHDA10) to 73.05 kDa (GmHDA5). The smallest pI was 4.61 (GmHDT6) and the largest was 9.6 (GmHDA10), with the majority (22 of 28) having a pI < 7, suggesting that the GmHDAC gene family favors acidic amino acids. The instability coefficients ranged from 24.39 (GmHDA14) to 52.1 (GmHDT1). Predictive analysis of subcellular localization showed that soybean GmHAT and GmHDAC exhibited diverse distribution patterns. Among them, 10 GmHAT proteins were mainly localized in the nucleus; 3 SIR2 subfamily members were specifically localized in chloroplasts; whereas GmSRT3 showed chloroplast and cytoplasmic co-localization patterns. Notably, 2 GmHAT and 8 GmHDAC belong to the nucleoplasmic co-localization, suggesting that 10 genes may catalyze both histone acetylation and non-histone acetylation. Taken together, the soybean GmHAT and GmHDAC family members showed significant diversity in physicochemical properties and subcellular localization, suggesting that they may be functionally differentiated in epigenetic regulation and participate in the dynamic balance of histone and non-histone acetylation/deacetylation through organelle-specific localization.

Based on the protein sequences and classification of A.thaliana and O.sativa HATs and HDACs, we performed phylogenetic analyses of soybean GmHATs and GmHDACs using the neighbor-joining method (Figure 1). The results indicated that the 12 GmHAT members were categorized into four subfamilies, with the GNAT and CBP subfamilies containing 4 members each, and the TAFII250 and MYST subfamilies 2 members each. The 28 GmHDACs were mainly categorized into 3 subfamilies, with 18 members of the RPD3/HDA1 subfamily, 4 members of the SIR2 subfamily, and 6 members of the HD2 subfamily. Among the 18 RPD3/HDA1 subfamily members, Class I had 10 members, Class II had 6 members, and Class III and Class IV each had 1 member.

Analysis of the phylogenetic tree for the GmHAT and GmHDAC families showed results consistent with Figure 1 (Figure 2). To analyze the structural composition of the GmHATs and GmHDACs, we constructed a structural map based on the soybean genome sequence and annotation files, including untranslated regions (UTRs), exons and introns (Figure 2). Results showed that all genes contained UTRs, with great variation in the number of exons and introns. The number of GmHATs exons ranged from 9 to 21 and the number of introns ranged from 9 to 20. While GmHDACs exons are between 3 and 17 and introns are between 2 and 17. We identified GmHAT and GmHDAC proteins through an online server (MEME). A total of 20 motifs were identified in 12 GmHATs and showed great variation (Figure 2). For GmHATs, harboring motifs ranged from 4 to 20. For these GmHATs, a total of 13 conserved structural domains were identified, and the number of GmHATs structural domains ranged from 1 to 6. The GNAT subfamily genes contain the Hat1_N and Bromo_gcn5_like structural domains; the CBP subfamily genes all contain the HAT_KAT11_superfamily structural domain; the MYST subfamily genes all contain the PLN00104 structural domain; and the TAFII250 subfamily all contain the Bromodomain. The motifs of the 28 GmHDACs proteins varied significantly, and a total of 20 were identified (Figure 2). The motifs of GmHDACs proteins range from 1 to 11. Only one motif was present in GmHDA10, whereas both GmSIR2 subfamilies contained motif 17 and motif 20. Conserved structural domain analysis revealed that both RPD3/HDA1 contain HDAC structural domains; NPL structural domains are present only in the plant-specific GmHD2 subfamily; and SIRT4 and SIRT7 are present only in the GmSIR2 subfamily. This result is highly consistent with the phylogenetic analysis.

Using the soybean genome and annotation files, we analyzed the positions of 12 GmHATs and 28 GmHDACs on the chromosomes, with all 40 members displayed in specific positions (Figure 3; Supplementary Figure 1). The 40 genes were unevenly distributed across the eighteen chromosomes (except for Gm16 and Gm20), and genes of the same subfamily were randomly distributed across the chromosomes. Gm01, Gm02, Gm09, Gm10, Gm13, Gm14, Gm15, and Gm18 contain 1 gene each; Gm03 and Gm08 contain 2 genes; Gm11 and Gm19 contain 3 genes each; Gm05, Gm06, Gm12, and Gm17 contain 4 genes each; Gm04 contains 5 genes. Above results suggest that GmHATs and GmHDACs may play important roles in the complex life processes of soybean.

Tandem duplication is a major driver of gene family expansion in plants, and chromosomal regions in the 200-kb range containing two or more genes are defined as tandem duplication events (Chen et al., 2017). A total of 2 tandem repeat events were observed in this study, including GmHDA12 and GmHDT2, and GmHDA13 and GmHDT3 (Supplementary Figure 1). Intra-species colinearity analysis revealed 11 segmental repeat events for GmHAC1, GmHAC2, GmHAC3, GmHAG1, GmHAG2, GmHAG3, GmHAG4, GmHAF1, GmHAF2, GmHAM1, and GmHAM2. Among them, GmHAC1, GmHAC2 and GmHAC3 showed correlation (Figure 3A). GmHDACs exhibited 25 fragment repeat events. Among them, GmHDA3, GmHDA9, GmHDA6, and GmHDA16 showed linear correlation; GmHDT1, GmHDT4, GmHDT5, and GmHDT6 showed linear correlation; GmHDT2, GmHDT4, GmHDT5, and GmHDT6 showed linear correlation; GmHDT3, GmHDT4, and GmHDT5, GmHDT6 showed linear correlation (Figure 3B). To further understand gene duplication events between GmHATs and GmHDACs species, we compared the homology of soybean with A.thaliana and C.aruetunum (both dicots) (Figure 3C), and O.sativa and Z.mays (both monocots) (Figure 3D). The results showed a higher degree of colinearity between the soybean and dicots genomes compared to the monocots plants. Additionally, we observed that most A.thaliana and C.aruetunum HATs and HDACs have corresponding direct homologs in soybean, with many genes having more than two direct homologs. The presence of orthologous gene pairs indicates conservation of gene function and the potential for similar biological processes between these species.

To investigate the potential regulatory mechanisms of GmHATs and GmHDACs during soybean development, we extracted the promoter sequence 2.0 kb upstream of the transcription start site and analyzed its cis-element using the PlantCARE website (Figure 4). The results show that common cis-element in promoter and enhancer regions such as TATA-box and CAAT-box are still the most widely distributed. Then, we selected 26 cis-element related to plant development, phytohormone response and abiotic stress. Light responsive elements and Low-temperature responsiveness (LTR) are present in all GmHATs and GmHDACs (Figure 4; Supplementary Figure 2). GmHATs all contain Abscisic acid responsiveness (ABRE). For GmHDACs, all but GmHDA8 and GmHDA15 contained more than one Ethylene response element (ERE). In addition, results indicated that the promoters of GmHATs and GmHDACs contain multiple stress-responsive elements, such as Drought-inducibility (MBS), Low-temperature responsiveness (LTR), Wound-responsive element (WRE3), Anaerobic induction (ARE), and Defense and stress responsiveness (STRE) (Figure 4; Supplementary Figure 2). Taken together, GmHATs and GmHDACs may play positive or negative regulatory roles in different abiotic stresses.

?>Interaction network analysis provides an effective way to analyze the biological functions of proteins and reveal their molecular mechanisms. Therefore, we predicted the interacting proteins of GmHAT and GmHDAC using the STRING database. A total of 17 nodes were identified in GmHATs, and 61 interacting proteins were observed in these nodes (Supplementary Figure 3A; Supplementary Table 3). Whereas a total of 30 nodes were identified in GmHDACs, 224 interacting proteins were observed in these nodes (Supplementary Figure 3B; Supplementary Table 3). Non-interaction of CBP and TAFII250 subfamily members in the GmHAT family. However, GmHATs interacted with C2 domain-containing proteins (I1N3I5, A0A0R0J7I7), Homeobox domain-containing proteins (I1LMP3, I1JI62) and Histone deacetylase domain-containing protein (I1MXC3). For GmHDACs, the HD2 subfamily (except for GmHDT4) did not interact with either of the other two subfamilies. Possible interaction between SIR2 and RPD3/HDA1 subfamily members and possible interaction with pathogenesis-related homeodomain protein (A0A368UH09, K7LC20, K7MK45, I1LVM1), HhH-GPD domain-containing protein (K7LYE7) and HhH-GPD domain-containing protein.

To further understand the functions of GmHAT and GmHDAC, we analyzed their expression levels in seedling, leaves, cotyledon, stem, pod, root, flower, and seed (Figure 5). The results showed that GmHATs and GmHDACs existed widely expressed and the expression levels were highly variable in different tissues. Among the 12 GmHATs and 28 GmHDACs, expression levels of 38 genes were detected in at least one tissue. In contrast, expression levels of GmHAC4 and GmHDA10 were not detected in any of the tissues, and these genes may be pseudogenes or have unique expression patterns not detected in our study. Some genes (HAC1/2/3, HAF1/2) were not differentially expressed in the tissues detected. However, many genes exhibit tissue-specific expression patterns. For instance, GmHAG2, GmHAG4, and GmHAM2 had the highest expression in seeds, and GmHDT1 had the highest expression in stem. GmHDT5 had the highest expression in flower. The differential gene expression levels among different tissues provide an important basis for screening candidate genes with more potential for in vivo studies.

The promoter regions of GmHAT and GmHDAC contain a large number of stress-related homeostatic regulatory elements, suggesting that they may play important roles in soybean adversity stress response (Figure 4). In this study, we used salt and drought stress as examples and randomly selected 16 genes out of 28 for analysis (Figure 6). The results indicated that in GmHATs, both GmHAG4 and GmHAM2 were significantly induced to be expressed by both salt and drought stresses (Figure 6A). GmHAG3 was not significantly induced by salt stress, but was significantly induced by drought treatment at 12 hours and 24 hours (Figure 6A). For GmHDACs, both salt and drought significantly induced the expression of the RPD3/HDA1 subfamily gene GmHDA16. At 12 hours and 24 hours post-stress treatment, GmHDA17 expression was also significantly induced. Conversely, GmHDA15 expression was significantly suppressed at 12 hours and 24 hours under salt treatment and at 12 hours under drought treatment (Figure 6B). Compared with other genes, the overall expression levels of SIR2 subfamily genes were low (Figure 6C). GmSRT1 expression was significantly induced by salt stress, while drought treatment alone induced it significantly after 12 hours. GmSRT2 showed opposite responses to salt and drought treatments at 6 hours: salt significantly induced its expression, while drought significantly suppressed it. GmSRT3 and GmSRT4 were significantly induced at 6 hours and 12 hours under salt and drought treatments. Meanwhile, GmSRT3 was significantly induced at 12 hours and inhibited at 24 hours under salt stress, while GmSRT4 showed no significant changes in response to salt stress (Figure 6C). The expression of other HD2 genes, except for GmHDT6, was induced up-regulated to varying degrees with treatment time (Figure 6D). Notably, the HD2 subfamily gene, GmHDT2, was significantly induced to be expressed by 24 hours of salt stress treatment, which was 61-fold higher than that of normal conditions, whereas drought stress treatment for 24 h up-regulated its expression by 43-fold.

The expression profiles of GmHAT and GmHDAC were significantly changed under salt and drought stress, so we further analyzed the histone acetylation patterns of soybean seedlings in response to salt and drought stress. We analyzed histones extracted from soybean seedlings treated with salt and drought stress for 12 hours by immunoblotting. The results showed that both salt and drought treatments resulted in a significant reduction in histone H3 and H4 acetylation levels (Figures 7A, B). Among these, H3K9ac showed no significant changes before and after stress treatments, but H3K18ac decreased markedly following salt and drought treatments. H4 acetylation analysis revealed no change in H4K5ac, but H4K8ac exhibited significant alterations following both stress treatments (Figures 7A, B). Interestingly, reduced acetylation showed a negative correlation with the high expression of HD2 family genes, particularly GmHDT2, GmHDT3, GmHDT4, and GmHDT5. Previous reports indicated that Arabidopsis RNAi-mediated HD2 family gene HDT1 (AtHD2A) leads to loss of H3K9 deacetylation (Lawrence et al., 2004). Overexpression of the HD2 family gene OsHDT701 in rice resulted in reduced H4K5/K16 acetylation levels (Ding et al., 2012). During drought stress response, HD2A and HD2B negatively regulate ABA biosynthesis by deacetylating H4K5ac on the key ABA synthesis gene NCED9, thereby suppressing its expression (Han et al., 2023). We therefore hypothesize that reduced acetylation levels, particularly at the H4 site, may be regulated by the aforementioned HD2 family genes, though this requires further validation.

Plant hormones act as signaling molecules, usually at very low concentrations, to regulate plant physiological processes through a complex network of signal transduction (Li et al., 2019a). This delicate coordination mechanism enables plants to flexibly adapt to internal developmental demands and external environmental changes. Cis-element analysis of the promoter regions of GmHAT and GmHDAC revealed the presence of ABA, GA, IAA, MeJA, and SA hormone-responsive elements, suggesting that they may be jointly involved in the regulation of soybean growth and development and stress response processes (Figure 4). To further investigate the roles of GmHAT and GmHDAC in the hormone response mechanism, we treated soybean seedlings with five hormones separately and sampled them at different treatment time periods to detect the expression patterns of the genes in Figure 6 (Figure 8). There were significant differences in the relative expression of all genes at different plant hormone treatment time periods, indicating that all of these genes responded to the five hormone treatments listed above. Although the expression pattern of each gene is not identical, some commonalities and individuality can still be observed. ABA, GA, and IAA treatments resulted in the down-regulation of the expression of most of the genes first and then sustained up-regulation (Figures 8A–C). All 16 genes showed significant up-regulation of expression after 24 hours of MeJA treatment (Figure 8D). GmHATs, RPD3/HDA1, and the two SIR2 subfamily genes showed a consistent trend of up-regulation at all times after SA treatment, but the expression of GmSRT3 and GmSRT4 was unchanged at all times (Figure 8E). In addition, the expression levels of GmHAG4 were consistently up-regulated after all four hormone treatments except ABA treatment (Figure 8). Notably, all six HD2 subfamily genes showed significant up-/down-regulation of their expression in the five hormone treatments, suggesting that they may be of particular importance in plant hormone response (Figure 8). Furthermore, GmHATs and GmHDACs play important roles in soybean in response to different plant hormones.