Identification and in silico characterization of soybean trihelix-GT and bHLH transcription factors involved in stress responses

Environmental stresses caused by either abiotic or biotic factors greatly affect agriculture. As for soybean [Glycine max (L.) Merril], one of the most important crop species in the world, the situation is not different. In order to deal with these stresses, plants have evolved a variety of sophisticated molecular mechanisms, to which the transcriptional regulation of target-genes by transcription factors is crucial. Even though the involvement of several transcription factor families has been widely reported in stress response, there still is a lot to be uncovered, especially in soybean. Therefore, the objective of this study was to investigate the role of bHLH and trihelix-GT transcription factors in soybean responses to environmental stresses. Gene annotation, data mining for stress response, and phylogenetic analysis of members from both families are presented herein. At least 45 bHLH (from subgroup 25) and 63 trihelix-GT putative genes reside in the soybean genome. Among them, at least 14 bHLH and 11 trihelix-GT seem to be involved in responses to abiotic/biotic stresses. Phylogenetic analysis successfully clustered these with members from other plant species. Nevertheless, bHLH and trihelix-GT genes encompass almost three times more members in soybean than in Arabidopsis or rice, with many of these grouping into new clades with no apparent near orthologs in the other analyzed species. Our results represent an important step towards unraveling the functional roles of plant bHLH and trihelix-GT transcription factors in response to environmental cues.


