The nucleic acid and protein databases of Arabidopsis, rice, soybean, apple, and grape were downloaded from the Ensemble Plants database.¹ The hidden Markov model (HMM) of TIFY domain (PF06200) was obtained from Pfam.² We then used the hmm-search program HMMER3.1 (Prince and Pickett, 2002; Xia et al., 2017) to identify the TIFY HMM for the TIFY proteins in the resulting protein databases. The 18 Arabidopsis AtTIFY protein sequences were obtained from TAIR³ and used to search the TIFY proteins from rice, soybean, apple, and grape protein databases by the BLASTp program of basic local alignment search tool (BLAST; Wang et al., 2017). We compared the results of the two methods to confirm TIFY candidate genes in these species. These candidate genes were identified in their domains with SMART⁴ and CDD⁵ to ensure that the TIFY domain was in sequence (Letunic et al., 2002; Marchler-Bauer et al., 2002). Finally, the ExPASy⁶ ProtParam tool was used to query the physical and chemical properties of the GmTIFYs (Appel et al., 1994; Wang et al., 2020a).

The TIFY protein sequences of Arabidopsis, rice, soybean, apple, and grape were compared using ClustalW in the MEGA-X software. The maximum likelihood (ML) method was used to construct a phylogenetic tree for analyzing the phylogenetic relationship between TIFYs (Kumar et al., 2016; Leng et al., 2021). The bootstrap method was used with 1,000 replicates. The methods and parameters were the same as the Jones–Taylor–Thornton (JTT) model, gamma-distributed rates (G), and the gamma parameter 1.

The position information of the GmTIFY family genes was extracted from the GFF3 file of the soybean genome. The location and distribution of GmTIFY family genes were visualized on the chromosomes using Map Gene 2 Chromosomal (Jiangtao et al., 2015; Wang et al., 2020a).⁷

For gene duplication analysis, the TBtools software was used to identify the duplication events of the soybean genome and GmTIFY genes. The collinearity pairs of GmTIFY genes were extracted and used to visualize a synteny map with the CIRCOS software (Lestari et al., 2013).

The TIFY coding sequences were aligned using ClustalW software. The alignment results were converted to PAML format using EasyCodeML and a tree file in Newick format was built using MEGA-X. The selection pressure was estimated using the branch model of EasyCodeML. The ratio of non-synonymous to synonymous substitution rates (ω) was determined by the free-ratio model and the two-ratio model among the branches of the TIFY tree file (Gao et al., 2019).

The motif information of the soybean GmTIFY proteins was analyzed using the MEME online tool (Bailey et al., 2009; Wang et al., 2019).⁸ The resulting files and soybean gene structure annotation files were imported into TBtools for visualization.

The 2,000 bp promoter sequences were submitted to the PlantCARE website⁹ to analyze the cis-acting elements of its family members (Guo et al., 2007; Su et al., 2020; Wang et al., 2020a). The resulting file was imported into GSDS¹⁰ for visualization.

The RNA-seq data of GmTIFY family members in different tissues and organs were downloaded from the Phytozome database.¹¹ The transcriptome data of several abiotic stresses were obtained from previous studies (NCBI SRA accession: PRJNA694374; Shi et al., 2018). TBtools software was used to visualize the expression levels of GmTIFYs.

The soybean variety Zhonghuang39 was used for this study. The soybeans were planted in a greenhouse in a mixture of humus and vermiculite (humus: vermiculite = 1:1). Seven-day-old soybean seedlings were treated with 10% PEG6000 and 250 mM NaCl, respectively. The samples were collected at 0, 0.5, 1, 2, 4, 8, 12, and 24 h after treatments (Li et al., 2017; Zhang et al., 2019).

An RNA plant extraction kit (Zhuangmeng, Beijing, China) was used to extract total RNA from soybean leaves and TranscriptR All-in-One First-Strand cDNA Synthesis SuperMix (TransGen Biotech, Beijing, China) was used for reverse transcription. The primers designed by Primer Premier 5.0 software were listed in Supplementary Table 1. The eukaryotic elongation factor 1-β (GmELF1b) was used as the internal control (Jian et al., 2008). An Applied Biosystems 7500 Real-Time PCR System was used to perform RT-qPCR. The 2^−ΔΔCT method was used to analyze the quantitative results analysis (Udvardi, 2008PC). Each experiment was performed with three biological replicates.

