Soybean genome, protein, complementary DNA (cDNA) sequences, and gene annotation files were obtained from NCBI¹ and the Phytozome database² (Finn et al., 2011; Fernandez-Pozo et al., 2015). The Hidden Markov Model (HMM) profile corresponding to the TLP Tub domain (PF01167) from the Pfam protein family database³ was used to identify potential TLPs in the soybean genome (G. max Wm82. a2.v1) using HMMER v3 (Eddy, 1998; Mistry et al., 2021). Finally, the presence of the Tub domain in each TLP protein sequence was confirmed with the SMART tool⁴ (Letunic et al., 2012) and Pfam database. The molecular weight and isoelectric point data for GmTLPs were calculated by ExPASY⁵ (Artimo et al., 2012). Subcellular localization was predicted with WoLF PSORT⁶.

The full-length amino acid sequences of TLP members in rice (OsTLPs), Arabidopsis (AtTLPs), cotton (GhTLPs), maize (ZmTLPs), apple (MdTLPs), poplar (PtTLPs), wheat (TaTLPs), tomato (SlTLPs), and the newly identified GmTLPs were obtained from NCBI and Phytozome, respectively, and aligned with default parameters using ClustalW (Chenna et al., 2003). An unrooted phylogenetic tree was constructed using the neighbor-joining (NJ) method in MEGAX (version 10.1.8) (Tang et al., 2021) with the following parameters: pairwise deletion; Poisson model; 1000 bootstrap replications.

The amino acid sequence of 22 TLP proteins of soybean aligned using DNAMAN (version 6.0.3).

Chromosomal location data for GmTLPs were obtained from the Phytozome database. Intron insertion sites were identified by comparing the coding sequence of each TLP gene with the corresponding full-length sequence using the Gene Structure Display Server (GSDS) 2.0⁷ (Hu et al., 2015). The conserved domain of the identified soybean GmTLP protein sequences were determined using MEME⁸ with the maximum number of motifs set to 10 (Bailey et al., 2009).

Protein sequences of the 22 GmTLPs were submitted to the SMART website⁹ to obtain data related to conserved protein domains, and GSDS 2.0 was used for visual analysis. Three-dimensional models of the Tub domain were built with SWISS-MODEL¹⁰ (Dong et al., 2019). Tub domain models were obtained for 20 GmTLPs with the protein sequence identity set to ≥30%.

Soybean gene expression files were downloaded from the Soybase website¹¹ to analyze the expression patterns of 22 GmTLPs members in different tissues at different developmental stages under normal conditions, including young-leaf, flower, pod, pod shell, seed, root, and nodule. In the database file provided by Soybase website, only 18 members’ tissue differential expression information were found for further analysis. Transcriptome data for GmTLPs members under various abiotic stresses from our previous studies (Wang et al., 2020). 22 GmTLPs members were used for searching in transcriptome data, and their expression levels under normal condition, ABA treatment, drought and salt stresses were analyzed. Finally, the relevant information of 21 members was obtained. TBtools (version 1.075) (Chen et al., 2020) was used for visualization and cluster analysis of GmTLP expression patterns.

GmTLP sequences obtained from the Phytozome database were extracted in batches with TBtools, and the 2000 bp upstream promoter sequences of the 22 GmTLP genes were obtained and submitted to the online program PlantCARE¹² to identify cis-acting elements. GSDS 2.0 was used for data visualization.

The soybean variety Zhonghuang39 was used for analysis of GmTLP gene expression in this study. Soybeans were grown in 1: 1 vermiculite: humus in a greenhouse with a 16/8 h light/dark cycle, day/night temperatures of 28/20°C, and a relative humidity of 70%. At 14 days, the seedlings at the four-leaf stage were stressed with drought or salt. Referring to previous research methods (Wang et al., 2021), soybean seedlings were removed from soil. For drought stress, the seedlings were placed on filter paper; for salt stress, the seedlings were immersed in 200 mM NaCl solution. The sampling time of drought or salt stress was 0, 0.5, 1, 2, 4, 8, 12, and 24 h. There were three biological replicates per treatment. After treatment, the leaves were frozen in liquid nitrogen and stored at –80°C before further analysis (Xu et al., 2008). These samples were used for qRT-PCR analysis of subsequent GmTLPs members.

Total RNA was extracted from soybean leaves using a plant RNA extraction kit following the manufacturer’s instructions (TIANGEN, Beijing, China). cDNA was synthesized using the PrimeScript™ RT Reagent Kit (TaKaRa, Shiga, Japan) following the manufacturer’s protocol. Primers (Supplementary Table 2) were designed using Primer Premier 5.0. The soybean Actin gene (U60506) was used as the internal control for quantitative real-time PCR (qRT-PCR). There were three technical replicates for each sample. Differential expression was determined from the relative gene expression data using the 2^{–Δ Δ CT} method (Le et al., 2011).