Introduction
Soybean [Glycine max (L.) Merril] is one of the most important crop species in the world. It is widely used for both human and animal consumption due to the high protein and oil contents of its grains. More recently, the potential for using soybean oil in renewable fuel production has also emerged (Programa Nacional de Produção e Uso de Biodiesel). Since it belongs to the Fabaceae family, soybean also takes part in the process of organic nitrogen fertilizer production through its symbiotic association with nitrogen-fixing bacteria (Gepts et al., 2005). Currently, soybean producers are primarily concerned with losses caused by drought stress, Asian Soybean Rust (ASR, caused by the fungus Phakopsora pachyrhizi) and soybean cyst nematode (SCN, caused by Heterodera glycines) (EM-BRAPA, 2007). Furthermore, the genetic variability found in soybean germplasm for those characteristics is restricted, which increases the vulnerability of this species to environmental stresses (Priolli et al., 2002;Miles et al., 2006).
As sessile organisms, higher plants are continuously exposed to a great variety of environmental stimuli. Because their survival depends on the ability to cope with those stimuli, plants have evolved a variety of sophisticated molecular mechanisms in response to environmental stresses. These generally involve alterations in gene expression, leading to changes in plant physiology, metabolism and developmental activities. Whether caused by abiotic (such as drought, salt and cold) or biotic factors (such as pathogens and insects), environmental stresses have serious adverse effects on agriculture. Therefore, a thorough understanding of the molecular mechanisms involved in plant stress tolerance has become pivotal for the development of new strategies and technologies related to the increasing demand on agricultural production (Rao et al., 2006;Yoshioda and Shinozaki, 2009).
Upon stimuli perception, responses of plants to environmental stresses comprise the activation of a multitude of interconnected signaling pathways (Singh et al., 2002). The phytohormones abscisic acid (ABA), ethylene (ET), jasmonic acid (JA) and salicylic acid (SA), aside from reactive oxygen species (ROS), are known to act as messenger molecules that trigger specific (but at times overlapping) pathways of this complex network, leading to the accumulation of stress-related gene products (Yoshioda and Shinozaki, 2009). Besides, a great number of studies have highlighted the importance of the transcriptional regulation of targetgenes through transcription factors in plant responses to environmental stresses (Zhou et al., 2008;Chen et al., 2009;Zhang et al., 2009). Transcription factors act by binding to cis-elements in the promoter regions of target-genes, thereby activating or repressing their expression. Transcriptional reprogramming is known to result in both spatially and temporally altered expression patterns of stressrelated genes. Thus, transcription factors are key players in fine-tuning stress responses at the molecular level (Singh et al., 2002;Eulgem, 2005).
A large part of a plant's genome is devoted to transcription. With the recent completion of the soybean genome sequencing and assembly, a comparative analysis of putative transcription factor-encoding genes found in both soybean and the model dicot Arabidopsis thaliana can be performed. In the leguminous plant (whose genome is six times larger than that of A. thaliana), over 5,600 transcription factors were identified, these corresponding to about 12% of the predicted protein-coding loci (Schmutz et al., 2010). In contrast, in the model plant the total number of transcription factors (~2,300) comprises only up to 7% of the predicted protein-coding loci (Singh et al., 2002). The overall distribution of these genes among known transcription-factor families is similar among the two genomes, although some families are relatively sparser or more abundant in soybean. Thus, even though the A. thaliana genome often serves general comparisons, differences in biological function between species might occur (Schmutz et al., 2010).
Basic helix-loop-helix (bHLH) proteins constitute one of the largest families of transcription factors. They are found in all three eukaryotic kingdoms and are involved in a myriad of regulatory processes. Members of this family share the bHLH signature domain, which consists of~60 amino acids comprising two distinct regions, a basic stretch at the N-terminus consisting of~15 amino-acids involved in DNA binding, and a C-terminal region of~40 aminoacids composed of two amphipathic a-helices, mainly consisting of hydrophobic residues linked by a variable loop (the "helix-loop-helix" region). This region is responsible for promoting protein-protein interactions through the formation of homo-and hetero-dimeric complexes (Toledo-Ortiz et al., 2003;Carretero-Paulet et al., 2010;Pires and Dolan, 2010). The Lc protein from Zea mays, reported as a transcriptional activator in the anthocyanin biosynthetic pathway (Ludwig et al., 1989), was the first plant bHLH member identified. The involvement of bHLH members in plant developmental processes (Szecsi et al., 2006;Menand et al., 2007), light perception (Liu et al., 2008), iron and phosphate homeostasis (Yi et al., 2005;Long et al., 2010;Zheng et al., 2010), and phytohormone signalling pathways (Abe et al., 1997;Friedrichsen et al., 2002;Lorenzo et al., 2004;Anderson et al., 2004;Fernandez-Calvo et al., 2011;Hiruma et al., 2011;Seo et al., 2011) has also been reported. In fact, Arabidopsis MYC2 is to date the most extensively characterized plant bHLH transcription factor, and it seems to be a global regulator of hormone signalling. MYC2 has been described as an activator of ABAmediated drought stress-response (Abe et al., 1997(Abe et al., , 2003. It also regulates JA/ET-induced genes, either as an activator in response to wounding, or as a suppressor in pathogen responses (Anderson et al., 2004;Lorenzo et al., 2004;Hiruma et al., 2011). In these cases, the activity of MYC2 is itself subject to regulation by JAZ proteins, in a SCF COI1 proteosome degradation -dependent pathway (Chini et al., 2007). Additionally, MYC2 seems to form homo-and heterodimers with two other closely-related bHLH proteins (MYC3 and MYC4), and their interaction is essential for full regulation of JA responses in Arabidopsis (Fernandez-Calvo et al., 2011).
Trihelix-GT factors constitute another family of plant-specific transcription factors. They are characterized by binding specificity for GT-elements present in the promoter region of many plant genes (Hiratsuka et al., 1994;Nagano et al., 2001) and are among the first transcription factors identified in plants ( McCarty and Chory, 2000). They share one or two trihelix (helix -loop -helix -loophelix) structures, each consisting of three putative a-helices, which are responsible for binding to DNA (Zhou, 1999). Dimerization of GT factors, or interaction between trihelix-GT and other transcription factors appear to play a major role in the regulatory function of this family (Zhou, 1999). In addition, recent studies demonstrated that posttranslational modifications may occur in at least some GTfactors, as shown for Arabidopsis light-responsive GT-1 (Maréchal et al., 1999;Nagata et al., 2010). Members of the trihelix-GT family were first described as being involved in the regulation of light-responsive genes (Green et al., 1987(Green et al., , 1988. Nevertheless, further studies in rice and Arabidopsis showed that some GT factors are not light-responsive at the transcriptional level (Dehesh et al., 1990;Kuhn et al., 1993). The involvement of this family in seed maturation (Gao et al., 2009), control of flower morphogenesis (Griffith et al., 1999;Brewer et al., 2004;, and response to environmental cues (O'Grady et al., 2001;Park et al., 2004;Wang et al., 2004;Xie et al., 2009;Fang et al., 2010) has also been reported.
In recent years, a growing number of transcription factors belonging to families, such as AP2, NAC and WRKY, have been connected to the responses of soybean against environmental stresses Pinheiro , 2009;Zhou , 2008). In addition, the involvement of two soybean trihelix-GT factors  and GmGT-2B (Glyma10g30300)] in abiotic stress toler-ance has recently been proposed, following heterologous expression in Arabidopsis (Xie , 2009). Nevertheless information regarding soybean bHLH and trihelix-GT members and their roles in this species remains scarce. In the present study we, therefore, aimed at identifying soybean bHLHand trihelix-GT-encoding genes, as well as investigating their involvement in response to environmental stresses. Given the dimension of the bHLH family in plants (with more than 600 members in Arabidopsis divided into 32 groups), we decided to focus on a single monophyletic group (subfamily 25, Carretero-Paulet et al., 2010), once we had found some interesting soybean candidates within the LGE Soybean Genome database (Nascimento et al., 2012) that belong to this group. At least 45 bHLH (from subgroup 25) and 63 trihelix-GT putative genes reside in the soybean genome. Among these, at least 14 bHLH and 11 trihelix-GT seem to be involved in responses to abiotic/biotic stresses. A phylogenetic analysis allowed us to successfully cluster these genes with members of bHLH and trihelix-GT proteins from other plant species. All together, our results represent an important step towards understanding the molecular mechanisms by which soybean responds to environmental cues.

Sequence identification and annotation
In order to identify putative soybean bHLH sequences, the TAIR (The Arabidopsis Information Resource) gene id from all 17 bHLH proteins belonging to subgroup 25 in Arabidopsis was used to search the soybean database in Phytozome and at JGI (Joint Genome Institute). Soybean peptide homologs for each A. thaliana sequence were identified from a BLASTP search with default parameters in Phytozome and redundant sequences were manually discarded. The protein sequences obtained were scanned for the existence of the bHLH domain using the SMART database. The software MEME (multiple EM for motif elicitation) version 4.4.0 was used for motif identification, using the following parameters: minimum and maximum motif width set to 6 and 50 amino acids, respectively, with any number of motif repetitions. Motif detection was restricted to a maximum of 10. Identified motifs were also compared with conserved compositions already described for bHLH sequences. In addition, the bHLH domain was manually delimited according to plant-specific boundaries, as determined by Toledo-Ortiz et al. (2003) and Carretero-Paulet et al. (2010). Classification of soybean sequences in subgroup 25 was accomplished by mismatch counting from the consensus established for A. thaliana (Carretero-Paulet et al., 2010). Sequences with more than 8 mismatches in conserved positions were discarded. Moreover, no mismatches were allowed at residues H 9 , E 13 and R 16 of the basic region, since these are crucial for DNA-binding activity, and a consensus among subgroup 25 sequences.
The identification of putative trihelix-GT protein sequences from soybean was accomplished as follows: the conserved trihelix sequence of previously reported soybean genes (O'Grady et al., 2001;Xie et al., 2009) along with motifs predicted for this family (Fang et al., 2010), were blasted (TBLASTN) against the soybean genome in Phytozome. All homologous sequences with an E-value of less than 0.0001 were scanned for the existence of the trihelix domain using SMART (domains with less significant scores than default cut-offs were also analyzed). Motif identification and comparison with conserved trihelix-GT compositions were performed using MEME. Sequences that did not fit these criteria were removed from the analysis.
To determine the intron-exon organization of all bHLH and trihelix-GT genes, the full length coding sequences were aligned with the corresponding genomic sequences available on Phytozome. Intron-exon maps of the genes were drawn using Fancy Gene v1.4 software.

Gene expression data mining
Expression profiles of the identified bHLH and trihelix-GT sequences in both biotic and abiotic situations were obtained by mining the LGE Soybean Genome database. A "gene" search was carried out using Phytozome's gene model codes and each gene had its 5' and 3' untranslated regions verified in Gbrowse. Gene expression was confirmed by database searches in NCBI ESTs and LGE superSAGE stress experiments with soybean leaves infected with Asian soybean rust (accession PI 561356, resistant) vs. uninfected leaves, and soybean roots subjected to drought (cultivar BR16, susceptible / cultivar Embrapa-48, tolerant) vs. untreated roots from both cultivars.

Phylogenetic analysis
The phylogenetic analysis of plant trihelix-GT factors was performed using protein sequences from A. thaliana, G. max, Medicago truncatula and Oryza sativa. For plant bHLH transcription factors, protein sequences from A. thaliana, G. max, O. sativa and Physcomitrella patens were used. In both cases, multiple sequence alignments were conducted with full-length protein sequences using the CLUSTALW tool (Thompson et al., 1994) implemented in MEGA ver. 4.0 (Tamura et al., 2007). The phylogenetic analysis was performed by two different and independent approaches, viz. the neighbor-joining (NJ) and Bayesian methods. The NJ method was performed within MEGA v4.0. Molecular distances of the aligned sequences were calculated according to the p-distance parameter, with gaps and missing data treated as pairwise deletions. Branch points were tested for significance by bootstrapping with 1000 replications. Bayesian analysis was conducted in MrBayes 3.1.2 software (Huelsenbeck et al., 2001;Ronquist and Huelsenbeck, 2003) with the mixed amino-acid substitution model + gamma + invariant sites. Two inde-pendent runs of 5,000,000 generations each, with two Metropolis-coupled Monte Carlo Markov chains (MCMCMC) were run in parallel, each one starting from a random tree. Markov chains were sampled every 100 generations and the first 25% of the trees were discarded as burn-in. The remaining ones were used to compute the majority rule consensus tree (MrBayes command allcompat), and the posterior probability of clades and branch lengths. The unrooted phylogenetic trees of trihelix-GT and bHLH proteins were visualized and edited using the software FigTree ver. 1.3.1.

Results and Discussion
Identification and analysis of soybean bHLH-encoding genes In the past few years several phylogenetic studies have emerged as attempts to perform the classification of bHLH proteins in plants (Heim et al., 2003;Toledo-Ortiz et al., 2003;Carretero-Paulet et al., 2010;Pires and Dolan, 2010). Nevertheless, the number of proposed subfamilies varies considerably among these studies. In the present one, the classification suggested by Carretero-Paulet et al. (2010) proposing the division of plant bHLH transcription factors into 32 subfamilies was used, since it represents the most recent and comprehensive study, so far.
From the BLASTP search at Phytozome, using all 17 Arabidopsis bHLH protein sequences from subgroup 25, 67 non-redundant homolog peptides were identified in the soybean genome. Seven of these were removed from the analysis as they did not contain any bHLH domain. Another 15 sequences were discarded after mismatch counting performed with their aligned domains. Using MEME, two other highly conserved motifs (with E-values of less than 1.7 e -851 ) were identified among the soybean subgroup 25 sequences. They are formed by residues right adjacent to the bHLH domain and had been previously reported (Heim et al., 2003;Li et al., 2006;Carretero-Paulet et al., 2010;Pires and Dolan, 2010). General characteristics related to the 45 remaining putative soybean bHLH genes are shown in Table 1. Remarkably, members of this subgroup were found spread throughout the 20 soybean chromosomes, with protein sequences ranging from 165 to 691 amino acids. Among the 45 annotated ORFs, 42 presented corresponding ESTs, suggesting that they are expressed genes and not pseudogenes. A complete overview of the gene expression results obtained for this group is presented in Figure 1. Differential expression in at least one of the stress situations/experiments available in LGE database was detected for 14 ORFs, four of these were differentially expressed in more than one situation and three respond to both abiotic and biotic stresses.
Lately, a growing number of studies accessing the functional role of specific plant bHLH transcription factors have been reported (Friedrichsen et al., 2002;Szécsi et al., 236 Osorio et al.  Liu et al., 2008;Chandler et al., 2009;Todd et al., 2010;Zheng et al., 2010). Nevertheless, a deeper (and broader) functional characterization of this family, focusing on the connection of members/subgroups to the biological processes they control, remains to be done. A first step in this direction has been recently taken by Carretero-Paulet et al. (2010) and Pires and Dolan (2010), where comprehensive information relating both classification and function of previously characterized plant bHLH transcription factors was assembled. More specifically, information regarding the function of subgroup 25 members is still scarce and concerns Arabidopsis members only. An alternative transcript of At1g59640 (ZCW32/BPE) seems to be involved in the control of petal size, whereas its counterpart is expressed ubiquitously (Szécsi et al., 2006). Furthermore, At4g34530 (CIB1) and At1g26260 (CIB5) were shown to interact with blue-light receptor CRY2 and promote floral initiation (Liu et al., 2008). Of most interest for this study, is the redundant role of At1g18400 (BEE1, Brassinosteroid Enhanced Expression1), At4g36540 (BEE2) and At1g73830 (BEE3) in brassinosteroids (BRs)/ABA antagonistic cross-talk during cell elongation (Friedrichsen et al., 2002). According to these authors, BEE1, 2 and 3 are early-response genes induced by BRs through the BRI1 receptor complex, and their expression is repressed by ABA through a yet unknown ABA receptor. Whether this pathway is also related to the ABA-dependent stress-responsive network, still requires further study. Moreover, Poppenberger et al. (2011) have demonstrated that At1g25330 (CESTA), a close homolog of BEE1 and BEE3 (Figure 2), is also involved in BR signaling, possibly by heterodimerization with its closest homologs. Remarkably, it has also been shown that lack of CESTA activity results in the misregulation of genes that are not only BR-responsive but also stress-responsive, such as Arabidopsis ERD5 (Early Responsive to Dehydration 5), TTL4 (Tetratricopetide-Repeat Thioredoxin-Like 4), WRKY18 and a putative LRR-disease resistance protein (Poppenberger et al., 2011), further suggesting that these pathways might indeed share common features.
As an attempt to predict gene function of the annotated genes, a comparison of their amino-acid sequences with subgroup 25 bHLH protein sequences from three other model plant species was carried out. Indeed, representative members from diverse taxonomic groups (P. patens, bryophytes; O. sativa, monocotyledonous; and A. thaliana, dicotyledonous) were included in the phylogenetic analysis in order to access the evolutionary features of this subgroup. The results obtained from the phylogenetic analysis proved to be consistent, since the clades formed were highly supported by posteriori probabilities (Figure 2, on left) and bootstrap (data not shown) analyses. Unlike previous phylogenetic reconstructions of the bHLH family that used the bHLH domain only, this study presents a tree reconstructed from full-length protein sequences. This adds accuracy and reliability to the tree resolution, since the short length of the bHLH domain (~60 amino-acids), along with its extremely high conservation within subgroups may compromise the reliability of the analysis (Amoutzias et al., 2004).
Patterns of intron distribution among bHLH-encoding genes from diverse species were shown to be conserved within subgroups and provide another criterion in phylogenetic analysis (Li et al., 2006;Carretero-Paulet et al., 2010). In this study, the overall intron-exon organization of bHLH subfamily 25-encoding genes from soybean and other three species was established (Figure 2, on right). Among 89 sequences, the number of introns ranged from 1 (Pp1s270_17v6) to up to 12 (LOC_Os03g12940), and in many cases, phylogenetically related proteins exhibited a closely related gene structure, corroborating the clustering results.
Since it is a basal species among land plants, the moss P. patens was added to this classification in order to help infer about this group's ancestral state (Rensing et al., 2008). Notably, all 12 members from P. patens grouped together into a clade, instead of grouping with the other plant species, indicating that the radiation within this subgroup has occurred independently in mosses and vascular plants, after the divergence of these taxonomic groups. The same result was obtained by Carretero-Paulet et al. (2010), even when a different method was applied [maximum likelihood (ML) analysis from bHLH-domain alignments]. Nevertheless, the chance that genes belonging to this subgroup might have independently evolved similar functions in both mosses and vascular plants should not be discarded, as suggested by Menand et al. (2007). In fact, while studying plant bHLH ancestry, Pires and Dolan (2010) concluded that the complex regulatory machinery that may be observed in modern plant lineages actually arose early in plant evolution.
The most striking feature that can be inferred from our phylogenetic analysis, which is in accordance with other previously published plant bHLH phylogenies (mentioned above), is the importance of gene duplication during the evolution of this family as a whole. Recurring events of single-gene duplications ("birth-and-death evolution"), combined with domain shuffling seem to rule bHLH evolution and diversification (Morgenstern and Atchley, 1999;Amoutzias et al., 2004;Nei and Rooney, 2005). Furthermore, whole genome duplication (WGD) events also seem to have had an active effect (as seen in the outer clades in Figure 2, on the left), and this seems to be even more intense in the soybean genome. According to our results, the subgroup in question encompasses almost three times more members in soybean than in Arabidopsis or rice (Table 1), with many of these grouping into new clades with no apparent near orthologs in the other analyzed species (Figure 2, in gray on the left side). Indeed, soybean suffered from two WGD events with an impressive retention of homologous blocks (Schmutz et al., 2010). Furthermore, specifically in the case of transcription factors (and other genes working in complex networks), duplications resulting from WGD events are vastly overretained, simply because they may be too costly to be removed, thus making functional redundancy a common feature among transcription factors, especially in plant species. Once retained, homologous duplicates might diverge in function or even subfunctionalize (Freeling, 2009), thus providing a source of evolutionary novelty in the form of new regulatory networks (Carretero-Paulet et al., 2010).
With all that in mind, an integrated analysis of both the expression profile ( Figure 1) and the phylogeny (Figure 2) presented herein provides a hint at the roles of subgroup 25 bHLH soybean genes. By focusing on soybeannear homologs shown in the tree (Figure 2 on left) we could see that for most of the paralogs whose expression has been detected, a divergent profile seems to prevail. An exception would be the cases of Glyma03g31510 and Glyma19g34360, which were both repressed during drought stress, with a broadly negative response in the latter, as its mRNA levels were down-regulated in both the susceptible and the tolerant cultivars analyzed. Moreover, the transcripts from Glyma19g32570 were up-regulated during ASR infection in the resistant genotype, whereas its counterpart Glyma03g29710 exhibited opposite differential expression. The near paralogs Glyma05g01590 and Glyma17g10290 also seem to be moving in different directions. Whereas the first seems to be up-regulated in response to fungal stress, the latter seems to be broadly down-regulated, in both susceptible and tolerant cultivars submitted to drought, as well as in P. pachyrhizi's infection. Furthermore, while Glyma15g33020 seems to be positively involved in soybean defense against ASR and during drought stress in tolerant Embrapa-48 cultivar, its nearest paralog (Glyma09g14380) was not differentially expressed in any of the situations assessed, and their near homolog Glyma17g08300 seems to be negatively involved in drought stress responses, since it was down-regulated in the same cultivar. Whether the examples mentioned above reflect functional divergence or subfunctionalization among duplicate homologs still requires further analysis.
Even though comparison of soybean genes with their orthologs in other species (such as Arabidopsis) is a tentative approach, and as such needs to be performed carefully. In this context it would be interesting to address the function of BEE orthologs in soybean, so as to determine whether they are similar to their Arabidopsis counterparts, and whether they somehow connected to stress responses. In this respect, special attention should be given to Glyma05g35060, which clustered together with the Arabidopsis BR-responsive genes, and whose transcripts turned out to be up-regulated in Embrapa-48 tolerant cultivar in response to drought. Previously reported bHLH genes were identified according to their accession/locus numbers, the other genes were designated according to their locus ID in Phytozome. A. thaliana (At); G. max (Glyma); O. sativa (LOC_Os) and P. patens (Pp). The graph on the right shows gene organization of full-length coding sequences from 89 plant bHLHs. Intron-exon maps were drawn using Fancy Gene v1.4 software, according to sequence data available in Phytozome.

Identification and analysis of soybean trihelix-GT encoding genes
The first isolated and described soybean GT-factor was GmGT-2 (Glyma02g09060), which binds to an element within the Aux28 promoter, and whose mRNA levels were down-regulated by light in a phytochrome-dependent manner (O'Grady et al., 2001). In a global approach using massive EST analysis, Tian et al. (2004) identified 13 putative trihelix genes in the soybean genome. Two of these [GmGT-2A (Glyma04g39400) and GmGT-2B (Glyma10g30300)] were cloned and had their roles in abio-tic stress tolerance described using transgenic Arabidopsis plants (Xie et al., 2009). The current annotation analysis indicates the occurrence of at least 63 GT-like genes in the soybean genome. 56 of these had their expression confirmed in the NCBI databases (Table 2). Unfortunately, since information available in Phytozome is not yet definitive, full-length cDNAs were not obtained for most sequences, so only gene-models were considered for this analysis. The 63 soybean trihelix-GT genes encode proteins with lengths ranging from 201 to 885 amino acids, distributed across most of the soybean chromosomes, ex- 240 Osorio et al. cept for chromosomes 5 and 14. There is an average of 3.5 GT-factor-encoding genes per chromosome, with the highest number of 9 genes found in chromosome 10, whereas a single member was detected in chromosomes 12 and 17, respectively. Three genes (Glyma09g19750, Glyma10g34610 and Glyma20g30630) with incorrect gene model predictions were manually curated.
Mining the LGE gene expression superSAGE experiments revealed that 11 soybean trihelix-GT genes were differentially expressed in the abiotic/biotic conditions tested (Figure 3). In accordance with our analyses, five trihelix-GT genes were up-regulated under drought in the tolerant cultivar (Embrapa-48), whereas only two genes were downregulated in this genotype. In the susceptible cultivar (BR16), Glyma10g34520 had its transcript levels increased in response to water deficit and the opposite situation occurred with Glyma10g36950. When plants were infected with P. pachyrhizi, only two genes displayed up-regulation of mRNA levels in response to biotic stress whereas two others seemed to be down-regulated. Interestingly, none of the soybean trihelix-GT previously reported as responsive to stress conditions and particularly to abiotic stress [GmGT-2A (Glyma04g39400) and GmGT-2B (Glyma10g30300)] were detected in the superSAGE experiments herein assessed. Divergence in experimental parameters and genotypes used might explain this unexpected result.
Transcript levels from Glyma01g35370 and Glyma20g30640 increased when plants were infected with ASR, while the opposite situation occurred with Glyma16g28240 and Glyma17g13780 mRNA levels. A rice GT-factor (OsRML1) was already reported to be upregulated in response to Magnaporthe grisea , which corroborates a connection between pathogen attack and trihelix-GT gene regulation. It is also possible that Glyma01g35370 may be involved in plant responses to both abiotic and biotic stresses, since the gene expression profile was modulated during water deficit and P. pachyrhizi infection.
The superSAGE experiments suggested that, at least in some cases, the same gene has variable transcript levels in different cultivars and/or in response to different stresses or agents. For example, when water deficit was imposed on soybean plants, Glyma10g36950 was down-regulated in the susceptible (BR16) and the tolerant (Embrapa-48) cultivars, whereas its transcript levels did not change in response to ASR. In another case, Glyma09g38050 was upregulated in response to drought stress in Embrapa-48, but no differences were detected in BR16. Furthermore, Glyma13g26550 was down-regulated in response to drought stress in the tolerant cultivar, whereas its expression in cultivar BR16 did not exhibit any alterations. In these cases, in addition to differential gene regulation, there may be other factors contributing to distinct regulatory function, such as post-translational modifications or variation in dimerization partners (Zhou, 1999). Modifications in individual cis-regulatory elements on trihelix-GT promoter regions of duplicated genes might lead to the processes of transcriptional neofunctionalization or subfunctionalization (Haberer et al., 2004), which may explain gene induction or repression without any counterpart response during the same stimuli. This seems to be the case for Glyma03g07590 and its nearest paralog Glyma01g29760, or for Glyma16g28240 and the phylogenetically related Glyma02g09050. Further studies focusing on identifying cis-elements, as well as performing promoter analyses to verify inducible expression patterns may clarify the involvement of duplicated genes in stressrelated responses.
A previous study regarding the phylogenetic analysis encompassing Arabidopsis and rice GT factors (Fang et al., 2010) showed that this family could be classified into three subfamilies (a , b and g), with unique composition of predicted motifs. Unfortunately, these results were not reproduced in our analysis, even when full-length protein sequences (Figure 4) or the trihelix domains alone were aligned (data not shown). An exception occurred with subfamily g, which had already been described as having low sequence similarity with the other reported GT factors. The introduction of soybean and M. truncatula sequences in the phylogeny might have affected the expected distribution within those subgroups. Besides, we also inserted into our tree the soybean gene AAK69274 described by Fang et al. (2010), which could neither be identified in the soybean genome nor detected in the expression database. According to our analysis, this unexpected result seems to indicate the 242 Osorio et al.  occurrence of an alternative splicing in Glyma19g37410 or Glyma03g34730, both considered to be phylogenetically closest to the unidentified gene locus. Hence, when taking into account the full-length protein sequence, the GT-factor family might be divided into two subgroups, in one of these subgroups a branch corresponded to the already described subfamily g (Figure 4, in gray). Despite the fact that subfamilies a and b were not distinguished, other probabilities supported our tree, especially when inner nodes were observed.
When gene organization among Arabidopsis and soybean sequences was compared (Figure 5), the number of introns ranged from zero (twenty three genes) up to 16 (At5g63420 and Glyma06g17980), and some phylogenetically close sequences showed the same gene structure. For example, the Arabidopsis At3g10040 and its soybean ortholog do not have intron, whereas At2g33550 and related members have two introns, with remarkable differences in intron size.
As observed for bHLH transcription factors, the soybean GT factor family encompasses almost three times more members than Arabidopsis or rice, a consequence of the WGD events that took place during plant evolution. In several cases, soybean paralogs clustered with one M. truncatula gene, indicating that these paralogs probably derived from a WGD event that occurred after the divergence of the two legume species. Similarly, Schmutz et al. (2010) refer to a Glycine-specific WGD event, estimated to have occurred about 13 million years ago. However, the possibility that extra M. truncatula orthologs might arise upon the completion of its genome sequencing should not be discarded.
Recently, the OsGTg subfamily was proposed to participate in the regulation of stress tolerance in rice (Fang et al., 2010). OsGTg -1 showed more specific expression pattern than their counterparts OsGTg-2 and OsGTg-3, which are supposedly redundant. None of them was responsive to light, but their transcript levels increased in response to salt and cold stresses, whereas OsGTg-1 was upregulated by ABA and SA stimulus. It is possible that some soybean members of this subfamily may act in response to stressor agents, but more studies are required in order to understand whether the pattern seen in rice GTg factors also occurs in soybean and M. truncatula. Our analysis, so far, does not indicate their involvement in an abiotic and/or biotic stress response. Moreover, soybean genes previously reported as involved in stress responses (Xie et al., 2009) together with other genes herein identified are dispersed in different tree branches, indicating that this family is in fact evolutionarily diversified.

Conclusion
The present study identified new members of soybean bHLH and trihelix-GT transcription factor families, some of which seem to be involved in responses to environmental stresses. It also emphasizes the role of duplication events in the expansion and evolution of soybean transcription factor families, indicating that exciting new layers of complexity might exist in this species' regulatory mechanisms, including biotic and abiotic stress responses.