Dissecting nutrient-related co-expression networks in phosphate starved poplars

Phosphorus (P) is an essential plant nutrient, but its availability is often limited in soil. Here, we studied changes in the transcriptome and in nutrient element concentrations in leaves and roots of poplars (Populus × canescens) in response to P deficiency. P starvation resulted in decreased concentrations of S and major cations (K, Mg, Ca), in increased concentrations of N, Zn and Al, while C, Fe and Mn were only little affected. In roots and leaves >4,000 and >9,000 genes were differently expressed upon P starvation. These genes clustered in eleven co-expression modules of which seven were correlated with distinct elements in the plant tissues. One module (4.7% of all differentially expressed genes) was strongly correlated with changes in the P concentration in the plant. In this module the GO term “response to P starvation” was enriched with phosphoenolpyruvate carboxylase kinases, phosphatases and pyrophosphatases as well as regulatory domains such as SPX, but no phosphate transporters. The P-related module was also enriched in genes of the functional category “galactolipid synthesis”. Galactolipids substitute phospholipids in membranes under P limitation. Two modules, one correlated with C and N and the other with biomass, S and Mg, were connected with the P-related module by co-expression. In these modules GO terms indicating “DNA modification” and “cell division” as well as “defense” and “RNA modification” and “signaling” were enriched; they contained phosphate transporters. Bark storage proteins were among the most strongly upregulated genes in the growth-related module suggesting that N, which could not be used for growth, accumulated in typical storage compounds. In conclusion, weighted gene coexpression network analysis revealed a hierarchical structure of gene clusters, which separated phosphate starvation responses correlated with P tissue concentrations from other gene modules, which most likely represented transcriptional adjustments related to down-stream nutritional changes and stress.