We constructed an expression vector labeled with green fluorescent protein (hGFP) for subcellular localization analysis. The full-length cDNA sequence of GmTLP8 was fused to the N-terminal hGFP protein driven by the CaMV35S promoter (Xu et al., 2007). The 35S:GFP vector was used as a control. A PEG4000-mediated method was used to transform the GmTLP8-GFP recombinant plasmid into Arabidopsis protoplasts (He et al., 2016). After incubation for 18–20 h in the dark at 22°C, the nucleus of GmTLP8-GFP protoplasts were specifically stained with 4’, 6-diamidino-2-phenylindole (DAPI). The fluorescence signal was observed using a confocal laser scanning microscope (Zeiss LSM 700, Oberkochen, Germany). There were three technical replicates for each group.

The transformation was conducted to produce soybean hairy roots that were characterized by overexpression of GmTLP8 (GmTLP8-OE), RNA interference of GmTLP8 (GmTLP8-RNAi), or with the empty pCAMBIA3301 vector (EV-Control) (Chen et al., 2021). The CDS of GmTLP8 was amplified without stop codon using gene-specific primer pairs, under the control of the CaMV35S promoter, GmTLP8 cDNA was ligated into the plant transformation vector pCAMBIA3301 to generate GmTLP8-overexpressing (GmTLP8-OE) vector. In order to construct the RNAi vector, a 564 bp fragment including the first intron sequence and its reverse complement was synthesized (Biomed, Beijing, China) and inserted into pCAMBIA3301 to generate the pCAMBIA3301-GmTLP8-RNAi (GmTLP8-RNAi) vector. The recombinant construct and the empty pCAMBIA3301 (EV-Control) vector were transferred into A. rhizogenes strain K599, as previously described, then injected into soybean (G. max cv. Zhonghuang39) hypocotyl for A. rhizogenes-mediated transformation of soybean hairy roots (Wang et al., 2009; Du et al., 2018).

The injected plants were placed in a high-humidity greenhouse until hairy roots were generated at the infected site and had grown to ∼5 cm in length. After cutting off the original tap root 0.5 cm below the infected site, the seedlings were transplanted into fertilized soil and cultivated in a greenhouse at 25°C with a 16/8 h light/dark photoperiod for 7 days (Yu et al., 2021). The qRT-PCR analysis of GmTLP8 expression in GmTLP8-OE, EV-control and GmTLP8-RNAi transgenic hairy root plants before processing (Supplementary Figure 2G). Each of the hairy root-related experiments was replicated at least three times independently. The primers of GmTLP8-3301-F and GmTLP8-3301-R were listed in Supplementary Table 2.

Transgenic hairy root composite soybean plants were used in drought and salt stress assays after 7 days of normal growth. For drought treatment, soybean plants were grown for 7 days without watering; for NaCl treatment, soybean plants were treated with 150 mM NaCl for 3 days. Drought and salt treatment experiments were conducted a minimum of three times. Both the treated and untreated soybean hairy roots were washed with water prior to RNA isolation and physiological/biochemical experiments.

Several physiological parameters were measured in transgenic GmTLP8-OE, EV-Control, and GmTLP8-RNAi lines after the drought and NaCl treatments, namely levels of proline (Pro), malondialdehyde (MDA), hydrogen peroxide (H₂O₂), superoxide anion (O₂^–), and chlorophyll. Measurements were taken in soybean leaves using appropriate assay kits (Cominbio, Suzhou, China) following the manufacturer’s instructions. All measurements were performed in three biological replicates.

Leaves from the three transgenic lines were stained with 3,3-diaminobenzidine and nitro blue tetrazolium after drought or salt stress treatment. The leaves were immersed in DAB solution or NBT staining solution (Solarbio, Beijing, China) for 18 or 14 h, respectively. Samples were then destained in a boiling solution of 3: 1 anhydrous ethanol: glycerol until the leaves were white (Du et al., 2018). Images were taken using a Canon 50D camera (Canon, Tokyo, Japan). There were three biological replicates for each plant line–treatment group combination.

Twenty-two GmTLP family members were identified in this study. The SMART and Pfam databases were used to confirm the presence of the conserved Tub domain in all of the putative TLP proteins. Twenty-two GmTLP genes were unevenly distributed across 13 chromosomes of soybean. According to their positions on chromosomes, we named them GmTLP1 to GmTLP22. The details of TLPs in soybean, such as the coding sequence (CDS) length, amino acid length (aa), molecular weight (MW), isoelectric point (pI), and subcellular location are shown in Table 1.