Pectinase and cellulase were used to lyse fresh Arabidopsis leaves and obtain Arabidopsis protoplasts. The gene coding regions were cloned into the 16318hGFP expression vector. The fusion expression vector GmTIFY10e-hGFP and GmTIFY10g-hGFP were transformed into Arabidopsis protoplasts mediated by PEG4000, respectively (He et al., 2016). After 18 h of incubation at 22°C in the dark, the GFP fluorescence signal was detected using a laser confocal microscope (Zeiss LSM 700, Germany; Riechmann et al., 2000).

Arabidopsis (Col-0) seeds were sterilized with 75% alcohol for 15 min. Sterilized Arabidopsis seeds were sprinkled on ½ MS medium and maintained at 4°C for 4 days, after which they were moved to a growth incubator at 22°C under a 16 h light and 8 h dark cycle. When the seedlings grew to four leaves, they were transferred to a mixture of humus and vermiculite for subsequent experiments (Riechmann et al., 2000; Du et al., 2018).

The coding regions of the GmTIFY genes were subcloned into the pCAMBIA1302 vector. The constructed pCAMBIA1302-GmTIFY10e and pCAMBIA1302-GmTIFY10g were transformed into Arabidopsis using the floral dip method, respectively (Clough and Bent, 1998). Positive lines were selected on ½ MS medium plates containing hygromycin (35 mg/L) and were further verified using PCR. The same method was used until transgenic three generation (T₃). The expression levels of transgenic lines were determined by RT-qPCR and three homozygous T₃ lines with the highest expression levels were used for the subsequent phenotypic analysis (Li et al., 2017).

For the experiment of root growth, 5-day-old seedings were transferred to MS medium and MS medium with 125 mM NaCl for another 7 days, after which the lengths of primary root and fresh weights were measured (Wang et al., 2019). For salt treatment, 5-day-old seedlings were transferred to the soil, and then, 21-day-old seedlings were treated with 250 mM NaCl for 14 days (Wang et al., 2019). All experiments contained three independent replicates.

To obtain the overexpression vector of GmTIFY genes, the coding regions of GmTIFY genes were ligated with the pCAMBIA3301 vector to obtain recombinant plasmids (Kereszt et al., 2007; Zhao et al., 2017). To obtain the RNA interference expression vector, a 546 bp interference fragment consisting of a 200 bp target fragment and its antisense sequence connected by 146 bp zeol dehydrogenase gene sequence was synthesized and inserted into pCAMBIA3301 (Wang et al., 2019).

The constructed overexpression vector, interference expression vector, and empty pCAMBIA3301 vector were transferred to the Bacillus rhizobacillus (B. rhizobacillus) strain K599, and the recombinant vector was transferred to the hypocotyl of soybean via the A. rhizogenes-mediated method (Wang et al., 2015). The injected plants were cultured under high humidity conditions in a greenhouse until hairy roots grew at the infected site. When the hairy root reached about 5 cm long, the hypocotyl was removed below 0.5–1 cm of the infection site. At the same time, the seedlings were transplanted in mixed soil and cultured in the greenhouse for 7 days (Wang et al., 2020a). The positive soybean plants were subjected to salt stress test.

The leaves and roots of plants were used to determine physiological indicators under salt stress. A Physiological Index Test Kit (Cominbio, Suzhou, China) was used to test the contents of malondialdehyde (MDA), proline (PRO), catalase (CAT), and peroxidase (POD) in leaves and roots (Shi et al., 2018). The whole leaves and roots were soaked in nitroblue tetrazolium (NBT) for overnight staining. After staining, the samples were soaked in a decolorizing solution (30% glycerol and 70% ethanol) to decolorize the sample until it turned white (Du et al., 2018). All experiments were performed in three biological replicates.

The co-expression genes were obtained from the Phytozome database. The enrichment analysis was performed by Database for Annotation, Visualization and Integrated Discovery (DAVID) online tools.¹² R software was used to visualize the results of the enrichment analysis.

One-way ANOVA test analysis was performed in Microsoft Excel 2007. Data were shown as means ± standard deviation (SD), with a p-value cutoff of 0.05 and 0.01. The method was used to analyze the RT-qPCR results and physiological indicators.