Introduction
Phosphorus (P) is an important plant nutrient and essential for plant growth and metabolism. But P is also one of the least available nutrients in soil [1]. Plants take up P only in its inorganic form, phosphate [2]. Free phosphate is present in very low concentrations in the soil solution a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 (< 10 μM; [3]) because it is bound to soil particles, to organic compounds, and forms precipitates with other soil elements like iron (Fe), aluminum (Al) and calcium (Ca) [1].
To acquire and translocate phosphate, plants possess a number of different phosphate transporters (PHTs). In Arabidopsis thaliana, plasma membrane located PHTs of family 1 are responsible for P uptake into the root and show strong up-regulation under P deprivation [4,5]. In poplar, the expression of orthologs of the AtPHT1 genes is also strongly increased upon P starvation [6][7][8][9]. Furthermore, the expression of distinct purple acid phosphatases is enhanced under P deficiency [6]. Phosphatases and organic acids, which are secreted by the roots, increase phosphate availability in soil and together with enhanced P uptake capacities, P acquisition is increased under P starvation [7,10]. Altogether adaptation to limiting P supply results in reduced growth, but higher P use efficiency [6,7]. P deficiency affects not only P uptake and utilization, but also other plant nutrients. For example, in Arabidopsis P starvation caused decreased potassium (K) and enhanced iron (Fe) concentration in the leaves [11][12][13]. Furthermore, reductions in calcium (Ca), magnesium (Mg) and manganese (Mn) were detected; it was speculated that these changes may indicate an adjustment of the ionic charge [11,13] or could be the result of reduced P availability in common P requiring steps of metabolism and transport [12]. There are obviously interrelations among processes that regulate plant nutrient balance because deprivation of tomato in either P or K or Fe resulted in rapid and partly overlapping transcriptional responses [14]. In Arabidopsis, a number of genes, which are highly regulated in response to P starvation [12, [15][16][17], are also required for the utilization of other elements, e.g., SULTR1;3 for sulfate transport and NAS1 for iron chelation by nicotianamine [18,19]. In poplar, reduced P availability leads to higher nitrogen (N) concentrations in roots and lower N concentrations in leaves with higher free amino acid concentrations and reduced activities of N-related enzymes (nitrate reductase (NR), glutamine oxoglutarate aminotransferase (GOGAT), glutamine dehydrogenase (GDH); [6,7]). However, comprehensive analyses how P deficiencies affect the tissue concentrations of essential nutrient elements in poplar and whether these changes are related to alterations in the transcriptome are lacking.
The main goal of this study was to characterize changes in nutrient element concentrations (P, S, C, N, K, Mg, Ca, Fe, Zn, Mn, and the non-essential element Al) in response to P starvation in poplar (Populus × canescens) and to differentiate between transcriptomic changes directly correlated with the plant P concentrations and transcriptomic changes linked to down-stream changes of other nutrient concentrations. For this purpose, we studied poplar transcriptomes under high, intermediate and low P availabilities in leaves and roots, determined weighted gene co-expression networks and identified gene modules that were correlated with nutrient element concentrations. The network showed a hierarchical structure with a top module related to P concentrations that was connected with two main secondary modules, one related to C and N and the other to Mg, S and biomass. The secondary modules were connected by co-expression to tertiary modules. The identified gene modules and their putative functions in the poplar P starvation response are discussed. As the consequence of P starvation the tissue concentrations of other nutrient elements also show massive modification. A significant outcome of our study is that the primary P starvation responses can be disentangled from down-stream gene regulation related to changes in other nutrient elements. three nutrient solutions. Long Ashton nutrient solution [20] with 641 μM P was used for control plants (high phosphate-HP). For intermediate phosphate (MP) availability, the phosphate concentration was reduced to 6.41 μM and potassium added as KCl (675.8 μM), for low phosphate (LP) availability to 0.0641 μM (additional 682.5 μM KCl).
After two months of growth with different P supplies, the poplars were harvested. Aliquots of the plant tissues (leaves, stem, coarse roots, fine roots < 2 mm diameter) were dried at 60˚C for seven days for determination of biomass (n = 10 per treatment): tissue biomass g ½ ¼ dry mass of aliquot ½gÂtotal tissue fresh mass ½g fresh mass of aliquot .
The first three leaves from the top (> 2 cm length) and aliquots of fine roots (< 2 mm diameter) were immediately shock frozen in liquid nitrogen and stored at -80˚C for RNAextraction.

RNA-extraction and microarray
RNA-extraction and microarray procedures at the Microarray Facility (MFT Services, Tübingen, Germany) were described by Kavka and Polle [7]. Three biological replicates (each consisting of two pooled tissue samples) were analyzed per treatment and tissue (fine roots and uppermost leaves).

Statistical analyses of microarrays
Microarray raw data were analyzed using the free statistic software R (version 2.14.2; [22]) after the protocol described by Janz et al. [23]. The R package "affy" [24] was used for normalization of the array probes ("rma" function) using Bioconductor [25]. The log 2 expression value of transcripts that were present ("mas5calls" function) on all replicate chips of at least one condition was used for further analyses. For the annotation of the microarrays, the best gene model for each Affymetrix ID (AffyID) in the annotation file of the Aspen Database was used [26]. Differentially expressed genes (DEGs) were identified by two methods: (a) Significance Analysis of Microarrays (SAM) was performed with R package "siggenes" [27] (FDR 0.05) and (b) the data were analyzed using linear models with R package "limma" [28] (adjusted p-value 0.05). Overlapping DEGs identified by both approaches were used for further analyses. The microarrays have previously been validated by quantitative real time polymerase chain reaction of phosphate transporter expression data [7]. Venn diagrams of DEGs were drawn using the online tool InteractiVenn (http://www.interactivenn.net/; [29]). Log-fold change ratios (logFCs) were calculated with the mean expression values for MP/HP and LP/ HP for roots and leaves.
A weighted correlation network analysis with all DEGs was conducted using R package "WGCNA" (Weighted Gene Co-expression Network Analysis; [30]). Block-wise clustering with a power of 24 was calculated and modules in the resulting dendrogram were merged at a height cut of 0.2. Pearson's correlations of eigengenes, which represent the expression profile of all genes in a given co-expression module, and mean values of element concentrations and biomass in fine roots and leaves of HP, MP and LP plants were calculated. Because the large number of DEGs did not allow a visualization of the co-expression network with the usual tools, gene numbers were reduced by using a p-value (limma) cut-off of 0.00001. Topology was recalculated for these genes by using the same power of 24 (command "TOMsimilarity-FromExpr"). Nodes (DEGs) and edges (co-expression connections between nodes, based on expression pattern) were displayed in Cytoscape (version 3.3.0; [31]) in group attributes layout to show DEGs belonging to one module.

Additional statistical analyses
ANOVA and Tukey's HSD were performed in R (version 2.14.2; [22]) to test for differences between treatments in biomass and element concentrations with p 0.05. Residuals were tested visually for normal distribution and homogeneity of variance and data transformed logarithmically (log 2 ), by square root or by the exponent "2" if needed. Differences between treatments at p 0.05 are shown with different letters in figures. Correlations of mean values were calculated using Pearson's test in R.

Data availability
Raw and normalized data of the microarrays were uploaded into the EMBL-EBI ArrayExpress database (http://www.ebi.ac.uk/arrayexpress/) under E-MTAB-3934. Data for biomass, nutrient elements concentrations and count tables with information on the transcript abundances of all genes are available from the Dryad Digital Repository: http://dx.doi.org/10.5061/dryad. 3t5q5.

P deprivation affects nutrient element concentrations
The concentrations of sulfur, nitrogen, potassium, magnesium, calcium, iron, zinc, manganese and aluminum were influenced by P starvation, but the differences were generally smaller than those for P (Fig 1, Table 1). Sulfur, magnesium and calcium concentrations declined in all tissues except the stem in response to MP or LP treatments (Fig 1, Table 1). Potassium was slightly reduced in all tissues except in fine roots, whereas iron, manganese, zinc and aluminum concentrations increased in most tissues in response to P deficiency (Fig 1, Table 1). Nitrogen concentrations were higher in MP than in HP poplars in all tissues. The carbon concentration did not change markedly under low P availability (Fig 1, Table 1).
Although the changes in element tissue concentrations were caused by differences in P availability in the growth medium, they did not correlate with the P concentrations in the different tissues (Table 2). For other nutrient elements, we detected three groups showing correlations: (i) S, K, Mg, and S with positive relationships, (ii) N, Mg, and Zn with positive relationships and (iii) C, Fe, Mn, and Al, where the latter three were positively related to each other and negative with C (Table 2). Furthermore, some elements were correlated with biomass (P: p = 0.005, r = 0.748, S: p<0.001, r = 0.884; K: p = 0.007, r = 0.731; Mg: p = 0.010, r = 0.706; Ca: p = 0.011, r = 0.700).

P deprivation leads to massive transcriptional reprograming
A total number of 12068 genes were differentially expressed in response to MP or LP treatment compared with HP (Fig 2). Common to all P treatments were 766 DEGs. In general, roots contained less DEGs than leaves although the decline in P was more pronounced in roots than in leaves (Figs 1 and 2). There were 1017 root-specific and 3347 leaf-specific DEGs in response to both MP and LP treatments (Fig 2). The complete list of DEGs, their annotation and foldchanges are shown in S1 Table.

Co-expression modules are correlated with nutrient element concentrations and biomass
To identify genes linked to P nutrition, co-expression modules were calculated and the eigengenes were related to the element concentrations in roots and leaves and the biomass of these tissues. The co-expression analysis resulted in eleven modules (named after colors), which contained 5 to 3632 DEGs (Fig 3, S1 Table). The eigengenes of four modules ("Purple", "Pink", "Red", "Grey") were unrelated to any one of the examined elements or biomass (Fig 3).
Module "Black" showed negative correlations with potassium, manganese and sulfur concentration, and "Magenta" positive correlations with manganese and aluminum and a negative correlation with calcium (Fig 3). Module "Yellow" was negatively correlated with iron and aluminum, and positively correlated with the zinc concentrations. Our analyses also revealed significant regulation of putative transporters with more than 100 annotated DEGs (S1 Table). The transport-related DEGs were spread across different modules without clear patterns related to element abundances.

Network analysis reveals a hierarchy of modules
Network analysis showed that module "Green", which was strongly correlated with P concentrations, was mainly connected with "Brown" and "Blue" by co-expression (Fig 4). Links of "Green" with other modules were rare ("Magenta", "Pink" and "Turquoise") or did not exist ("Red", "Black", "Purple" and "Yellow"). Strong co-expressions were found between genes in module "Blue" with those in modules "Yellow", "Pink", and "Brown" (Fig 4). Module "Brown" was further strongly connected to module "Turquoise", whose genes were co-expressed with genes in modules "Purple", "Magenta" and "Black". The network structure, thus, indicated coexpression of P-related genes with the C-and N-related genes ("Blue") on the one hand and the S-, Mg-and biomass-related genes on the other hand (" Brown", Fig 3, Fig 4). All other modules were apparently not strongly connected with the P-related module, and therefore, the elements correlated with these modules were neither strongly connected with changes in P concentrations.
In addition to "Biological Process" we analyzed the enrichment of GO terms for "Molecular Function" and "Cellular Component" in the co-expression modules. In agreement with "Biological Process" we found GO terms related to the respiratory chain and mitochondrium enriched in module "Green" (S2 Table).
The comparison of the most strongly up-or downregulated DEGs in poplar with a similar study in Arabidopsis [12] revealed a huge overlap of the identified P-responsive orthologs (S3 Table).

Regulation of P responses
Here, we present the transcriptional profiles of poplar roots and leaves in response to phosphate starvation. Leaves of P-starved poplars showed more DEGs than roots. This finding agrees with the P starvation response in Arabidopsis [34] and suggests that the response in leaves is more complex than in roots. Among 12068 DEGs, we identified one gene cluster (module "Green") with 565 genes, which was correlated with the tissue P concentrations. Whether cause-effect-relationships exist between DEGs in module "Green" and the tissue P concentration is unknown, but the module contains orthologs of Arabidopsis transcription factors known to regulate P sensing and signaling [35,36]. It is, therefore, likely that this module is involved in the surveillance of the P status of the cell and in crucial metabolic adjustments because it was also enriched in functions indicating enhanced degradation activities, which may be required to achieve higher P recycling in the plant tissues. Nutrient-related co-expression networks in poplar Among the DEGs in roots and leaves several members of the SPX transcription factor family (highly expressed: Potri.006G069500, Potri.006G253400, moderately enhanced: Potri.014G061200, Potri.018G028200, Potri.006G069500) were detected, which are key regulatory factors of P homeostasis in plants [35,36]. It has already been demonstrated for several plant species that the SPX protein interacts with PHR (PHOSPHATE STARVATION RESPONSE REGULATOR), e.g., AtPHR1 in Arabidopsis, OsPHR2 in rice [37,38]. PHR enables the transcription of P responsive genes by binding to the P1BS element in the upstream region of those genes [37,38]. Arabidopsis mutants with suppressed PHR1 or a homolog transcription factor PHL did not recruit genes enriched in P1BS elements any more while its overexpression resulted in a P starvation response in the presence of sufficient phosphate [39]. SPX modulate the PHR responses: under high P conditions, SPX proteins bind PHR and thereby prevent the activation of phosphate responsive genes in a P dependent manner [37,38]. Upregulation of the expression of SPX genes under low P conditions allows plants to react rapidly to P resupply because the binding of PHR stops the enhanced expression of PHR target genes [37,38]. Because the poplar SPX genes in module "Green" showed P concentration related transcript abundances, it is likely that they fulfill similar functions as their homologs in Arabidopsis and rice [37,38], where they link P perception and signaling.
Further differentially expressed transcription factors known to respond to P starvation were found in other co-expression modules. For example, OsWRKY74 [40] and AtWRKY75 are upregulated during P starvation and regulate P starvation induced genes, whereas AtWRKY6 and AtWRKY42 bind to the PHO1 promoter under P sufficient conditions, inhibiting the transcription of PHO1 (PHOSPHATE1, involved in Pi xylem loading) [35]. Here, the poplar WRKY ortholog Potri.001G097200 (module "Brown") was upregulated upon P starvation in both roots and leaves, whereas others (e.g. Potri.003G138600, Potri.006G263600, Potri.016G128300 from module "Blue", Potri.001G044500 from module "Green", Potri.003G169100, Potri.003G182200, Potri.014G096200, Potri.014G096200 from module "Pink") were downregulated upon P starvation, especially in roots. Different expression patterns for different members of the large WRKY family suggested distinct functions for adaption to low P.
Ethylene-responsive element binding factors (ERFs) were especially abundant in the Nand C-related co-expression module. Upon P deprivation they were downregulated in roots and showed no or slight up-regulation in leaves indicating an influence of ethylene signaling on root responses. This speculation is supported by the observation that low P conditions in Arabidopsis integrate ethylene response with enhanced root hair formation and primary root growth arrest [36]. Furthermore, in Arabidopsis zinc finger proteins (ZATs) such as AtZAT6 have functions in root development and P acquisition [4,36]. Here poplar ZAT orthologs from module "Blue" (Potri.009G027700, Potri.001G235800, Potri.014G017300, Potri.001G295500, Potri.008G051200, Potri.002G119300) were responsive to P deprivation with a strong downregulation in roots upon P starvation suggesting conserved functions of these genes in poplar and Arabidopsis.
An important result of our study was that the identified transcription factors were assigned to different hierarchical levels of the transcriptional network: those being involved in perception and signaling were present in top level (module "Green"), while others, which control long distance transport and morphological adjustments were found in secondary positions (modules "Brown" and "Blue"). Change in transcript abundance of P-responsive genes in roots relative to that in leaves for the modules "Green" (a), "Brown" (b) and "Blue" (c). Log 2 -fold change of MP/HP and LP/HP DEGs in roots were plotted to the corresponding log 2 -fold changes in leaves. Abbreviations of the genes names are listed in S3 Table. doi:10.1371/journal.pone.0171958.g006 Nutrient-related co-expression networks in poplar Transcriptional regulation of P acquisition and P uptake Under P deficiency the expression of genes for phosphate transporters (PHTs) and Purple Acid Phosphatases (PAPs) is upregulated [6][7][8]. Here, upregulated PAPs were present in the P-related module "Green" (Potri.005G233400, Potri.010G158200) but also in the secondary modules "Blue" (Potri.012G042200, Potri.015G031400, Potri.010G158200, Potri.005G233400, Potri.003G030700) and "Brown" (Potri.008G139100). PAPs can be secreted or have intracellular localization. The latter PAPs may contribute to increase the internal P reuse of the plant and are involved in leaf senescence due to their ability to form peroxides [41]. Here, most PAPs were found in module "Blue", in which genes for stress responses and DNA/RNA remodeling clustered as major GO terms suggesting that these invoked enhanced P recycling.
An unexpected result was that none of the PHTs was assigned to module "Green", although various PHTs respond strongly to P deficiency [7][8][9]. Here, PtPHT1;12 and PtPHT1;9 (Potri. 001G318500, Potri.002G005500) were detected in module "Brown", whose eigengenes were strongly correlated with biomass, Mg and S and in which GO terms for cell cycle, chromosome organization and DNA were enriched. The orthologs of PtPHT1;9 in Arabidopsis, AtPHT1;8 and AtPHT1;9, are crucial for P nutrition because they maintain high-affinity P uptake into the roots under P starvation [42]. Also the orthologs in rice, OsPHT1;9 and OsPHT1;10 function in P uptake [43]. Poplar PtPHT1;12 and its Arabidopsis ortholog (AtPHT1;5) show enhanced expression under P deficiency and in senescing leaves [8,44] [44]. Therefore, PtPHT1;12 may function in the redistribution of P inside the plant. Furthermore, in the module "Turquoise", which forms a tertiary level in the co-expression network, further PHTs were upregulated. Module "Turquoise" was strongly negatively correlated with Fe and Mn, which suggests links between the latter elements and P translocation. Earlier studies support this suggestion because enhanced P uptake prevented over-accumulation and toxicity of Mn [45]. P release from plasma membranes replacing phospholipids by galactolipids is a common adaptation to P starvation [46,47]. In our study this response was evident at the transcriptome level. GO term enrichment indicated a correlation between membrane lipid remodeling and P concentrations in module "Green". These results agree with studies in Arabidopsis, where at least parts of membrane lipid remodeling depend on P signaling and auxin/cytokinin crosstalk [48]. For example, the Arabidopsis SQD1 (SULFOQUINOVOSYLDIACYLGLYCEROL1) as well as the rice ortholog are involved in sulfolipid synthesis and regulated by the phosphate content in plant tissues [49,50]. Here, a putative SQD (Potri.016G112600) was one of the most strongly upregulated genes in module "Green" and correlated with the P concentration supporting similar functions in the poplar and Arabidopsis P starvation response.

Metabolism of mineral nutrients, sulfur, nitrogen and carbon upon P starvation
The concentrations of basic cations declined 1.5 to 2-fold (Ca, Mg, K) upon P starvation in poplar leaves, but did not fall below threshold values reported for sufficient nutrition in leaves of young P. × canescens (Ca: 15.9 mg g -1 , Mg: 2 mg g -1 , K: 7.6 mg g -1 [51]). This was also true for sulfur, whose foliar concentrations declined 3-fold under P deficiency (threshold: 3 mg g -1 [51]). Together these elements were positively correlated among each and also with the co-expression module "Turquoise" in which GO terms for sulfur metabolism and organic acid biosynthesis were enriched (S1 Fig). Production and exudation of organic acids are known responses to P deficiency because those metabolites can exchange anions such as Pi from soil particles, thereby increasing Pi availability [52]. The present data imply that the production of organic acids is not a direct response to declining cellular P concentrations but a down-stream consequence that may perhaps also require the perception of general changes in the cellular ion balance.
Cross talk exists between the P sensing and Zn as well as Fe nutrition involving the PHR1 gene regulation [53,54]. In our study increased Fe concentrations were accompanied by enhanced Al and Mn concentrations, whose over-accumulation may damage plant tissues [55] and correlated here with decreased tissue carbon concentrations. Furthermore, in our study the concentration of Zn and Fe were increased or relatively stable and significantly correlated with the eigengenes of module "Turquoise" (negative for Fe) and "Yellow" (positive for Zn, negative for Fe). In contrast to poplar, in Arabidopsis and rice the Zn concentrations are reduced or unaffected during P starvation, whereas the Fe concentrations increased considerably upon P starvation [11,[56][57][58][59][60]. In Arabidopsis, Fe as well as P homeostasis are regulated via PHR1 that induces the iron storage protein FERRITIN1 [53]. Upon P starvation, iron is stored with ferritin in chloroplasts, which prevents precipitation with phosphate [57].
Our data indicate that chloroplastic metabolism must have been massively affected because module "Yellow" was enriched in the GO terms "photosynthesis", "pigment biosynthesis", "starch and carbohydrate metabolism" (S1 Fig). However, no significant up-regulation of ferritin was detected in leaves. In roots, one ferritin gene (Potri.016G124900) was induced, whereas two others (Potri.006G103900; Potri.010G184500) were decreased (S1 Table). In agreement with these observations we also detected up-regulation of an iron transport protein (Potri.015G117900) in roots, whose Arabidopsis orthologs were expressed highly in the outer layers of the roots under Fe-deficiency [61]. In contrast, a putative vacuolar iron transporter (Potri.010G104200) showed decreased expression in poplar roots. These observations suggest that P deficiency modifies iron uptake and distribution in roots, probably similar to that in Arabidopsis [57]. How different species deal with Fe under P starvation may also depend on their ability to produce organic acids which can sequester cations. For example, in strawberry roots the Fe concentration is strongly reduced upon P starvation [62]. Both P and Fe starvation lead to changes in root carboxylate metabolism and an increased exudation of citrate [62]. This suggests similar pathways for P and Fe signaling and independent regulation of Fe metabolism during P starvation in different plants.
In poplar, carbohydrates accumulated in roots upon P starvation [7] in agreement with other studies [63] but the total tissue carbon concentrations were not significantly affected (this study). One reason for this observation could be a shift in metabolic use of carbon-bearing compounds. For example, in P starved poplars, we found a very high up-regulation of PPC (phosphoenolpyruvatcarboxylase) and the PPC activating kinase PPCK. PPC catalyzes the production of oxaloacetate from phosphoenolpyruvate and bicarbonate releasing Pi. Oxaloacetate can then be secreted or further metabolized into other organic acids and used for energy production.
Furthermore, P starvation has consequences for N nutrition. Both carbon and nitrogen clustered with module "Blue", which is characterized by defense responses including respiratory burst (ethylene, phosphorelay signal transduction system) and histone modification suggesting massive transcriptional reprogramming. Unlike most of the mineral nutrients, tissue N concentrations increased, which may happen when growth is inhibited but N uptake continues. Here, we noted a strong induction of BSP (Bark Storage Protein) gene expression. In poplar, BSP expression increases in fall resulting in BSP protein accumulation and N storage during winter [64,65]. Our result suggests that the strong growth reduction imposed by P deficiency led to N accumulation and as a further consequence induced the formation of BSPs to store N. This speculation is supported by earlier results showing that shoot growth influenced BSP promoter activity [66]. Seasonal N cycling in perennial species also involves protein degradation [67]. BSP induction may, thus, be related to protein degradation, which was found here among the directly P correlated transcriptional responses.

Conclusions
P is often the least available nutrient in soil. Here, we demonstrate that P starvation leads not only to P depletion in the tissues but also to massive modifications in the concentrations of other nutrient elements. Under these conditions, in roots as well as in leaves several 1000 genes were differentially expressed. Based on this finding we conclude that a massive acclimation of poplar metabolism is required to cope with this critical situation. A central goal of this study was to find out whether it was possible to separate DEGs related to P starvation and those responding to down-stream events. We successfully identified a hierarchical network structure in which coexpressed genes clustered in modules that correlated with distinct nutrient element concentrations in roots and leaves. The genes in the primary module which comprised only a small fraction (about 5%) of all DEGs highlighted a tight connection between tissue P concentrations and P recycling due to membrane and nucleotide remodeling. The co-expression of genes for metabolic adaptation to P starvation with master switches for the regulation of P uptake such as SPX proposes cross-talk between these processes. Our results further imply similar regulation of these processes in roots and leaves because with few exceptions the DEGs in this module showed the same P starvation response for different tissues. We therefore conclude that the primary module contained mainly those genes, whose transcript abundance is regulated by the tissue P availability.
Surprisingly, PHTs were not present in the P correlated module but occurred in downstream modules. DEGs that represent other typical P starvation responses such as the production of organic acids were neither directly correlated with P, but occurred in a down-stream module that was correlated with sulfur. We, thus, infer that under prolonged conditions of P starvation acclimation of tissue metabolism to P deficiency is the primary requirement and acclimation of P acquisition a down-stream response.
We found correlative links between N enhancement, BSP up-regulation, and growth depression for P starved poplars. Since BSPs are storage proteins with strong seasonal regulation, our findings also highlight that P starvation can interfere with tree-specific physiological processes. Overall, the modular transcriptional responses to P starvation uncovered new links between nutrients and the molecular regulation of biological functions that can be exploited for a better understanding of the plant nutrient balance.   Table. List of most strongly differentially expressed genes in module "Green", "Brown" and "Blue" shown in Fig 6,