Among the 22 GmTLPs, the protein length ranged from 183 (GmTLP5) to 454 amino acids (GmTLP21). The minimum protein MW was 20.8 kDa (GmTLP5), and the maximum was 50.6 kDa (GmTLP21). The pI ranged from 7.54 in GmTLP10 to 9.65 in GmTLP8. Twelve of the proteins were predicted to be located in the nucleus, five in the cytosol, four in the mitochondria, and two in the chloroplast (Table 1), with GmTLP1 predicted to be located in either the nucleus or cytosol.

A physical location map of the GmTLPs was drawn using physical location data from the soybean genome. The 22 GmTLP genes were distributed across 13 chromosomes, which were chromosome 1, 2, 7, 8, 10, 11, 12, 13, 14, 15, 16, 17, 20, respectively. There were three genes on chromosomes 2 and 13, two genes on chromosomes 7, 12, 15, 16, and 17, and only one gene each on chromosomes 1, 8, 10, 11, 14, and 20 (Figure 1).

To reveal the phylogenetic relationships between TLPs in different plant species, an unrooted phylogenetic tree was constructed by comparing the amino acid sequences for all of the known TLP members in several species, totaling 132 proteins. There were 11 from Arabidopsis (Lai et al., 2004), 15 from maize (Chen et al., 2016), 14 from rice (Liu, 2008), four from wheat (Hong et al., 2015), 11 from tomato (Zhang et al., 2020), nine from apple (Xu et al., 2016), 11 from poplar (Yang et al., 2008), and 35 from cotton (Li Z. et al., 2021) in addition to the 22 putative TLPs identified in soybean. Phylogenetic tree was divided into five groups based on protein homology, and there were one, two, five, six, and eight GmTLP members in groups I, II, III, IV, and V, respectively (Figure 2).

The results of multiple sequence alignment showed that the positions of F-box domain and Tub domain in GmTLP protein sequence were located in the front and rear segments of the sequence (Figure 3).

We analyzed the gene structure of the 22 GmTLPs using GSDS 2.0 online to determine the intron and exon distribution of each (Figure 4A). A total of 15 genes contained three introns and four exons each, and the other seven genes contained four introns and five exons.

A total of 10 conserved motifs (E ≤ 0.01) were analyzed using the MEME website to explore conservation and diversity of soybean TLPs. Among the 22 GmTLP family members, 19 contained all 10 motifs, with GmTLP4 containing two copies of motif 3 and GmTLP20 containing two copies of motif 5. GmTLP5, GmTLP10, and GmTLP14 contained six, five, and two motifs, respectively (Figure 4C). Consensus sequences for putative motifs are shown in Supplementary Figure 1.

From the Pfam database, we found that two conserved domains in GmTLPs were Tub (PF01167) and F-box (PF00646). We then analyzed the conserved protein domains using both the SMART website (Figure 4D) and homology modeling in SWISS-MODEL (Figure 5).

The results of conserved domain analysis showed that 20 of the GmTLPs contained one F-box domain in the N-terminal region and one Tub domain in the C-terminal region. GmTLP4 had a Tub domain in the C-terminal region but no F-box domain in the N-terminal region, and GmTLP20 contained one F-box domain in the N-terminal region and two Tub domains in the C-terminal region. The analysis also revealed one ANK domain in the C-terminal region of GmTLP10, the function of which was not clear (Figure 4D).

The homology modeling was a useful tool for the prediction of protein structure, and protein structural information was often more valuable than sequence data alone in determining protein function. We generated three-dimensional (3D) models of the Tub domains for 20 of the GmTLPs. These models showed that the Tub domain of each GmTLP was closed by a β-barrel with 12 anti-parallel strands and a central hydrophobic α-helix, which is a typical structure for a Tub domain (Figure 5).

Cis-acting regulatory elements play an important role in modulating gene expression. To understand transcriptional regulation of GmTLPs, we identified cis-acting elements within the promoter region of each GmTLP gene, defined as the 2000 bp region upstream of the start codon (Figure 4B and Supplementary Table 1). Results showed that most cis-acting elements in GmTLP promoters were involved in hormone or stress responses.

The main hormone-related cis-acting elements identified were ABA response element (ABRE), TATC-box (gibberellin), AuxRR core (auxin), TCA element (salicylic acid), and the TGACG motif (methyl jasmonate). Among the GmTLP promoters, 14 genes contained ABRE, four contained a TATC-box, six contained the AuxRR core, 11 contained a TCA element, and 12 contained a TGACG motif. This indicated that the GmTLPs may be involved in hormone-related responses.

