Characterization on the P-associated and agronomic traits as well as associated molecular processes in wheat under Pi deprivation condition

Phosphorus (P) acts as one of the essential macronutrients and plays critical roles in regulating growth, development, and yield formation capacity of the crop plants. Elucidating the physiological and molecular processes underlying plant P deprivation responses benefit the crop cultivation with high P use efficiency (PUE). In this study, we investigated the P-associated and agronomic traits as well as the transcriptome profile under contrasting inorganic phosphorus (Pi) levels combined with deficit irrigation, using two wheat cultivars contrasted by PUE (i.e., high PUE Shixin 828 and P deprivation sensitive Jimai 518). The deficient-P (DP) treatment decreased P accumulative amounts, photosynthetic function, and biomass in plants at various growth stages and reduced yields with respect to sufficient-P (SP) condition. Although the two cultivars were comparable on growth and P-associated traits as well as yields under SP, Shixin 828 displayed better growth traits, more P accumulation, and higher yields than Jimai 518 under DP, suggesting that the high PUE cultivars with enhanced P uptake positively affects photosynthetic function, biomass production, and productivity of plants under P deprivation conditions. A large quantity of genes with differential expression patterns, including 2948 upregulated and 1844 downregulated, were identified based on RNA-seq analysis in the Shixin 828 plants treated with P deprivation. These DE genes were shown to be associated with biological process (i.e., metabolic process and cellular process etc.), cellular components (cell body and organelle etc.), and molecular function (binding and catalytic activity), and phytohormone signaling pathways. Transgene analysis on TaZFP1, a gene in ZFP transcription factor family that exhibited upregulated expression in P-deprived plants, validated its role in enhancing P accumulation and P starvation adaptation of plants. Our results suggested that plant P deprivation response is underlying modulation of various physiological processes via modifying gene transcription at global level. RNA-seq sequencing combined with physiological and transgene analyses dissected the molecular processes related to P-associated and agronomic traits in wheat cultivated under Pi deprivation condition.


