Phosphorus stress in common bean: root transcript and metabolic responses.

Phosphorus (P) is an essential element for plant growth. Crop production of common bean (Phaseolus vulgaris), the most important legume for human consumption, is often limited by low P in the soil. Functional genomics were used to investigate global gene expression and metabolic responses of bean plants grown under P-deficient and P-sufficient conditions. P-deficient plants showed enhanced root to shoot ratio accompanied by reduced leaf area and net photosynthesis rates. Transcript profiling was performed through hybridization of nylon filter arrays spotted with cDNAs of 2,212 unigenes from a P deficiency root cDNA library. A total of 126 genes, representing different functional categories, showed significant differential expression in response to P: 62% of these were induced in P-deficient roots. A set of 372 bean transcription factor (TF) genes, coding for proteins with Inter-Pro domains characteristic or diagnostic for TF, were identified from The Institute of Genomic Research/Dana Farber Cancer Institute Common Bean Gene Index. Using real-time reverse transcription-polymerase chain reaction analysis, 17 TF genes were differentially expressed in P-deficient roots; four TF genes, including MYB TFs, were induced. Nonbiased metabolite profiling was used to assess the degree to which changes in gene expression in P-deficient roots affect overall metabolism. Stress-related metabolites such as polyols accumulated in P-deficient roots as well as sugars, which are known to be essential for P stress gene induction. Candidate genes have been identified that may contribute to root adaptation to P deficiency and be useful for improvement of common bean.