The abiotic stress cis-acting elements identified were as follows: anaerobic inducing element (ARE), low temperature response element (LTR), a MYB binding site involved in drought induction (MBS), a drought and salt response element (MYC), and a defense and stress response element (TC-rich element). MYC, which was previously reported to be involved in drought and salt stress-induced responses, was revealed to be distributed in all of the GmTLP promoter sequences. In addition, 20 GmTLP promoters contained ARE, seven contained LTR, seven contained MBS, and six contained TC-rich elements. The presence of these cis-acting elements related to abiotic stresses indicated that GmTLPs are abiotic stress-responsive.

Among the 22 GmTLPs, 14 contained ABREs, indicating that these genes can be regulated by ABA. Drought and salt response element MYCs were present in all 22 GmTLP promoter sequences, strongly suggesting that GmTLP members are involved in the responses to those stressors (Zuo et al., 2020). In total, cis-acting element analysis indicated that most members of the GmTLP family may be regulated by ABA in response to drought and salt stresses.

To understand expression patterns of TLPs during the growth and development of soybean and throughout different plant tissues, publicly available transcriptome sequencing data from the SoyBase database were analyzed. For 18 GmTLP members, we analyzed gene expression levels in different plant tissues including young-leaf, flower, pod, pod shell, seed, root, and nodule. The results showed that GmTLP7 was extremely high expressed in all seven tissues. GmTLP6, 10, 22, 11, and 21 were expressed at extremely low levels or not expressed in seven tissues. GmTLP2, 4, 9, 12, and 14 were expressed in some tissues, but not in others. GmTLP13, 17, 3, 5, 20, 15, and 16 were expressed in all seven tissues, with extremely high expression in some tissues and extremely low expression in others. Phylogenetic analysis divided the 18 GmTLP members into different groups, and members within each group shared similarities at the expression level in all seven tissues (Figure 6). The results showed a great deal of spatiotemporal difference in GmTLP expression levels.

We used a previously published transcriptome sequencing database to quantify the expression of GmTLPs under normal condition, ABA treatment, drought and salt stresses (Wang et al., 2020), and screened 21 GmTLP members (Figure 7). The results showed that GmTLP21 was up-regulated and GmTLP14 was down-regulated under ABA treatment. Under drought stress, six genes were up-regulated, and GmTLP3, 8, 11, 13, 19, and 22, four genes were down-regulated, and GmTLP14, 15, 18, and 20, respectively; Under salt stress, GmTLP8 was up-regulated and GmTLP14 was down-regulated. The P-value of data in abiotic stress expression profiles is shown in Supplementary Table 3.

According to the expression profiles of GmTLPs under different abiotic stresses, five genes (GmTLP8, 11, 13, 19, and 22) up-regulated expression under drought stress and one gene (GmTLP8) under salt stress were selected for qRT-PCR analysis to further verify their relative expression levels under drought and NaCl treatments. The selected genes were GmTLP8, 11, 13, 19, and 22. Under drought treatment (Figure 8A), GmTLP8 expression peaked at 8 h (with an 8.5-fold increase compared to 0 h), GmTLP11 at 24 h (3.7-fold), GmTLP13 at 1 h (3.7-fold), GmTLP19 at 24 h (4.9-fold), and GmTLP22 at 8 h (4.8-fold). Under salt treatment (Figure 8B), GmTLP8 expression peaked at 8 h (7.3-fold), GmTLP11 at 2 h (4.5-fold), GmTLP13 at 24 h (1.8-fold), GmTLP19 at 12 h (4.4-fold), and GmTLP22 at 12 h (4.5-fold). These results showed that GmTLP8 was the most highly expressed in response to drought and salt treatments, and it was therefore selected for further study.

To determine the subcellular localization of GmTLP8, the open reading frame (ORF) sequence (excluding the termination codon of GmTLP8) was fused with the N-terminal of the humanized green fluorescent protein (hGFP) reporter and co-transformed into Arabidopsis protoplasts. A 35S:hGFP as the control, the fluorescence signal in the cells was detected by confocal laser scanning microscopy. The fluorescence of GmTLP8 was detected in the nucleus and cytoplasm, while the fluorescence of the control 35S:hGFP was observed in the whole cell. DAPI staining also showed the localization of GmTLP8 in the nucleus (Figure 9). It suggests that GmTLP8 act as a transcription factor in the nucleus (Li S. et al., 2021).