Introduction
Phosphorus (P) acts as one of the essential inorganic nutrients for plant growth and development, impacting largely on productivities of the crop plants. Given the nature of P nutrition to be prone to immobilization in the growth media, application of the P fertilizers has been used as an effective technique for improving plant biomass production and yield formation capacity for the cereal crops. However, unsuitable management on P fertilizer has resulted in lowered crop P  (PUE) in cultivated fields, such as those of North China as well as other regions with similar ecology, leading to the deteriorated environmental quality aside from increased production cost (Holford et al. 1997;Hinsinger et al. 2001). Therefore, further improvement for P nutrition management benefits the elevation of crop PUE and promotion of the sustainable agriculture worldwide. The P uptake and internal translocation across tissues are accomplished through modified complicate biological processes, such as those associated with P acquisition, cellular P remobilization, and P nutrition recycling during late growth stage. Under Pi deprivation condition, plants are acclimated to the stressor via a subset of evolved mechanisms at molecular level. For example, a quantity of genes functional in diverse physiological processes, such as signaling transduction, transcriptional regulation, primary and secondary metabolisms, protein synthesis, cellular growth processes, and stress defensiveness, modify transcription efficiency under P starvation and are involved in plant low-Pi stress response (Rausch et al. 2002;Raghothama et al. 2005;Jain et al. 2007;Wang et al. 2011;Zhang et al. 2019). Further elucidation of the related molecular processes can help dissect the mechanisms underlying plant Pi starvation acclimation and benefit the breeding of high PUE crop cultivars.
Winter wheat-summer maize cultivation across whole year has been an important cropping system in North China, contributing greatly to the crop production and regional food security. Currently, enhancement of water and nutrient use efficiencies for crops has been an urgent issue to improve the regional sustainable agriculture, given the limitation of water resource and the increase of inorganic fertilizer invest (Xue et al. 2003(Xue et al. , 2006Zhao et al. 2006). Thus far, although the investigations concentrated on elucidation of plant P uptake, P translocation across tissues, agronomic traits under P treatments and the related physiological processes and biochemical pathways have been extensively investigated in wheat (Leikam et al. 1983;Fiedler et al. 1989;Luther et al. 2010), the growth and P-associated traits in wheat cultivars with contrast PUE and the molecular processes contributing to high PUE behavior in wheat under P-deprived combined with water-saving conditions, are needed to be further characterized. In this study, we used two wheat cultivars sharing different P deprivation response, including the high PUE cultivar Shixin 828 and P deprivation sensitive Jimai 518 (about 20 % biomass increase of seedlings in former with respect to that in latter cultivated under hydroponic condition supplemented with low P (0.05 mM Pi), our unpublished data) to address following issues: (i) elucidation of varietal variation on the plant P-associated and growth traits under Pi deprivation combined by deficit irrigation; (ii) characterization of the transcriptome profile underlying high PUE and the associated biochemical pathways in wheat treated with P-and water-saving conditions; (iii) functional characterization of TaZFP1, a gene in the transcription factor (TF) family of zinc finger protein (ZFP), in mediating plant tolerance to low-Pi stress. Our investigation provides insights into plant P deprivation response and PUE-associated molecular processes, and benefits to breeding for high PUE cultivars in wheat cultivated under P-saving and water-limited conditions.

Experimental design
Field experiments were conducted at the Experimental Station of Hebei Agricultural University, Xinji city, China, during the 2017-2018 and 2018-2019 growth seasons. The climate in this region is specified by a typical temperate continental monsoon, with precipitation concentrated at summer season. The meteorological factors during two spring growth seasons are shown in Table S1. The surface soil for experiments was loamy containing following nutrients: organic matter 17.56 g/kg, available N 70.38 mg/kg, available P 20.41 mg/kg, and exchangeable K 121.04 mg/ kg, with a pH 7.66. The treatments were arranged based on a split plot design with triplicates, with plot area of 35 square meters (7 m in length and 5 m in width). Among the treatments, P input level was set up as main plot whereas cultivar Shixin 828, a high PUE cultivar and Jimai 518, a P deprivation sensitive as sub-plot. The P input treatments were set up by two levels, including sufficient-P with P amount 150 kg P 2 O 5 /ha (control, SP) and deficient-P with P amount 60 kg P 2 O 5 /ha (P deprivation, DP). SP and DP were established by using superphosphate (containing 15% P 2 O 5 ) applied prior to seeds-sown, with amount 1000 kg/ha and 400 kg/ha, respectively. Additionally, in total of 225 kg N/ha (urea as N source, with half applied as basal mode and another half topdressed at jointing together with irrigation). In addition, potassium chloride with amount 120 kg/ha K 2 O were applied for all of the plots together with P nutrition application. Deficit irrigation management was set up by two irritations, which were performed prior to seeds-sown (75 mm) and at jointing stage (60 mm, controlled by water amount analyzer). The seeds-sown dates were arranged on October 12 and 14 during the 2017-2018 and 2018-2019 growth seasons, respectively, with seed amounts to generate a 3750-thousand seedling population per hectare. During the growth cycles, the chemical herbicide Tribenuron was used at erection stage (i.e., March 23 and 25 during two seasons, respectively) to control field weeds following the manufacturer's suggestion. cultivation practices conducted were similar to the conventional ones used by the local farmers.

Assay of P-associated traits, photosynthetic parameters, and agronomic traits
At growth stages of jointing, flowering (0 d), mid-filling (20 d after flowering), and maturity, representative plants in each treatment were sampled and subjected to assay of plant biomass, P concentrations, and P accumulative amounts. Among them, plant biomass was obtained based on the oven-dried samples following conventional approach; P concentrations were assessed following the method described by Guo et al. (2011); P accumulative amounts in plants were calculated by multiplying plant biomass and P concentration. Additionally, several photosynthetic parameters, including photosynthetic rate (Pn), stomatal conductance (gs), intercellular CO 2 concentration (Ci), photosystem II photochemical efficiency (ΨPSII), and non photochemical-quenching coefficient (NPQ), were measured as reported by Guo et al. (2013). At maturity, the spikes in two square meters were counted in each plot based on which to calculate population spike numbers. The spike seeds were separately threshed from each plot at maturity (i.e., June 12 and 14 at 2018 and 2019 seasons, respectively) using a mini harvest machine. After air dried, they were weighed and subjected to calculation of grain yields. Grain weights were obtained based on the weight of the air-dried seeds. Spike seed numbers were calculated based on the total seeds counted from thirty representative spikes.

Obtainment of the transcripts datasets upon P deprivation
At flowering stage, representative flag leaves of Shixin 828 plants treated with both SP and DP were collected and subjected to the high-throughput RNA-seq analysis. Total RNA in the samples was extracted using the Plant RNA Extraction Kit (Invitrogen, USA) and the quality of RNA was analyzed based on an agarose gel electrophoresis. In total of 2 µg of high quality RNA derived from the SP-and DP-treated leaves was separately subjected to transcriptome detection using an Illumina sequencing analyzer (Huanuo biotechnology Co. Ltd., China). The transcripts generated that shared high quality were obtained after removal of adapter sequences, low-quality reads, and the reads less than 20 bp using the tool TopHat version 2.1.0, using following default parameters: "-g 1" (uniquely mapped), "-N 1" (allow one mismatch), "-max-intron-length 10", "-max-segment-intron 10", "-fusion-read-mismatches 1", "-fusion-multireads 1" and "-fusion-anchor-length 25" (Kim et al. 2013). Datasets derived from the transcript reads were further assembled using the paired end assembly method implemented in Cufflinks program (Trapnell et al. 2009). Among them, the reads mapped to the genomic regions of reference genome (cv. Chinese spring) were further subjected to expression analysis.
The genes with differentially expressed (DE) pattern under DP were defined based on the transcript counts in DP plants with respect to those in SP ones, using the default parameters of Chi-square test as described by Benjamini and Yekutieli (2001). Of which, the thresholds with P ≤ 0.05 and a |log 2 RPKM ratio| ≥ 1 were used for defining the statistical significance of the gene transcripts.

Gene ontology (GO) annotation and kyoto encyclopedia of genes and genomes (KEGG) analysis for the DE genes
The DE genes identified in wheat plants under DP condition were subjected to GO annotation using the Web Gene Ontology Annotation Plot (WEGO) software. Moreover, the functional categorizations of them were defined as described previously (Ye et al. 2006). KEGG ontology analysis on the DE genes was conducted using the KOBAS 2.0 tool (KEGG Orthology Based Annotation System, v2.0) as reported by Xie et al. (2011).

Transgene analysis on TaZFP1
TaZFP1, a gene of the transcription factor of zinc finger protein (ZFP) family that showed upregulated transcripts under DP, was selected and subjected to functional characterization for its role in mediating plant P deprivation acclimation. h was amplified based on RT-PCR using gene specific primers (Table S2), then integrated into NcoI/BstEII restriction sites in binary vector pCAMBIA3301 at position downstream of the CaMV35S promoter. The target gene was genetically transformed into T. aestivum (cv. Shixin 828) based on an A. tumefaciens-mediated approach as described previously (Guo et al. 2013). Line 2 and Line 3, two T3 lines with more TaZFP1 transcripts together with wild type (WT) were subjected to P input treatments by culturing in standard MS solution (1.2 mM Pi, SP) and modified MS solution containing low P (0.1 mM Pi, DP). Five weeks after the treatments, the P concentrations, biomass, P accumulative amounts, and the photosynthetic parameters in transgenic and WT plants were assessed to be similar to those as mentioned above.

Statistical analysis
Averages of plant biomass, P concentrations, photosynthetic parameters, P accumulative amounts, grain yields, RNAseq datasets, and the qRT-PCR data were all derived from results of three replications. Standard errors on averages and significant differences among the averages were analyzed by using the Statistical Analysis System software (SPSS16.0 statistical software 2, SAS Corporation). All of the data were subjected to statistical calculation using analysis of variance (ANOVA) and to significant test based on the least significant difference (LSD) at P < 0.05 probability level.

Growth and P-associated traits, photosynthetic parameters, and yields
The growth and P-associated traits as well as the photosynthetic parameters at various growth stages during two growth seasons are shown in Fig. 1a and h. Compared with sufficient-P (SP), deficient-P treatment (DP) decreased the biomass, P concentrations, P accumulative amounts of plants (Fig. 1a, c), deteriorated behaviors of Pn, gs, Ci, ΨPSII, and enhanced NPQ (Fig. 1d, h) in tested cultivars managed by deficit irrigation. Although Shixin 828 (high NUE cultivar) was comparable on above growth, P-associated and photosynthetic traits under SP conditions, it was more improved on the traits mentioned than Jimai 518 under DP (Fig. 1a, h). The yields in the cultivars under contrasting P input treatments are consistent with above growth, P-associated, and the photosynthetic traits, with higher shown in Shixin 828 than Jimai 518 under DP (Table 1). These results suggested that the high PUE cultivars possess enhanced yield formation capacity under P deprivation condition, which is largely associated with the improved P accumulation, photosynthetic function, and biomass production of plants.

Transcript datasets generated under P deprivation condition
The transcript datasets in Shixin 828 treated with both SP and DP were established based on RNA-seq analyses. Among the raw base pairs yielded over 58 millions each, more than 51 million were identified to be clean reads. Of which, 88.36-90.22% were mapped to the reference genome (c.v. Chinese spring), with 75.33-76.42% were suggested to be unique and 14.75-15.08% to share a nature of multi mapping characterization ( Table 2). The large transcript datasets obtained in this study suggested the feasibility of RNA-seq approach for detection of transcriptome in wheat under modified P input conditions (Table 3).

The DE genes identified under P deprivation conditions
A large quantity of genes in transcript datasets were differentially expressed (DE) under DP, including 2948 upregulated and 1844 downregulated in the DP-treated plants (Datasets 1 and 2). The scatter plot depicting the expression patterns of the DE genes is shown in Fig. 2. The large set of DE genes detected by the RNA-seq analysis indicates that the plant P deprivation response is comprehensively determined by the modified transcription of genes at global level.

Expression patterns of randomly selected DE genes
Ten of DE genes were randomly selected from the modified transcript datasets under DP and subjected to expression level evaluation. As expected, the five genes with upregulated expression pattern elevated the transcripts under DP with respect to those under SP, with comparable folds of modified expression as shown in the transcript datasets (Fig. 3a, e). Likewise, the five genes with downregulated pattern lowered the expression levels in plants treated with DP compared with those with SP, all of them being similar on transcript changes as shown in the RNA-seq analysis (Fig. 3f, j). These results together validated the reproducibility of the transcriptome results generated from high throughput RNA-seq analysis.

Functional classifications of the DE genes
The DE genes were suggested to be categorized into three functional groups based on GO assignment analysis, including 'biological process', 'cellular components', and 'molecular function' (Fig. 4a). Among them, the genes in the 'biological process' group are enriched by following processes: metabolic process, cellular process, single-organism process, response to stimulus, cellular component organization or biogenesis, biological regulation, regulation of biological process, localization, developmental process, and multicellular organisamal process; the genes in the 'cellular components' group are overrepresented by following constituents: cell, cell part, organelle, organelle part, membrane, macromolecular complex, membrane part, and extracellular region; the genes in 'molecular function' group are closely associated with molecule binding and catalytic activity (Fig. 4a). These results suggested that plant P deprivation acclimation is underlying modulation of numerous modified biological processes regulated by the DE genes, which act in a coordination manner in mediating plant P starvation response.
The Pi deprivation-associated biological processes overrepresented by the DE genes were defined based on KEGG analysis. Results indicated that they are associated with photosynthesis-antenna proteins, carbon fixation in photosynthetic organisms, glyoxylate and dicarboxylate metabolism, phenylalanine metabolism, porphyrin and chlorophyll metabolism, alanine, aspartate and glutamate metabolism, cysteine and methionine metabolism, arginine biosynthetsis, DNA replication, taurine and hypotaurine metabolism, and fructose and mannose metabolism, etc. (Fig. 4b). Therefore, these processes and related biochemical pathways mentioned exert essential roles in modulating plant P starvation adaptation in the high PUE wheat cultivars, through improvement of P uptake, internal P translocation across tissues, photosynthetic function, and biomass production of plants once challenged by Pi depriration combined by deficit irrigation.

Phytohormone signaling-associated genes from DE -ones
The phytohormones, such as auxin, cytokinin, gibberellin, abscisic acid, ethylene, salicylic acid, and jasmonic acid, act as critical regulators and are involved in modulating wide range of physiological processes. Among the DE genes identified, a large set of them were shown to be phytohormone signaling-associated ones, including those involving responses to auxin, cytokinin, gibberellin, abscisic acid (ABA), ethylene, salicylic acid (SA), and jasmonic acid (Table 2). These results suggested the crucial roles of the phytohormone signals in plant P deprivation response through their roles in co-mediation of diverse stress-associated biological processes.

Function of TaZFP1 in regulating P starvation tolerance
Two lines overexpressing TaZFP1 (i.e., Line 2 and Line 3) with more target transcripts (Fig. S1) together with wild type (WT) were cultivated under two P contrasting treatments, to address the target function in regulating plant Pi deprivation response. Under SP condition, Line 2 and Line 3 were comparable on biomass, P concentrations, P accumulative amounts, and photosynthetic traits with the WT plants (Fig. 5a,  h). Under DP, however, the transgenic lines displayed more improvement on biomass, P concentrations, P accumulative amounts, and photosynthetic traits (i.e., Pn, gs, Ci, ΨPSII, and NPQ) of plants than WT plants (Fig. 5a, h). These results suggested that TaZFP1 acts as one of crucial regulators in plant P deprivation adaptation.

Discussion
The cereal cultivars display genetic variation on P uptake, growth, and agronomic traits of plants upon Pi deprivation, which impacts largely on the plant productivity under the stress conditions (Marschner et al. 2005;Seguel et al. 2017). In this study, the biomass, P-associated and photosynthetic traits, and yields in two contrasting PUE cultivars were shown to be in agreement with the previous findings. Although the two cultivars (i.e., high PUE cultivar Shixin 828 and Jimai 518, P deprivation-sensitive) examined were similar on the growth and P-associated traits under sufficient-P condition (SP), Shixin 828 was shown to be much better on above traits than Jimai 518 under the P deprivation combined by deficit irrigation conditions, though the modified precipitation pattern during two growth seasons (136.46 mm in 2017-2018 and 111.5 mm in 2018-2019) could impact Pi uptake of plants challenged with P and water deficit stresses. These results validated the previous reports and indicated that the high PUE cultivars can improve the P uptake and the yield formation of plants under P deprivation treatment . Therefore, adoption of the high PUE cultivars can act as one of effective strategies for improving productivity of winter wheat managed with P-and water-saving cultivations. Plant responses to abiotic stresses are accomplished through the modulation of transcription efficiency of various functional group genes (Sham et al. 2014;Calixto et al. 2018;Zhang et al. 2019). In this study, large quantities of genes were identified to be differentially expressed (i.e., 2848 upregulated and 1844 downregulated) in plants treated with DP based on RNA-seq analyses. These results suggested the complicate nature that plants respond to the P deprivation condition. Previously, it has been confirmed that analyses on DE genes using GO annotation and KEGG can effectively dissect the molecular processes underlying plant stress acclimation (Rehman et al. 2018). In this study, we found that the DE genes identified under DP were categorized into three functional groups based on GO characterization, including 'biological process', 'cellular components', and 'molecular function'. Additionally, the DE genes were suggested to be overrepresented by distinct biochemical pathways based on KEGG analysis, which are associated with photosynthesis-antenna proteins, carbon fixation in photosynthetic organisms, glyoxylate and dicarboxylate metabolism, phenylalanine metabolism, porphyrin and chlorophyll metabolism, alanine, aspartate and glutamate metabolism, cysteine and methionine metabolism, arginine biosynthetsis, DNA replication, taurine and hypotaurine metabolism, and fructose and mannose metabolism, etc. Therefore, the biochemical pathways related to photosynthesis (i.e., antenna proteins, carbon fixation, and chlorophyll metabolism etc.) impact behaviors of the photosynthetic function of wheat challenged by the DP condition. In addition, the pathways related to diverse primary and secondary mechanisms, such as glyoxylate and dicarboxylate metabolism, phenylalanine metabolism, porphyrin and chlorophyll metabolism, alanine, aspartate and glutamate metabolism, and those associated with N metabolism, including alanine, aspartate and glutamate metabolism, cysteine and methionine metabolism, arginine biosynthesis, were shown to be enriched in the high PUE cultivars treated with DP. These results together suggested that plant response to P deprivation is accomplished by integration of various biochemical and molecular pathways.
Phytohormones such as auxin, gibberellin, and cytokinin, are involved in modulation of diverse physiological processes, mediating the plant growth and stress defensiveness upon environmental stressors (Bhargava et al. 2013;Olatunji et al. 2017;Gras et al. 2018). Previously, phytohormones were recorded to act with central regulator in plant stress responses. For example, a set of genes involving auxin signaling (Jain et al. 2009;Wang et al. 2010;Blomster et al. 2011), cytokinin-mediated ROS homeostasis (Liu et al. 2002;Zavaleta-Mancera et al. 2007;Mýtinová et al. 2010), gibberellin metabolism, (Shan et al. 2014), and ABA signaling pathways (Zhang et al. 2014;Ma et al. 2016;Zang et al. 2016), play important roles in mediating plant responses to abiotic stresses. In this study, we found a large set of the DE genes identified in high PUE cultivars treated by DP to be related to phytohormone signaling, such as auxin response, cytokinin metabolism, gibberellin signaling, ABA response, cellular response to ethylene stimulus, SA signaling, and jasmonic acid signaling (Table 2). These results suggested the complicate cross-talk among the phyhormone signaling pathways modulated by DP that further impact on plant P deprivation adaptation. Further characterization of the phyhormone signaling pathways can deepen understanding of the high PUE mechanisms in wheat cultivated under P-and water-saving conditions. ZFP TFs constitute one of the large TF families in plant species (Miller et al. 1985;Klug 2010;Han et al. 2014). ZFP3, ZAT10, and ZAT12, three genes in the ZFP family, play critical roles in regulating plant responses to various abiotic stresses, such as high salinity (Han et al. 2014) and osmotic stress (Mittler et al. 2006) by enhancing osmolyte contents and alleviating electrolyte leakage, via regulation of stress-responsive genes at transcriptional level. Additionally, ZAT6 improves phosphate homeostasis due to regulation of root development and Pi uptake (Devaiah et al. 2007). In this study, TaZFP1, a gene with significantly upregulated in expression in RNA-seq datasets, was selected for functional characterization in mediating P starvation tolerance. Results indicated that overexpression of the target gene conferred plants elevated P concentrations, biomass, and P accumulative amounts under DP. Therefore, TaZFP1 is suggested to be one of the valuable indices for PUE evaluation in wheat cultivars and a useful target for molecular breeding of high PUE cultivars in wheat cultivated under P deprivation and deficit irrigation conditions.

Conclusions
The wheat cultivars examined (i.e., Shixin 828 and Jimai 518) were varied on growth and P-associated traits, photosynthetic function, and yields under P deprivation combined by deficit irrigation conditions. High PUE cultivar Shixin 828 displayed improved P accumulation, photosynthetic parameters, and agronomic traits with respect to Jimai 518 under deficient-P treatment (DP), suggesting its potential for high PUE wheat cultivation managed by P-and water-saving practices. Large quantities of genes were identified to be differentially expressed in high PUE cultivar Shixin 828 challenged by P deprivation, which were associated with diverse biological processes (i.e., metabolic process and cellular process etc.), cellular components (cell body and organelle etc.), and molecular function (binding and catalytic activity). The DE genes are largely enriched in the categories associated with phytohormone signaling pathways, suggesting the critical role of modified hormone signals under P deprivation in plant P deprivation response. TaZFP1, a transcription factor gene in the ZFP family, is functional in positively regulating plant P accumulation, biomass production under DP and conferred plants low-Pi stress tolerance. It is valuable for evaluation of wheat PUE and for molecular breeding of high PUE cultivars in T. aestivum cultivated under low-P and deficit irrigation conditions.

Fig. 5
Plant growth, P-associated traits, and photosynthetic parameters in the transgenic lines overexpressing TaZFP1 under P input treatments. a biomass, b P concentrations. c P accumulative amounts; d Pn; e gs; f Ci; g ΨPSII; h NPQ. WT, wild type. Line 2 and Line 3, two TaZFP1 overexpression lines. Data are shown by averages derived from triplicate results and symbol * indicate to be significant between the transgenic lines and WT under same P treatment (P < 0.05) ◂