The BLASTp program and hmmsearch program were used to search for TIFY genes in databases of Arabidopsis, rice, soybean, apple, and grape. We then compared the results of two programs and identified the TIFY domain using SMART and CDD to confirm TIFY members in five species. Finally, 38 GmTIFY genes were identified in the soybean genome. The number of TIFY genes in Arabidopsis (18), rice (20), grape (19), and apple (30) are consistent with previous reports (Ye et al., 2009; Bai et al., 2011; Li et al., 2014). The GmTIFY genes were named according to their relationship with Arabidopsis and their location on the chromosomes (Table 1). The protein lengths, molecular weights (MV), and isoelectric points (pI) are provided in Table 1. In 38 GmTIFY genes, the coding sequences range from 133 (GmTIFY5d) to 415 aa (GmTIFY8a); the MW ranges from 15.20 (GmTIFY5d) to 43.69 kD (GmTIFY8a) and the pI ranges from 4.63 to 9.89.

To investigate the phylogenetic relationships of GmTIFYs, we constructed a phylogenetic tree using TIFY genes from Arabidopsis, rice, soybean, grape, and apple (Figure 1). The TIFY genes were divided into four subfamilies of TIFY, JAZ, ZIM, and PPD according to specific domains (Zhang et al., 2012). The JAZ subfamily contains the largest number of TIFY genes including 12 AtTIFY genes, 11 VvTIFY genes, 22 GmTIFY genes, 22 MdTIFY genes, and 16 OsTIFY genes. The ZIM subfamily contains three AtTIFY genes, five VvTIFY genes, nine GmTIFY genes, two MdTIFY genes, and two OsTIFY genes. The PPD and TIFY subfamilies have no genes in monocotyledon plants which is consistent with the previous research results (Zhang et al., 2012; Li et al., 2014). There are two AtTIFY genes, two VvTIFY genes, two GmTIFY genes, and two MdTIFY genes in PPD subfamily and one AtTIFY genes, one VvTIFY genes, five GmTIFY genes, and four MdTIFY genes in TIFY subfamily. Further analysis revealed that the number of TIFY genes in soybeans and apples was almost twice than that of Arabidopsis and grapes. This may result from more events of chromosomes doubling, fusion and rearrangement occurred in soybean and apple evolution. Interestingly, we found the expression levels of GmTIFYs were significantly upregulated in our previous transcriptome data of drought, salt, and ABA treatments which were all TIFY10 and TIFY11 homologous genes. The transgenic Arabidopsis lines of AtTIFY10a and AtTIFY10b showed higher salt tolerance than WT plants (Zhu et al., 2014). We speculated that the GmTIFY10 and GmTIFY11 homologous genes may be involved in responding to salt stress.

In the soybean genome, 38 GmTIFY genes are irregularly distributed on 16 chromosomes (Figure 2). The duplication events analysis demonstrated that nine GmTIFY genes were singletons and 29 GmTIFY genes experienced gene duplication events, including eight genes with segmental duplication and 21 genes with dispersed duplication (Figure 3). The most duplication events occurred in the JAZ subfamily (29 times). Only one duplication event occurred in the PPD subfamily. Further analysis demonstrated that large segmental chromosome duplication events occurred between chromosomes 4/6 and chromosomes 10/20. The most duplication events were identified in the TIFY10 homologous genes in the JAZ subfamily (15 times).

To determine the significance of GmTIFY family genes during evolution, the EasyCodeML software was used to test the selective pressure of GmTIFY genes including purification selection, positive selection, and negative selection (Table 2). Since the TIFY1 and TIFY2 homologous genes are always grouped together in the evolutionary tree, we calculated their total ω. The positive selection of GmTIFY1, GmTIFY2, GmTIFY3, GmTIFY7, and GmTIFY10 homologous genes exceeded 1, which indicated that these TIFY homologous genes experienced positive selection during their evolutionary history.