The stress-tolerant effect of GmTLP8 in soybean was explored using transgenic soybean hairy root composite plants. The hairy roots of GmTLP8-OE, EV-Control and GmTLP8-RNAi transgenic lines were used to analyze the relative expression level of GmTLP8. qRT-PCR analysis showed that the expression level of GmTLP8-OE transgenic hairy roots was significantly higher than that in EV-Control, and the expression level of GmTLP8-RNAi transgenic hairy roots was lower than that in EV-Control (Supplementary Figure 2G). Under normal growth conditions, no significant differences were observed between GmTLP8-OE, the EV-Control, and GmTLP8-RNAi lines (Figure 10A). However, after exposure to drought (Figure 10B) and salt (Figure 10C) treatments, there were significant phenotypic differences between GmTLP8-OE, EV-Control, and GmTLP8-RNAi plants. Compared with EV-Control, GmTLP8-RNAi plants showed more severe leaf dehydration and wilting stress phenotype, whereas GmTLP8-OE showed fewer rolled leaves and a delayed leaf wilting phenotype. The survival rates of GmTLP8-OE, EV-Control, and GmTLP8-RNAi lines under drought stress were 93, 67, and 40%, respectively; these survival rates were comparable to those of salt-stressed plants (Figure 10H).

Proline (Pro), malondialdehyde (MDA), hydrogen peroxide (H₂O₂), and superoxide anion (O₂^–) levels are important indicators of the effects of abiotic stresses on plant growth (Leng et al., 2021). Proline is a protective agent against osmotic stress; MDA reflects the degree of lipid oxidative damage; H₂O₂ and O₂^– play immune and signal transduction roles, although excessive accumulation may lead to cell membrane damage (Du et al., 2018; Zhang et al., 2019). Chlorophyll levels are an important indicator of plant photosynthetic capacity (Tanaka and Tanaka, 2006). To further analyze the potential physiological mechanism of GmTLP8 in plant stress tolerance, we measured the levels of Pro, MDA, H₂O₂, O₂^–, and chlorophyll in the leaves of GmTLP8-OE, EV-Control, and GmTLP8-RNAi plants under normal growth conditions and under drought or salt stress (Figures 10F,G,I–K). Levels of Pro and chlorophyll were higher in GmTLP8-OE compared with EV-Control, whereas levels of MDA, H₂O₂, and O₂^– were lower. In contrast, the GmTLP8-RNAi lines had lower Pro and chlorophyll levels but higher MDA, H₂O₂, and O₂^– levels than EV-Control.

H₂O₂ and O₂^–, produced by the reactive oxygen species (ROS) pathway in leaf cells under abiotic stress, were measured to assess the degree of damage in leaf cells (Cui et al., 2019). This was done using DAB and NBT to stain the leaves of GmTLP8-OE, EV-Control, and GmTLP8-RNAi plants (Figures 10D,E). Under normal conditions, leaves from the GmTLP8-OE, EV-Control, and GmTLP8-RNAi lines showed minimal staining, with no significant difference between lines. Under drought and salt stresses, compared with EV-Control, GmTLP8-OE leaves showed shallow staining, whereas GmTLP8-RNAi showed deeper staining. These results demonstrated that the GmTLP8-OE line had lower levels of leaf damage and the GmTLP8-RNAi line had more severe leaf damage compared to EV-Control in response to exogenous abiotic stresses. The results of staining leaves with DAB and NBT were consistent with the physiological indexes of H₂O₂ and O₂^– contents.

To analyze the potential stress tolerance mechanism of GmTLP8, genes known to be involved in drought and salt stress responses were selected, namely GmDREB1 (Kasuga et al., 1999), GmDREB2 (Chen et al., 2007), GmNAC11 (An et al., 2018), GmNCED3 (Pandey and Gautam, 2020), GmSOS1 (Ma et al., 2020), and GmWRKY27 (Wang et al., 2015) (Figures 11A–F,H–M). Expression of these genes in the hairy roots of GmTLP8-OE, EV-Control, and GmTLP8-RNAi transgenic soybean lines were measured via qRT-PCR. Plants were drought-treated by withholding water for 7 days or salt-treated with 150 mM NaCl for 3 days. Under normal growth conditions, the selected stress responsive genes were expressed at lower levels in all three plant lines compared to plants that had been exposed to drought or salt stress (Supplementary Figures 2A–F and Figures 11A–F,H–M). In plants that had been stressed, compared with EV-Control, the six stress-related genes were significantly up-regulated in GmTLP8-OE plants and down-regulated in GmTLP8-RNAi plants. These results suggest that overexpression of GmTLP8 may activate expression of downstream drought- and salt-response genes.