The genome data were downloaded from the Soybean database (https://www.soybase.org/dlpages/). The AtTPSs and OsTPSs protein sequences were downloaded from TAIR and RGAP, respectively. The hidden Markov Model (HMM) file of PF01397, PF03936 and PF19086 were downloaded from InterPro database (Paysan et al., 2023). Utilizing HMMER 3.0, we screened for TPS proteins within soybean (E-value <= 1e-5, similarity > 50%) (Mistry et al., 2013). Additionally, we employed the BLASTP method (Camacho et al., 2009) to search GmTPS protein sequences using AtTPS and OsTPS proteins as references (E-value <= 1e-5, similarity > 50%). Subsequently, the identified candidate TPS protein sequences underwent domain verification, adopting the analysis approach outlined by Xu et al. (2017). We filtered the longest transcript using the R package seqfinder (https://github.com/yueliu1115/seqfinder).

We employed the Muscle software (Edgar, 2004) for multiple sequence alignment of GmTPS, AtTPS and OsTPS proteins, and constructed a phylogenetic tree using IQ-TREE (Nguyen et al., 2015). The tree was visualized using the R package ggtree (Yu et al., 2017).

The 2 kb promoter region sequences upstream of the GmTPS gene were extracted using a Python program and submitted to the PlantCare (https://bioinformatics.psb.ugent.be/webtools/plantcare/html/) database for cis-element prediction (Lescot et al., 2002). All results were visualized in R software.

The chromosomal distribution of GmTPS genes was derived from the soybean genome annotation information. MCScanX software (Wang et al., 2012) was utilized for gene duplication and colinearity analysis, identifying duplication types such as tandem (TD), and whole genome (WGD). For interspecies collinearity analysis and visualization, JCVI software (Tang et al., 2008) was employed. JCVI software was used for interspecies collinearity analysis and visualization (Tang et al., 2008). ClustalW software was used to align the protein sequences and CDS sequences of TPS genes with gene duplication (Thompson et al., 2003). KaKs_Calculator software was used to calculate the synonymous substitution rate (synonymous, Ks), nonsynonymous substitution rate (nonsynonymous, Ka) and evolutionary ratio (Ka/Ks) between TPS genes duplicate gene pairs (Zhang, 2022).

The GmTPS proteins interaction network were predicted based on the AraNet2 database (Lee et al., 2015).

The transcriptome data of soybean under different tissues and development stages from the Soybean database (https://www.soybase.org/dlpages/). The expression data were visualized using the R package Pheatmap (Kolde and Kolde, 2015).

The CDS sequence of the Glyma.07G192800 and Glyma.15G263300 were cloned into the pGBKT7 vector (BD-TPS); the CDS sequence of the Glyma.12G140600, Glyma.09G029400, Glyma.20G074400 and Glyma.01G153300 were cloned into the pGADT7 vector (AD-TPSs or AD-p450s). Yeast transformants with empty pGBKT7 and AD-TPSs or AD-P450s; yeast transformants with empty pGADT7 and BD-TPSs were used as the negative control. The positive control: AD-T + BD-53. Yeast transformants with AD-TPSs, BD-TPSs and AD-p450s were used to identify the TPSs interact with other TPS proteins or P450 proteins.

Here, a total of 36 TPS members in soybean were identified by HMMER and BLASTP methods ( Table 1 ). The length of the GmTPSs protein sequence ranged from 144 (scaffold_311) to 835 (Glyma.08G163900.1); the molecular weight ranged from 16.3 (scaffold_311) to 95.6 KDa (Glyma.13G183600.1); the PI ranged from 4.27 (Glyma.13G304700.1) to 8.45 (Glyma.20G248300.1) ( Table 1 ). It’s worth noting that 86% of GmTPS have isoelectric points less than 7 ( Table 1 ), it is suggest that most GmTPS genes may are acidic proteins.

To deeper understanding the evolutionary dynamics of GmTPSs, we constructed a phylogenetic tree that comprehensively encompasses 36 GmTPSs, along with 33 AtTPSs and 53 OsTPSs. The TPSs members could be grouped into Group1, Group2, Group3, Group4 and Group5 ( Figure 1 ). Group5 contained the largest number of 35 TPSs, while Group3 contained the smallest number of 10 TPSs. Group1, Group2, and Group4 contained the number of 19, 27 and 32 TPSs, respectively. Interestingly, no GmTPSs and AtTPSs member was found in the subgroup of Group5, and Group1 and Group4 contain only GmTPS and AtTPS members ( Figure 1 ). These observations provide insights into the evolution of the TPS gene family.

To comprehend the diversity of GmTPS genes, we analyzed the gene structures, and conserved domain and conserved motif in GmTPS, the results are presented in Figure 2 . Most of the GmTPS members contained two conserved domains, the Terpene_synth domain at the N-terminal and the Terpene_synth_C domain at the C-terminal ( Figure 2 ). It is worth noting that Glyma.U032900 and Glyma.07G192800 consist only of the Terpene_synth_C domain; Glyma.19G156800 consist only of the Terpene_synth domain. Interestingly, Glyma.06G291800 contain two Terpene_synth_C domains ( Figure 2 ). Except for the Terpene_synth_C domains, Glyma.12G179500 also contains Terpene_syn_C_2 domain.

On the other hand, in the same subgroups, we found the conserved motifs of GmTPS were similar, although variations existed among different subfamilies ( Figure 2 ). For example, Group4 Subgroup members Glyma.12G138600, Glyma.12G140600, Glyma.12G138800, Glyma.12G216200 and Glyma.13G285200 contain conserved motif 1, motif 2, motif 3, motif 5, motif 6, motif 7, motif 10, motif 11, motif 14 and motif 15; Glyma.12G197500 and Glyma.12G197400 contain conserved motif 1, motif 2, motif 3, motif 4, motif 6, motif 7, motif 9, motif 10 motif 11 and motif 14 ( Figure 2 ). These diverse motifs reflect the functional diversity of GmTPS proteins. On the other hand, we found that most TBS genes contain multiple introns, except for Glyma.U032900 and Glyma.13G304800, which contain only one intron ( Figure 2 ). The variations of motifs and gene structure may contribute to the diverse biological functions of GmTPSs.

According to GFF files, we analyzed the gene distribution of 36 GmTPS. The 36 GmTPS genes were distributed on 10 chromosomes, while no GmTPS genes were distributed on chromosomes Chr1, Chr2, Chr4, Chr5, Chr10, Chr11, Chr14, Chr16, Chr17 and Chr18 ( Figure 3 ). Notably, the most GmTPS genes were distributed at chr12 and chr13, with 10 and 19 respectively, forming gene clusters. Both Chr9 and Chr15 contain one GmTPS members, while both Chr3, Chr6, Chr8, Chr19 and Chr20 contain two GmTPS members ( Figure 3 ).

TD and WGD drive the expansion of the gene family (Freeling, 2009; Panchy et al., 2016). Therefore, we explored the duplication events of TPS genes in soybean, Arabidopsis and rice. In this study, 16 WGD gene pairs and 18 TD gene pairs were confirmed in soybean, Arabidopsis and rice ( Figures 3A–C ; Supplementary Table 1 ). Overall, in soybean, Arabidopsis and rice, 10 (27.78%), 18 (60%) and 16 (30.19%) TPS genes were confirmed to be TD, and 23 (63.89%), 2 (6%) and 0 (0%) TPS genes were found to be WGD, respectively ( Figures 3A–C ; Supplementary Table 1 ). These data show that both TD and WGD contributive to the expansion of TPS genes in soybean, and WGD played a major role. However, in Arabidopsis, TD and WGD both promoted the expansion of TPS genes, and TD plays a leading role. Interestingly, in rice, only TD replication events were found. These data suggest that TPS genes expand differently in soybean, Arabidopsis and rice.

To deeper investigate the homology of the TPS gene family in G. max, we conducted a comparative analysis of TPS gene collinearity between G. max and two model organisms, Arabidopsis and rice. Our findings revealed that 2 AtTPS and 1 OsTPS were homologous gene pairs with GmTPS ( Figure 4 ). To gain insights into the evolutionary pressures of TBS genes in soybean, Arabidopsis and rice, we employed DnaSP software to calculate Ka/Ks ratios. Here, we found that the Ka/Ks value of all TPS duplication gene pairs is less than 1 in soybean, Arabidopsis and rice ( Supplementary Table 1 ). These results suggest that GmTPS, AtTPS and OsTPS genes were under purifying selection.

To unravel the may regulatory mechanisms of GmTPS genes, we analyzed their 2k promoter regions, uncovering a diverse array of 21 cis-elements ( Figure 5 ). These elements encompass various functional categories, including light response-related cis-elements (CAAT-box, Box-4, e.g); plant growth and developmental; phytohormone (ABRE, CGTCA-motif, TGACG-motif and TCA-element) and stress response related cis-elements (ARE, as-1, WUN-motif, MBS and TC-rich repeats) ( Figure 5 ). Interestingly, Group4 GmTPS genes contain ABRE, CGTCA-motif and TGACG-motif cis-elements ( Figure 5 ), hinting that GmTPS genes may be involved in ABA and JA signaling pathways. In addition, Glyma.13G285100 contains CGTCA-motif, TGACG-motif and TCA-element, suggesting that it may antagonistically participate in SA and JA signaling pathways ( Figure 5 ). It’s worth noting that Glyma.19G157000 contains a large number of light, plant growth and developmental and stress response cis-elements, while no hormone response cis-elements are found. Furthermore, our analysis revealed that Glyma.07G192800, Glyma.12G138100, Glyma.15G263300, and Glyma.09G122500 contain varying numbers of MBS cis-elements, indicative of potential roles in drought signaling pathways ( Figure 5 ). The above data indicates that GmTPS may have complex regulatory functions.

To gain deeper insights into the functional roles and regulatory intricacies within the GmTPS gene family, we leveraged the AraNet2 database to analyze and predict a protein-protein interaction network. This analysis revealed extensive interconnectivity among the GmTPS proteins, with nearly all members engaging in interactions ( Figure 6 ). Interestingly, we found that some GmTPS proteins can interact with P450 proteins, suggesting potential functional crosstalk or coordinated activities between these two protein classes ( Figure 6 ).

In addition, we selected Glyma.07G192800 and Glyma.15G263300 proteins to verify whether TPS proteins can interact with other TPS proteins or P450 proteins. To this end, we used the Y2H assay. As data shown in Figure 7 , we found that Glyma.07G192800 can interact with Glyma.12G140600, and Glyma.15G263300 can interact with Glyma.20G074400 ( Figure 7 ). Moreover, we found that Glyma.07G192800 can interact with P450 proteins Glyma.09G029400, and Glyma.15G263300 can interact with P450 proteins Glyma.01G153300 ( Figure 7 ). These data indicate that TPS protein may form heterodimer to function, or may form complex with P450 protein to function.

To further explore the functions of GmTPS, we analyzed its expression pattern based on RNA-seq data. Glyma.06g45780, Glyma.12g16940, Glyma.12g32370 and Glyma.07g30700, expressed higher in flower than other tissues. We speculate that these GmTPS genes may synergistically regulate soybean growth and development ( Figure 8 ). And Glyma.12g16940 and Glyma.12g32370 are expressed only in flower. Additionally, these GmTPSs no expression of Glyma.13g32380, Glyma.17g05500, Glyma.12g16990, Glyma.13g38050, Glyma.06g45780, Glyma.12g16940, Glyma.12g32370, Glyma.07g30700, Glyma.07g30710 and Glyma.20g18280 are detected in the organs of seed and nodule. On the other hand, Glyma.08g17470 and Glyma.15g41670 are widely expressed in root, young leaf, pod shell, flower, seed and nodule ( Figure 8 ). The result indicate that GmTPS genes may involve in diverse aspects of plant growth and developmental processes.