The structural characteristics of GmTIFY family genescan be obtained by analyzing the phylogenetic tree, motifs, and positions of exons and introns (Figure 4). These results showed that genes belonging to the same phylogenetic group have similar motifs and exon/intron structures. According to previous studies of TIFY genes, GmTIFY family genes were divided into four subfamilies (TIFY, JAZ, PPD, and ZML) with different structural features. Five conserved motifs were found based on the analysis of TIFY protein sequences (Figure 4B). Motif 1 was identified as TIFY domain and distributed in all genes. Motif 2 (Jas domain) was contained by the JAZ and PPD subfamily genes. Motif 3 and 4 were identified as the sequence of CCT domain and the ZnF_GATA domain in ZIM subfamily, respectively. The PPD subfamily contains motif 1, 2, and 5, while motif 2 of the PPD subfamily lacks a PY motif which different from the original Jas domain. There was no PPD motif sequence information in the Pfam database. The conservative motif 5 domain sequence was constructed with PPD subfamily genes which were previously annotated in Arabidopsis, apples, and grapes. The motif 5 was confirmed with sequences from the report of earliest defined PPD genes (Bai et al., 2011). Conserved sequences were submitted to the MEME online tool to generate the domain logo (Figure 4D). Analysis of the promoters through PlantCARE website revealed that GmTIFY family genes contain many abiotic stresses responsive cis-elements, such as ABA-responsive element (ABRE), MYB banding site (MBS), and methyl jasmonate-responsive element (MeJA element) (Supplementary Figure 1).

To obtain the expression profiles of soybean GmTIFY genes in different tissues, RNA-seq data were downloaded from the Phytozome database and visualized by using the TBtools software (Figure 5). Results demonstrated that GmTIFY genes had significant expression differences in multiple tissues. The expression levels of most genes are relatively low in all tissues or are not even expressed. The homologous genes of TIFY3, TIFY5, TIFY10, and TIFY11 in the JAZ subfamily showed higher expression in the roots, stems, and leaves.

We analyzed the transcriptome data of GmTIFY genes under drought, salt and ABA treatment (Figure 6). Most TIFY10 and TIFY11 homologous genes significantly responded to salt treatment. From these homologous genes, six significantly upregulated candidate genes (GmTIFY10a, GmTIFY10e, GmTIFY10f, GmTIFY10g, GmTIFY11a, and GmTIFY11b) were selected and confirmed their expression level under drought and salt stress by RT-qPCR (Figure 7). The six genes were not significantly changed after drought treatment. After salt stress, the expression levels of six genes were significantly upregulated (>10-fold), especially GmTIFY10e and GmTIFY10g, which reached the peak at 1 h.

The results of the phylogenetic, duplication events, selective pressure, and expression profiles of GmTIFY genes indicated that GmTIFY10e and GmTIFY10g could play an important role in soybean. We selected GmTIFY10e and GmTIFY10g for further study, both of which significantly responded to salt stress. To determine the subcellular localization of GmTIFY10e and GmTIFY10g, we transformed the recombinant 16318hGFP vector linked to the GmTIFY10e and GmTIFY10g into Arabidopsis protoplasts using the PEG4000-mediated method, respectively (Figure 8). Both GmTIFY10e-hGFP and GmTIFY10g-hGFP fusion proteins are located in the nucleus.

To investigate the function of GmTIFY10e and GmTIFY10g in plant salt tolerance, we obtained transgenic Arabidopsis lines with high expression levels of GmTIFY10e and GmTIFY10g, respectively. For root length assay, the 5-day-old Arabidopsis plants were transferred to MS medium containing 125 mM NaCl and we calculated the primary root lengths and fresh weights. Under normal conditions, there was no significant difference between the WT plants and transgenic lines. However, the growth of the WT plants was significantly repressed when treated with 125 mM NaCl and the transgenic lines showed better growth than WT plants, with longer root lengths and heavier fresh weights (Figures 9A,C,D). To verify the salt tolerance of plants in soil, 21-day-old plants were subjected to salt stress (Figure 9B). There is no significant difference between the transgenic and WT plants under normal growth conditions. After 14 days of 250 mM NaCl treatment, the transgenic lines showed higher salt tolerance than WT, when WT is seriously wilted or even dead.