Common beans (Phaseolus vulgaris) are the world's most important grain legume for direct human consumption; they comprise 50% of the grain legumes consumed worldwide (Broughton et al., 2003;Graham et al., 2003). In several countries of Central and South America, beans are staple crops serving as the primary source of protein in the diet. Environmental factors, such as low soil nitrogen (N) and phosphorus (P) levels, and acid soil conditions are important constraints for bean production in most of the areas where this crop is grown (Graham et al., 2003). In bean, symbiotic N fixation rates, seed protein level, and tolerance to P deficiency are low in comparison to other legumes (Broughton et al., 2003). P is an essential element required for plant growth and development. Besides N, P is the most limiting nutrient for plant growth, and it is a common limiting factor for crop production in arable soils. Plants have evolved general strategies for P acquisition and use in limiting environments that include: mycorrhizal symbioses, decreased growth rate, remobilization of internal inorganic phosphate (P i ), modification of carbon (C) metabolism bypassing P-requiring steps, increased production and secretion of phosphatases, exudation of organic acids, modification of root architecture, expansion of root surface area, and enhanced expression of P i transporters (for review, see Raghothama, 1999;Smith, 2001;Vance et al., 2003;Plaxton, 2004).
In contrast to disease-resistance traits, where resistance may be due to a single dominant or recessive gene, enhancing tolerance to P stress requires multiple genes and involves several different mechanisms. In recent years, macro/microarray technologies have provided valuable information on global changes in gene expression in response to P starvation in several plant species and organs, including white lupin (Lupinus albus) proteoid roots , rice (Oryza sativa) leaves and roots (Wasaki et al., 2003(Wasaki et al., , 2006, and Arabidopsis (Arabidopsis thaliana) roots, shoots, and leaves (Hammond et al., 2003;Wu et al., 2003;Misson et al., 2005;Mü ller et al., 2007).
Although macro/microarray studies have identified genes differentially regulated by P starvation, little is known about the regulation of gene expression changes. Transcription factors (TFs) are master control proteins in all living cells, regulating gene expression in response to different stimuli (Riechmann, 2002;Czechowski et al., 2004). Chen et al. (2002) reported that Arabidopsis TF gene expression is regulated in a cell type-or tissue-specific manner and in response to specific environmental biotic and abiotic stresses. Mü ller et al. (2007) reported that specific TFs are induced in Arabidopsis P-starved leaves. These studies have opened new possibilities to elucidate the sensing, signaling, and regulatory pathways of the P deficiency response in plants.
Despite the agronomic importance of beans, there is little information on global gene expression of bean tissues in response to P deficiency. In previous work, we attempted to identify candidate P stress-induced genes in beans using an in silico approach that clustered bean ESTs with previously identified P stressinduced genes across three other legume species and Arabidopsis (Graham et al., 2006). Here, we undertook a three-step approach to identify genes important to P deficiency in common bean. First, macroarray technology was used for transcript profiling of P-deficient bean roots with the aim of identifying those genes, gene networks, and signaling pathways that are important for the plant response to P deficiency. Second, we identified bean TFs and used quantitative reverse transcription (RT)-PCR to assess TF gene expression in P-deficient bean roots, with the aim of identifying TFs that regulate the differential expression of genes during P stress. Third, we performed nonbiased metabolite profiling of bean roots using gas chromatography coupled to mass spectrometry (GC-MS) to correlate metabolic differentiation orchestrated by global changes in gene transcription as response to P starvation. The overall goal of this research is to identify candidate genes that may be useful to bean improvement and that will contribute to understanding common bean adaptation to P deficiency.

Phenotypic Characterization
The long-term P deficiency treatment used in this work consisted of growing common bean plants in pots under controlled environments for 3 weeks using 200-fold lower phosphate concentration as compared to P-sufficient (1P) control plants. Control plants accumulated higher concentrations of soluble P i . P i content in 1P leaves was 2.6-and 13-fold higher than in 1P stems and 1P roots, respectively (Fig. 1A). Compared to 1P plants, a drastic reduction (2-23-fold lower) in P i content was observed in plants grown under P-deficient conditions (Fig. 1A). P i content in P-deficient plants was similar in leaf, stem, and root tissues (Fig. 1A). Typical P stress responses were observed (Raghothama, 1999;Gilbert et al., 2000;Ma et al., 2003), including a 4-fold reduction in leaf area and 1.5-fold higher dry weight root to shoot ratio (Fig. 1,B and C). The latter response was due to arrested shoot growth and proliferation of lateral roots and root hairs of P-deficient plants.
Content of photosynthetic pigments such as chlorophyll a and b and carotenes was similar in plants under 2P and 1P treatments (data not shown). However, P-deficient plants showed significant inhibition of net photosynthetic rate (P n ) regardless of internal CO 2 (C i ) concentration (Fig. 1D). In contrast, P-deficient plants showed 50% lower P n at ambient CO 2 concentration (350 mmol mol 21 ), reflecting lower carboxylation efficiency. In addition, P-stressed plants showed 60% of the maximum P n of 1P plants, which is consistent with changes associated with increasingly larger limitations of P n by Rubisco and ribulose 1,5-bisphosphate regeneration as leaf P i declines (Fig. 1D). However, stomatal conductance and resistance was not altered in P-deficient plants (data not shown).

Macroarray Analysis of Root Response to P Deficiency
Macroarray analyses were performed to evaluate gene expression from P-deficient roots of bean plants as compared to control P-sufficient roots. Nylon filter arrays were spotted with ESTs that represented a 2,212 bean unigene set consisting of 1,194 singletons and 1,018 contigs derived from the 2P roots cDNA library from bean 'Negro Jamapa 81' previously reported (Ramírez et al., 2005;Graham et al., 2006).
Total RNA was isolated from plants grown under similar conditions as described for each treatment (2P and 1P). Ten nylon filter arrays were hybridized with first-strand cDNA synthesized from four independent sources of total RNA. From the 10 hybridizations, six replicates with high determination coefficients (r 2 $ 0.8) were chosen for analysis of differential gene expression. A total of 126 cDNAs showed significant (P # 0.05) differential expression (Tables I and II).
Tables I and II list the genes that were significantly induced or repressed, respectively, in P-deficient roots. To aid in annotation, cDNAs were assigned to tentative consensus sequences (TCs; Institute of Genomic Research [TIGR]/Dana Farber Cancer Institute [DFCI] Common Bean Gene Index, v. 1.0) when possible. The TC or EST sequences were then compared (BLASTX, E , 10-4;Altschul et al., 1997) to the Uniprot protein database (Apweiler et al., 2004) to assign putative function. Based on information available in the literature, sequences were then assigned to functional categories. Table I shows the genes (78) that were induced 2-fold or more in P-deficient roots, classified in nine functional categories. The ''unknown function'' category included those genes with similarity to hypothetical proteins with unknown function and those for which no BLAST hit was found. The two most abundant functional categories, accounting for 23% of genes each, were the regulation/signal transduction category and those coding for genes that participate in secondary metabolism pathways and/or are related to several stress/defense plant responses. Ten genes (13%) were classified as membrane proteins or proteins that participate in transport, both extracellular and intracellular. Six genes (8%) were classified in cell structure, cell cycle, or developmental functions. Nineteen genes (24%) were classified in different metabolic pathways: P i cycling, C and N metabolism, amino acid/protein synthesis or degradation, and lipid metabolism. Finally, 9% of genes had no known function. Table II lists the functional classification of the genes (48) that were repressed in 2P roots as compared to control roots. The most abundant category was the amino acid/protein metabolism with 11 genes (23%). Only five genes participating in metabolic C/N pathways were identified (10%), and no genes involved in P i cycling were identified. Nine (19%) and seven (15%) genes were classified in the transport/membrane protein and cell structure/cell cycle/development categories, respectively. Only 8% and 6% of the repressed genes participate in regulation/signal transduction and secondary metabolism/defense pathways, respectively.
It was evident that a number of genes from within a single functional category could either be induced (Table I) or repressed (Table II). We found that 10 P deficiency-induced genes identified by the macroarray analysis had been previously proposed by Graham et al. (2006) as candidate P stress-induced genes in bean (Table I). Graham et al. (2006) identified candidate P stress-induced genes of bean by statistical analysis of contigs overrepresented with ESTs from P-stressed tissues and by clustering candidates with P stress-induced genes identified from a variety of plant species, including Arabidopsis, lupin, soybean, Medicago truncatula, and bean. As expected, none of the 2P-repressed genes identified by macroarrays ( Table   Figure 1. Effect of P deficiency on common bean. A, Soluble P i content in different plant organs. B, Leaf area from fully expanded leaves. C, Root to shoot dry weight ratio. D, P n rate as a function of changing C i . Plants were grown for 3 weeks under P-deficient (black bars or circles) or in P-sufficient conditions (white bars or circles). Values are mean 6 SE from 12 determinations: three independent experiments with four replicates per experiment. II) were included in the Graham et al. (2006) analysis, which only evaluated induced genes.

Expression Analyses of Selected Genes by RT-PCR
Nine ESTs selected from both Tables I and II (18 total) were chosen to assess whether macroarray expression data could be confirmed by an alternate method. We performed semiquantitative RT-PCR on ESTs representing at least four functional categories designated in Tables I and II. As shown in Figure 2, all 18 genes tested for expression by RT-PCR gave results confirming their expression obtained with macroarray experiments. From the P deficiency stress-induced genes, UDP-Glc-6-dehydrogenase, senescence-related dihydroorotate dehydrogenase, glycolipid transfer protein, and hypothetical protein were the most highly induced genes in their particular categories, as measured by macroarrays. These genes showed enhanced expression by RT-PCR ( Fig. 2A). Likewise, from the P deficiency-repressed genes in Table II, isocitrate dehydrogenase, SAM-decarboxylase, multidrug resistance protein, and caffeine-induced death protein were among the most highly repressed genes detected by macroarray analysis (Fig. 2B), and these genes showed reduced expression in P deficient as compared to P sufficient when evaluated by RT-PCR.

TF Transcript Profiling by Real-Time RT-PCR
The TIGR/DFCI Common Bean Gene Index contains 9,484 total unigenes (2,906 TCs and 6,578 singletons) comprised of 21,290 input EST sequences. The first step in our work was to define the set of bean EST/TC sequences in the TIGR/DFCI Common Bean Gene Index (www.tigr.org; http://compbio.dfci.harvard. edu/tgi/plant.html) coding for proteins with Inter-Pro domains diagnostic or characteristic of TF genes. A total of 372 sequences, corresponding to 4% of the bean unigene set, was identified using 41 of the preselected Genes reported as bean candidate P stress-induced genes through clustering analysis across five or four plant species by Graham et al. (2006). c Annotation according to TF genes identified in this work (Table III; supplemental data). TF diagnostic Inter-Pro domains. This constitutes the whole set of TF genes used for our real-time RT-PCR analyses. Most likely, some of the genes are not true TFs; however, they were included because they contain DNA-binding and other domains that are characteristic of TF proteins. Based on the classification of Arabidopsis TF gene families (Riechmann, 2002; http://range.gsc.riken.ip/rart; http://daft.cbi.pku.edu. cn), bean TF genes were grouped into 47 families (Fig. 3).
Although TF classes in bean were restricted to those identified from cDNA libraries, a general correspondence was found between the most abundant TF families in beans and those from Arabidopsis (Riechmann, 2002), such as the MYB superfamily with 46 gene members (12%), C2H2(Zn) (10%), and AP2/EREBP (8%; Fig. 3). However, in our dataset, we found that CCAAT and bHLH families were equally abundant in our bean TF gene set (Fig. 3), while in Arabidopsis the bHLH Figure 2. Verification of macroarray results by RT-PCR analysis. Selected genes identified as induced (A) or repressed (B) in P-deficient roots were evaluated. The ubiquitin gene was included as control for uniform RT-PCR conditions (bottom). The primer sequences and reaction conditions used are presented in Table V.  (372) were grouped in 47 different families with different Inter-Pro domains according to TF gene families reported for Arabidopsis (Riechmann, 2002; http://range.gsc.riken.ip/rart; http://daft.cbi.pku.edu. cn). The identity of each TF gene family with three or more members is shown. Twelve gene families with two members each (2/TF fam) are: TAZ, MBF1, ARID, Nin-like, Dof-type C2C2(Zn), S1Falike, YABBY C2C2(Zn), BES1, K-box, Histonelike/CBFA_NFYB_topo, Auxin_resp, and Lambda_ DNA_bd. Eleven gene families with one member each (1/TF fam) are: FHA, LIM-domain, E2F/DP, Jumonji JmjN, SBP, SHAQKYF_MYB_bd, ZF_HD, SRS, POX, EIL, and Euk_TF_DNA_bd. family is around 3-fold more abundant than the CCAAT family (Riechmann, 2002). Other families of bean TF genes consisted of between one and 12 genes (Fig. 3).
We performed TF profiling based on real-time RT-PCR to determine differential expression of bean TF genes that might be involved in gene expression response to P deficiency. There were three biological replicates of 2Pand 1P-treated roots. In each RT-PCR run, the phosphatase gene (TC201) was included as a P-deficient marker. This marker gene, known to be induced in P-deficient roots (Ramírez et al., 2005), showed an average expression ratio 2P to 1P of 18.48 (P 5 0.005), confirming the P-deficient status of the roots. From the 372 TF genes, 46 (12%) were differentially expressed (P # 0.05) in 2P-treated roots, 10 were induced, and the rest were repressed in 2P roots. Table  III shows those TF genes that were induced (four) or repressed (13) 2-fold or more in P-deficient roots. To annotate the P-regulated TFs, the TC sequences were blasted (BLASTX, E , 10 24 ; Altschul et al., 1997) against the Uniprot protein database (Apweiler et al., 2004; Table III).
Most of the TF genes induced in 2P roots belong to the MYB superfamily (Table III). The induction of Arabidopsis MYB TF genes in response to different biotic stresses (Chen et al., 2002) and to P starvation (Mü ller et al., 2007) has been shown previously. It has been demonstrated that the Arabidopsis PHR1 and PHR2 genes, which resemble the PSR1 gene from Chlamydomonas reinhardtii and belong to the TF MYB superfamily, are crucial for P starvation signaling Todd et al., 2004). Our BLAST analysis revealed that the deduced translated amino acid sequence of MYB TF TC2883, induced 2-fold in 2P roots (Table III), showed 59% amino acid identity to PHR1 (BLASTX E value 5 4.1E 239 ). The C2C2(Zn) TF family was the most highly represented among the repressed TF genes, and members from eight other TF gene families were also repressed (2-fold or more) in 2P roots (Table III).

Metabolome Analyses
To assess the degree to which changes in plant gene expression in P-deficient bean roots affect overall metabolism, we performed nonbiased metabolite profiling of bean roots using GC-MS. The complete information of the 81 metabolites and mass spectral metabolite tags (MSTs) detected in bean roots subjected to both treatments (2P and 1P) is provided as supplemental data. Table IV shows the retention time index (RI) value and RI SD of those metabolites and MSTs (42) with 2P to 1P response ratios 1.5-fold or more and those with lower ratios but highly significant (P # 0.05). The metabolites thus identified were in agreement with previous analyses (Desbrosses et al., 2005), mostly primary metabolites belonging to the compound classes: amino acids, organic acids, polyhydroxy acids, fatty acids, sugars, sugar phosphates, polyols, and other nitrogenous compounds. Most of the metabolites showed a response ratio higher than 1, indicating an increase in P-deficient roots; only eight metabolites were decreased in P-stressed roots (Table IV). Most of the amino acids were increased in P-stressed roots; in addition, the polyols and sugars showed high and significantly different 2P to 1P response ratios (Table IV).
Quantitative data for the metabolites listed in Table  IV were used for independent component analysis a Whenever the ratio was lower than 1 (genes repressed in P-deficiency), the inverse of the ratio was estimated and the sign was changed. (ICA) to identify major differences in metabolite composition in P-deficient and normal roots. ICA of metabolite response ratios of all 81 metabolites in 12 samples from P-deficient roots and 12 samples of P-sufficient roots allowed nonbiased partitioning into two sample groups showing gradual differentiation of individual plants from a P-sufficient metabolite phenotype (Fig. 4). The score plots (Fig. 4) show a clear separation between 2P and 1P samples, though some overlap in the samples can be seen, which probably indicates a P deficiency but not total P starvation in bean roots. The response ratio of average 2P root response compared with average 1P root response is listed (t test significance of P # 0.05 is indicated by bold format of the response ratio). For ratios lower than 1, the inverse of the ratio was estimated and the sign was changed. b Represents the sum of two or more metabolites.

DISCUSSION
In this report, we have advanced the fundamental understanding of common bean root gene expression and plant adaptation to P deficiency by: (1) identifying differential patterns of gene expression in P-stressed roots through macroarray analysis; (2) identifying 372 TFs and evaluating their expression profile by quantitative RT-PCR; and (3) complementing gene expression analysis with unbiased metabolomic profiling. Transcript expression patterns revealed by hybridization of nylon filter arrays spotted with some 4,000 ESTs from bean 2P roots cDNA library (Ramírez et al., 2005) resulted in a total of 126 differentially expressed genes with 2-fold or more induction or repression in 2P roots (Tables I and II). In addition, transcript profiling of 372 TF genes identified from the bean gene index (TIGR/DFCI) resulted in 17 differentially expressed (2-fold or more) bean TF genes in 2P versus 1P bean roots (Table III). Nonbiased metabolite profiling using GC-MS technology led to the identification of 64 metabolites and 17 MSTs from bean roots, 42 of which showed $1.5-fold and/or significantly different 2P to 1P response ratios (Table IV). ICA analysis from the 81 identified metabolites revealed a gradual differentiation of individual plants from a P-sufficient metabolic phenotype (Fig. 4). Our results reveal a suite of responses ranging from changes in growth and development to altered gene expression and metabolic profile that may contribute to adaptation of common bean roots to P deficiency.
An overriding question regarding our macroarray experiments is: are genes designated as having enhanced expression during P stress in actuality responding to low P, or do they show enhanced expression due to root developmental effects? Several pieces of evidence suggest that a great many bean genes are responding to P stress. First, of the 50 TCs listed as induced during P stress in Table I, more than 80% have the majority of ESTs derived from a P stress root library. In fact, 11 of the 50 have 100% of their ESTs derived from the P stress root library. Second, semiquantitative RT-PCR of several P stress-induced genes (Fig. 2) show enhanced expression in P-stressed roots. Furthermore, an in silico statistical analysis of ESTs overrepresented in P stress libraries in legumes and Arabidopsis gene indices, similar to that described by Graham et al. (2006), showed that at least 50% of the TCs in Table I would be predicted to be highly expressed under P stress. Unfortunately, similar statistical methods cannot be applied to singletons or to underexpressed (underrepresented EST) TCs. However, semiquantitative RT-PCR of a number of underexpressed genes in Table II showed that they had reduced expression in P-stressed roots as compared to P-sufficient roots (Fig. 2).
As an initial step in responding to P deficiency, plants must sense that nutrient stress is occurring and transduce appropriate signals into processes that facilitate adaptation. Although the genes affecting P stress signal recognition and transduction in legumes are unknown, studies in Arabidopsis and rice have implicated MYB (PHR1), WRKY (WRKY75), and bHLH (OsPTF1) TFs in the P-signaling process Yi et al., 2005;Devaiah et al., 2007). Recently, the interaction of miRNA 399 with ubiquitin-conjugating enzyme (UBC) has also been demonstrated to play a key role in the P stress response of Arabidopsis (Fujii et al., 2005;Miura et al., 2005;Chiou et al., 2006). Our array study as well as those of others (Hammond et al., 2003;Uhde-Stone et al., 2003;Wu et al., 2003;Misson et al., 2005;Morcuende et al., 2007;Mü ller et al., 2007) have detected a plethora of putative signaling and regulatory genes that could be involved in P stress signaling. We found some 39 genes (Tables I-III) that may contribute to P stress signal transduction/regulation in common bean roots. As in Arabidopsis, we found representatives of MYB, UBC, and bHLH gene families as either up-or down-regulated in expression. We detected three members of the MYB superfamily that were induced in P-deficient roots (Table III). Of these, TC 2883 had the highest similarity (93%) to Arabidopsis PHR1, a MYB gene implicated in the P deficiency response process. Three Arabidopsis genes have recently been documented to be involved in signal transduction and regulation of P acquisition/homeostasis. These genes encode WRKY75, PHO2 (an E2 conjugase), and SIZ1 (a SUMO E3 ligase; Miura et al., 2005;Aung et al., 2006;Bari et al., 2006;Devaiah et al., 2007). Common bean TCs 2419, 1095, and 2445 have high similarity to WRKY75, PHO2, and SIZ1, respectively. Although none of the bean TCs cited above had enhanced expression in P-stressed roots of common bean, their similarity of Arabidopsis P-signaling genes suggests that comparable bean TCs may play similar roles in bean. Noticeable is that all of the ESTs comprising TC 2419 and 1095 are derived from the P-stressed root library. Interestingly, TC 1622, which is up-regulated in P-deficient bean roots, encodes a putative UBCligase related to a pepper CaPUB1 that has been implicated in resistance to abiotic stress (Cho et al., 2006).
Studies with white lupin Liu et al., 2005) and Arabidopsis (Nacry et al., 2005;Karthikeyan et al., 2006;Mü ller et al., 2007) have shown that sugars and P stress signaling are closely interrelated. Rychter and Randall (1994) found that within 2 weeks of P stress, common bean partitioned more sugars to roots than nonstressed plants. The enhanced expression of P stress-induced genes requires the presence of available sugars. Deprivation of sugars by either shading or stem girdling blocks the expression of P stress-induced genes (Liu et al., 2005). Our metabolic analysis of bean P-stressed roots provides additional support for the role of sugars in P stress. Several sugars (Table IV) were more abundant in P stress roots as compared to P-sufficient roots, suggesting that sugars may be partitioned preferentially to P-stressed roots to support the expression of P stress-induced genes. It is noteworthy that PRL1-associated protein, encoded by CV543658, which has enhanced expression in P stress bean roots, is known to interact with SNF1 to derepress Glc metabolism, stimulate starch accumulation, and inhibit root elongation (Bhalerao et al., 1999). These traits are characteristic of P-stressed plants.
It is also worthwhile to note the reduced amounts of organic acids in P-stressed roots as compared to P-sufficient roots (Table IV). It is well known that P-stressed legume roots release organic acids as a P-adaptive mechanism (Johnson et al., 1996;Neumann and Römheld, 1999;Shen et al., 2002;Dong et al., 2004). Release of organic acids into the rhizosphere enhances P i solubilization, making P more available. The reduced amount of organic acids in P-deficient roots more than likely reflects exudation from the root into the rhizosphere. The altered organic acid content of P-stressed roots is also reflected in the reduced expression of TC 1864 isocitrate dehydrogenase-ICD (Table II). This enzyme is a key regulatory enzyme in the tricarboxylic acid cycle. Reduced expression of ICD would lead to a buildup of malate acids that could be available for exudation into the rhizosphere. Almost 23% of the genes showing enhanced expression in P-stressed bean roots encode proteins having roles in either stress/defense or secondary metabolism (Table I). Hammond et al. (2003) have shown that many genes that respond to P stress in Arabidopsis shoots also respond to other environmental challenges, including salinity, wounding, pathogen attack, anoxia, and other nutrient stresses. In bean, P stress results in the induction of oxidative responses, including increased lipid peroxidation, elevated peroxide levels, and increased catalase and peroxidase activity (Juszczuk et al., 2001). Our results confirm and extend this observation by showing that genes encoding proteins in several aspects of oxidative stress have enhanced expression in P-stressed roots. Moreover, several genes implicated in plant response to diseases, such as PR and NBS-LRR resistance, are up-regulated in bean P-stressed roots along with genes involved in phenylpropanoid synthesis (Table I). Similar patterns of gene activation have been noted for plants undergoing potassium, zinc, iron, and N deficiency stress Armengaud et al., 2004;Shin et al., 2005;van de Mortel et al., 2006).
Because enhanced P i transporter gene expression is highly indicative of the P i stress response (Raghothama, 1999;Smith, 2001), it was surprising that we did not find any P i transporter to be highly expressed in P-stressed common bean roots. In fact, we found only a single P i transporter EST derived from the P-stressed root library. The lack of P i transporter ESTs in the root library could reflect that the library was made from roots of 21-d-old P-stressed plants. Perhaps earlier sampling dates would have yielded more P i transporters. On the other hand, we did detect enhanced expression of other types of transporters, including a putative aquaporin, an ATP-binding cassette transporter, and an acetylglucosamine transporter (Table I).
As demonstrated in Figure 1C and previously shown in numerous studies, the root to shoot ratio increased in P-stressed plants as compared to P-sufficient plants.
The ratio change was due in part to proliferation of lateral roots in P-stressed plants. Modified root architecture in response to P stress has been noted previously in common bean (Rychter and Randall, 1994;Lynch, 1995;Ge et al., 2000;Liao et al., 2001;Lynch and Brown, 2001) and Arabidopsis (López-Bucio et al., 2003;Ma et al., 2003;Wu et al., 2003). Phosphate starvation was recently shown to induce determinant root development programs in Arabidopsis (Sánchez-Calderó n et al., 2005). Recently, quantitative trait loci for root architecture traits that correlate with P acquisition have been identified in bean, strengthening the importance of root structure for low P adaptation (Beebe et al., 2006). Modification of root architecture in response to P deficiency results from the interplay between internal balance of the phytohormones auxin, cytokinin, and ethylene (Gilbert et al., 2000;Williamson et al., 2001;Al-Ghazi et al., 2003;Ló pez-Bucio et al., 2003;Ma et al., 2003;Karthikeyan et al., 2006). As one might expect, we found several genes in bean roots related to phytohormone biosynthesis and activity to be responsive to P. Accompanying increased lateral root growth, genes involved in cell wall synthesis and growth were responsive to P.
Reduced shoot growth accompanied by reduced photosynthetic rate (Fig. 1) was symptomatic of P stress in bean. Phosphate content and photosynthesis are related in several ways, and alteration of photosynthesis as a result of P starvation has been shown for several plant species, including common bean (Rychter and Randall, 1994;Mikulska et al., 1998). It has been shown that tobacco plants grown under P deficiency have reduced photosynthate demand in sink organs, resulting in carbohydrate accumulation and decrease in net photosynthesis (Pieters et al., 2001). Our data support the proposition of Morcuende et al. (2007) that repression of photosynthesis may be a secondary response linked to lower demand of photosynthate and higher sugar levels during P limitation.
The results from this work provide an abundance of candidate genes with diverse function that are postulated to play important roles in adaptation of common bean plants to P deficiency. These newly identified genes may be of utility in marker-assisted selection for P-efficient genotypes. The identified candidate genes expand the current information available on the regulation and signaling pathways during P deficiency in plants. In future studies, we propose to define the precise roles of selected candidate genes using reverse genetics approaches.

Plant Material and Growth Conditions
The common bean (Phaseolus vulgaris) Mesoamerican 'Negro Jamapa 81' was used in this study. Plants were grown in controlled-environment (26°C-28°C, 16-h photoperiod) greenhouses at Centro de Ciencias Genómicas/ Universidad Nacional Autó noma de México (Cuernavaca, México) and Max Planck Institute of Plant Molecular Physiology (Golm, Germany), or in growth chambers at the University of Minnesota (St. Paul). Surface-sterilized seeds were germinated at 30°C for 3 d over sterile wet filter paper and then planted in pots with vermiculite or coarse quartz sand. Pots were watered 3 d per week with the plant nutrient solution reported by Summerfield et al. (1977). For 2P conditions, K 2 HPO 4 concentration of the plant nutrient solution was reduced from 1 mM to 5 mM. In 2P conditions, cotyledons from each plant were cut 1 week after planting. Plants were grown for 3 weeks before harvesting. Roots for RNA isolation were harvested directly into liquid nitrogen and stored at 280°C.

Soluble P i Concentration
Soluble P i content was determined in leaves, stems, and roots from plants grown for 3 weeks in 2P or 1P conditions using the colorimetric assay reported by Taussky and Shorr (1953). For each assay, tissues were harvested, weighed, and immediately homogenized in 10 N TCA. For each determination, 12 replicates were analyzed. These were derived from three independent experiments with plants grown in similar conditions with four replicate assays from each treatment (2P roots or 1P roots) per experiment.

Photosynthesis and Photosynthetic Pigments Content
The relationship between CO 2 assimilation rate (P n ), increasing C i , and stomatal conductance and resistance was determined using a portable photosynthesis system (LI-6200 Primer; LI-COR) in 2Pversus 1P-treated plants. The measurements from mature bean trifolia were undertaken in a greenhouse maintaining leaf temperature and photosynthetically active photon flux density at 25°C and 1,600 mmol m 22 s 21 , respectively. Each point represents the average of 12 determinations from three independent experiments with plants grown in similar conditions and four replicate assays from each treatment (2P roots or 1P roots) per experiment. The CO 2 assimilation rate was adjusted to each leaf area value.

EST Sequencing and Annotation
Because the macroarrays used in this study were spotted prior to sequencing, 65 of the spotted clones had poor quality sequence and were not included in sequence-based analyses (Ramírez et al., 2005;TIGR/DFCI, Quackenbush et al., 2001) or submitted to GenBank. To include these clones in our analyses, the clones were resequenced. DNA sequencing was performed at the Advanced Genetic Analysis Center (University of Minnesota) and at the Center for Genomic Sciences/Universidad Nacional Autó noma de México (Cuernavaca, México). The new sequences were submitted to GenBank (accession nos. EH791054-EH791109, EH792671-EH792678, and EH795233).
To assign newly sequenced ESTs to existing TCs in the TIGR/DFCI Common Bean Gene Index, the EST sequences were compared to the TCs using TBLASTX (Altschul et al., 1997). To confirm the placement of the EST with the putative matching TIGR/DFCI TC, the overlap of both sequences was checked using the SeqManII program in the DNASTAR software package. Sequence matching indicated that the analyzed bean sequence belonged to the TC, as indicated in Tables I and II. To annotate the sequences described in Tables I to  III, all sequences were cross referenced with the TIGR/DFCI Common Bean Gene Index to find the corresponding TCs or singletons. TC or singleton sequences were compared to the Uniprot (version July 2006;Apweiler et al., 2004) protein database using BLASTX and an E-value cutoff of E , 10 24 .

Nylon Filter Arrays, Hybridization, and Data Analysis
The preparation of a cDNA library from roots from P-deficient bean 'Negro Jamapa 81' plants and the sequence of ESTs (4,329) have been reported (Ramírez et al., 2005;Graham et al., 2006). For macroarray preparation, the cDNA portion of each root EST was amplified by PCR using standard T3 and T7 primers. The PCR products were spotted onto Gene Screen Plus membranes (NEN Life Science Products) using the Q-bot (Genetix) automated spotting system.
Total RNA was isolated from 4 g frozen roots (as reported by Chang et al. [1993]) from plants grown under similar 2P or 1P conditions in four independent experiments. Radiolabeled cDNA probes were synthesized from total RNA (30 mg) by RT, as reported (Ramírez et al., 2005). Hybridization and washing conditions of nylon filters were performed as reported (Ramírez et al., 2005). Ten independent nylon filter arrays were hybridized with cDNA from each treatment.
Hybridized filters were exposed to phosphor screens for 5 d, and the fluorescent intensity of each spot was quantified as reported (Ramírez et al., 2005;Tesfaye et al., 2006). To work with highly reproducible experiments, linear regression analysis was performed for each pair of membrane replicas; only those replicas for which the linear model could explain at least 80% of the variation (r 2 $ 0.8) were considered. This process yielded a total of six wellcorrelated replicas for each treatment: 2P roots and 1P roots, respectively. Array data were normalized and quantified using GeneSpring (version 7.2; Silicon Genetics), as provided by the Supercomputing Institute at the University of Minnesota. t tests were performed with a P-value cutoff of #0.05.

RT-PCR Analysis for Verification of Array Analysis
Total RNA for RT-PCR was isolated from 3 g frozen roots using the RNeasy isolation kit (Qiagen). Quantification of transcripts was performed using twostep RT-PCR following the manufacturer's directions (Ambion and Invitrogen) using poly thymine deoxynucleotide primer. The sequences of oligonucleotide primers and conditions used in RT-PCR reactions are shown in Table V. RT-PCR products were resolved in 1% (w/v) agarose gels in Tris-acetate-EDTA buffer, along with a 1-kb DNA-standard ladder (Invitrogen). Amplification of ubiquitin gene was used as control for uniform PCR conditions.

TF Gene Selection and RT-PCR Primer Design
Genes (EST/TC) coding for proteins specifically involved in transcriptional regulation were selected from the TIGR/DFCI Common Bean Gene Index (www.tigr.org). For protein domain prediction, Inter-Pro Release 11 (www.ebi.ac.uk/interpro) was used. The text of all Inter-Pro database entries was searched for the specific strings ''*transcription*'', ''*DNA*binding*'', and ''*zinc*finger*'' using the SRS search tool (www.ebi.ac.uk/interpro/search. html). The identified domains were assembled in a list. The list was supplemented by Inter-Pro domains that are components of the Gene Ontology (GO) branches ''Transcription factor activity'' (GO:0003700), ''Transcriptional activator activity'' (GO:00165643), ''Transcriptional repressor activity'' (GO:0016564), and ''Two-component response regulator activity'' (GO:0000156). The GO-Inter-Pro mappings were found using the QuickGO browser on the Inter-Pro page (www.ebi.ac.uk/ego/). In total, 1,533 domains of proteins potentially involved in transcriptional regulation were selected.
Subsequently, all common bean sequences were analyzed for the occurrence of these domains using Inter-ProScan (www.ebi.ac.uk/Inter-ProScan). In 372 sequences, 41 of the preselected domains were found. The Inter-Pro descriptions of these domains were evaluated to select the domains of proteins that are involved in transcriptional regulation.
RT-PCR primers were generated for the 372 TF genes with TIGR's Primer Design Pipeline, which was designed with the aims of high throughput and specificity. The pipeline iterates through three phases: design, specificity, and selection.
First, the design phase queried every region of the target TF sequences with sliding windows to generate primer set candidates that fit the experimental requirements. Each sliding window was 250 bp across and stepped 50 bp along the target sequence per iteration. The experimental requirements were enforced by the following MIT Primer3 (Rozen and Skaletsky, 2000) parameters: PRIMER_MIN_TM 58, PRIMER_OPT_TM 60, PRIMER_MAX_TM 62, PRIMER_SELF_ANY 6, PRIMER_SELF_END 2, PRIMER_MAX_POLY_X 3, and PRIMER_PRODUCT_SIZE_RANGE#60 to 150#.
Next, the specificity phase aligned primer candidates via WU-Blast (W. Gish, 1996Gish, -2004http://blast.wustl.edu) to the TIGR/DFCI Common Bean Gene Index. The selection phase then discarded primer candidates that registered possible secondary hits, defined as specificity alignments that achieved 80% or greater identity over the length of the primer and included at least one of the terminal ends of the primer in the alignment. The remaining, qualifying primer sets were further prioritized by self-complementarity and poly-X characteristics to achieve selection of the most preferred primers for every target.
The primer design pipeline was implemented in object-oriented Perl modules and supported by a relational MySQL database. Sequences of primer pairs used to amplify each TF gene are shown as supplemental data.

Real-Time RT-PCR Conditions and Analysis of Bean TF
Total RNA for real-time RT-PCR was isolated from 400 mg frozen roots based on the protocol reported by Heim et al. (1993). Three biological replicas were isolated for each treatment (2P and 1P roots), extracting RNA from different sets of plants grown in similar conditions. RNA concentration was measured in NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies), and 10 mg of total RNA was digested with TURBO DNase (catalog no. 1907, Ambion), adding ribonuclease inhibitor (catalog no. N211B, Promega) and following the manufacturer's directions. Absence of genomic DNA contamination was subsequently confirmed by real-time PCR amplification, using primers designed for the bean E2 UBC9 reference gene (TC362; primers: F, 5#-GCTCTCCATTTGCTCCCTGTT-3#; R, 5#-TGAGCAATTTCAGGCAC-CAA-3#). cDNA was synthesized using SuperScriptIII reverse transcriptase (Invitrogen), according to manufacturer's instructions. The efficiency of cDNA synthesis was assessed by real-time PCR amplification of control UBC9 gene. Only cDNA preparations that yielded similar C T values (i.e. 19 6 1) for the reference gene were used for comparing TF transcript levels.
Quantitative determinations of relative transcript levels of TF genes using RT-PCR were carried out at the Max Planck Institute of Molecular Plant Physiology (Golm, Germany) according to Czechowski et al. (2004). Each real-time RT-PCR run for the whole set of TF genes (372) plus reference (housekeeping) and marker genes was performed in a 384-well plate. The bean phosphatase gene (TC201), which is highly induced in P-deficient roots (Ramírez et al., 2005), was included as a P deficiency marker in every real-time PCR run. The primers used for phosphatase gene amplification are: F, 5#-GCCCAAGTTT-GAGGCTGAAAG-3#; R, 5#-TCAAGTCCCACACCGGAAAGT-3#. TF gene expression was normalized to that of UBC9, which was the most constant of the four housekeeping genes included in each PCR run. 2P/1P average expression ratios were obtained from the equation ð1 1 EÞ DDC T , where DDC T represents DC T(2P) 2 DC T(1P) , and E is the PCR efficiency. Student's t test was performed with a P-value cutoff at 0.05.

Plant Metabolite Extraction
Plant metabolite extraction from root samples of 2Pand 1P-treated bean plants and GC-MS metabolite profiling was done as reported previously for Lotus japonicus (Colebatch et al., 2004;Desbrosses et al., 2005). Twelve replicate samples each of roots from plants grown under 1P and 2P conditions were harvested from pods, rinsed with tap water, dried on filter paper, and shock frozen in liquid nitrogen. Frozen samples of 60 to 150 mg fresh weight (FW) were ground 3 min in 2-mL micro vials with a clean stainless steel metal ball (5-mm diameter) using a Retschball mill set at 20 cycles s 21 . Grinding components of the mill were cooled with liquid nitrogen to keep samples deep frozen. Frozen powder was extracted with hot MeOH/CHCl 3 and the fraction of polar metabolites prepared by liquid partitioning into water and further processed as described (Desbrosses et al., 2005).

GC-Time of Flight-MS Metabolite Profiling
GC-time of flight (TOF)-MS profiling was performed using a FactorFour VF-5ms capillary column, 30 m long, 0.25 mm i.d., 0.25 mm film thickness with a 10-m EZ-guard precolumn (Varian BV), and an Agilent 6890N gas chromatograph with splitless injection and electronic pressure control (Agilent) mounted to a Pegasus III TOF mass spectrometer (LECO Instrumente). Details of GC-TOF-MS adaptation of the original profiling method (Desbrosses et al., 2005) are described by Wagner et al. (2003) and Erban et al. (2006). Metabolites were quantified after mass spectral deconvolution (ChromaTOF software version 1.00, Pegasus driver 1.61; LECO) of at least three mass fragments for each analyte. Peak height representing arbitrary mass spectral ion currents of each mass fragment was normalized using the amount of the sample FW and ribitol for internal standardization of volume variations to obtain normalized responses (per gram FW) and response ratios as described (Colebatch et al., 2004;Desbrosses et al., 2005).

Identification of Metabolites within GC-MS Metabolite Profiles
Metabolites were identified using the NIST05 mass spectral search and comparison software (National Institute of Standards and Technology; http:// www.nist.gov/srd/mslist.htm) and the mass spectral and RI collection  of the Golm Metabolome Database . Mass spectral matching was manually supervised, and matches were accepted with thresholds of match .650 (with maximum match equal to 1,000) and RI deviation ,1.0%. Information on the polar metabolites, using the corresponding mass spectral identifiers can be found at http://csbdb.mpimpgolm.mpg.de/csbdb/gmd/msri/gmd_smq.html. Metabolites are characterized by chemical abstracts system identifiers and compound codes issued by the Kyoto Encyclopedia of Genes and Genomes (Kanehisa et al., 2004). Metabolites were identified by standard substances or by MSTs. The term MST is used for repeatedly occurring yet nonidentified compounds, which can be recognized by mass spectrum and RI, as defined earlier (Colebatch et al., 2004;Desbrosses et al., 2005). MSTs are characterized and named by best mass spectral match to compounds identified by NIST05 or Golm Metabolome Database using match value and hit name (Table IV). The response ratio 2P to 1P for each metabolite/MST was calculated dividing the average metabolite concentration from 12 samples from roots of P-deficient plants over the average metabolite concentration from 12 samples from roots of control plants (Table IV).

ICA and Statistical Analysis
ICA (Scholz et al., 2004) was applied to metabolite profiles (as compiled in supplemental data). Data were normalized by calculation of response ratios using the median of each metabolite as denominator and subsequently subjected to logarithmic transformation. Missing value substitution was as described earlier (Scholz et al., 2005). Statistical testing was performed using the Student's t test. Logarithmic transformation of response ratios was applied for approximation of required Gaussian normal distribution of metabolite profiling data.
Sequence data from this article can be found in the GenBank/EMBL data libraries under accession numbers EH791054 to EH791109, EH792671 to EH792678, and EH795233.

Supplemental Data
The following materials are available in the online version of this article.
Supplemental Table S1. Root transcript levels of all the genes in the common bean macroarray.
Supplemental Table S2. Complete list of common bean TF gene and primer sequences.
Supplemental Table S3. Root transcript levels of all common bean TF genes determined by real-time RT-PCR.
Supplemental Table S4. Complete metabolic profile response from common bean roots.