Salt stress can reduce the scavenging function of reactive oxygen species (ROS) in cells (Raza et al., 2021). Accumulation of ROS will lead to membrane lipid peroxidation, forming MDA, and activating oxygen enzymatic scavenging system including CAT, POD (Mhamdi and Van Breusegem, 2018; Wang et al., 2021a; Youssef et al., 2021). At the same time, plant cells will accumulate a large number of PRO under stress to maintain normal cell swelling pressure, prevent excessive water loss of protoplasm and enhance the adaptability of plants to adversity (Wang et al., 2021b). We measured the MDA, PRO, CAT, and POD contents of transgenic and WT lines under normal and salt treatment conditions (Figures 9E–H). Compared with WT plants, the contents of PRO, CAT, and POD in the transgenic lines were significantly increased, while the MDA contents of the transgenic lines were significantly reduced. These results all indicated that GmTIFY10e and GmTIFY10g play a role in improving the tolerance of salt stress.

To verify the function of GmTIFY10e and GmTIFY10g in response to salt stress in soybean, we obtained overexpressing plants (OE), empty vector plants (EV) and RNAi plants with A. rhizogenes-mediated transformation of soybean hairy roots and treated with 250 mM NaCl. RT-qPCR analysis demonstrated that the expression levels of OE plants were significantly higher than that of the EV, while the expression levels of RNAi plants were lower than that of the EV (Supplementary Figure 3). Phenotypic identification showed that there was no significant difference among the OE plants, RNAi plants, or EV plants under normal growth conditions. After 5 days treatment of salt, the leaves of the RNAi plants turned yellow and were more wilted than the EV plants. In the OE lines, only a few bottom leaves showed yellowing and wilting, which indicated significantly increased salt tolerance (Figures 10A,B). It indicated that overexpression of GmTIFY10e and GmTIFY10g improved the salt stress tolerance in soybean.

Reactive oxygen species is the most important signal substance for plants to respond to abiotic stresses. A large amount of ROS will accumulate in plants under abiotic stresses. NBT staining can reveal the levels of ROS accumulation in plants (Figures 10C,D). We measured the levels of ROS accumulation of plants roots. After salt treatment, the staining levels of OE plants were significantly lighter than that of EV plants and the RNAi plants showed more significant NBT stains. The same results were obtained from staining the leaves of OE plants, RNAi plants, and EV plants (Supplementary Figures 4A,F).

We further measured the contents of MDA, PRO, CAT, and POD in the roots of OE plants, RNAi plants and EV plants (Figures 10E–H). Under normal growth conditions, the MDA, PRO, CAT, and POD contents in OE plants and RNAi plants were not significantly different than in EV plants. After salt treatment, the PRO, POD, and CAT contents in OE plants were significantly higher than in EV plants and RNAi plants, while the MDA content in OE plants was lower than in EV plants and RNAi plants. We obtained similar results by measuring the MDA, PRO, CAT, and POD contents of the leaves (Supplementary Figures 4C–F). These results further confirmed that overexpression of GmTIFY10e and GmTIFY10g can enhance salt stress tolerance in soybean, which is consistent with the results of transgenic Arabidopsis.

Previous transcriptome data demonstrated that the expression levels of GmTIFY10e and GmTIFY10g were upregulated under ABA treatment. Enrichment analysis of GmTIFY10e and GmTIFY10g co-expression genes demonstrated that GmTIFY10e and GmTIFY10g can participate in the JA signal pathway and the ABA signal pathway (Supplementary Figure 5). In previous studies, the JAZ proteins were identified to participate in the ABA-dependent signal pathway through the target protein MYC/MYB transcription factor or ABA signal receptor PYL4 (Fu et al., 2017; Luo et al., 2020). To study the possible mechanisms regulated by GmTIFY10e and GmTIFY10g during salt stress responses, we analyzed the expression levels of GmCAT1, GmPOD, and GmERF115 which were genes in the ABA-mediated osmotic stress signals, and the GmSnRK2, GmPP2C, and GmMYC2 which were key genes in the ABA signal transduction pathway (Figure 11). Under normal growth conditions, the expression levels of GmCAT1, GmPOD, and GmERF115 genes were higher in OE plants and were lower in RNAi plants compared to EV plants, while GmSnRK2, GmPP2C, and GmMYC2 were higher in RNAi plants. For salt stress, the expression levels of GmCAT1, GmPOD, and GmERF115 were significantly upregulated in OE plants and downregulated in RNAi plants, while GmSnRK2, GmPP2C, and GmMYC2 significantly upregulated in RNAi plants and downregulated in OE plants. These results indicated that GmTIFY10e and GmTIFY10g may be involved in responding to salt stress through the ABA regulation pathway.