Elucidating the unknown transcriptional responses and PHR1 mediated biotic and abiotic stress tolerance during phosphorus-limitation

Phosphorus (P) limitation in the majority of world soils is a major constraint for plant growth and crop productivity. RNA sequencing was used to discover novel P-responsive gene transcripts (PRGT) in leaves and roots of Arabidopsis. Hisat StringTie and Cufflinks TopHat transcript assembler were used to analyze reads and identify 1,074 PRGTs with a >5-fold altered abundance during P-limitation. Interestingly, 60% of these transcripts were not previously reported. Among the novel PRGT, 106 were from unannotated genes, and some were among the most P-responsive, including At2g36727 which encodes a novel microRNA. Annotated novel PRGTs encode for transcription factors, microRNAs, small signaling peptides, long non-coding RNAs, defense-related proteins, and transporters, along with proteins involved in many biological processes. We identified several genes that undergo alternative splicing during P-limitation, including a novel miR399 resistant splice variant of PHOSPHATE2 (PHO2.2). Several novel P-responsive genes were regulated by PHOSPHATE STARVATION RESPONSE1 (PHR1), PHR1-LIKE 1 (PHL1) and PHO2. We discovered that P-limited plants show increased resistance to pathogens and drought stress mediated by PHR1-PHL1. Identification of novel P-responsive transcripts and the discovery of the influence of P-limitation on biotic and abiotic stress adds a significant component to our understanding of plant P-signaling. Highlight Phosphorus limitation elicits the expression of several novel genes including many previously unannotated genes, noncoding RNAs, small peptides and alternatively spliced RNAs, and leads to enhanced disease and drought tolerance.


INTRODUCTION
Phosphorus (P) is an essential element required for plant growth. P is an indispensable component of cellular intermediates in central and energy metabolism, signaling molecules, phospholipids and nucleic acids (Stigter and Plaxton, 2015). P is often the most limiting element for plant growth and development (Holford, 1997;Darcy et al., 2018), and application of P-fertilizers is required to enhance soil P availability and optimize crop yields. However, P-fertilizer use poses a sustainability challenge. Intensive application of P enables high-yield agriculture but also contributes to eutrophication of aquatic systems, while inadequate access to affordable P fertilizer limits crop yields in many developing countries. Moreover, mined rock P, the primary source of Pfertilizer, is a finite resource subject to large price fluctuations and geopolitical risks (Cordell and White, 2014). P sustainability is vital for food security and improvements in plant P efficiency are required through better acquisition or use with the ultimate goal of crop plant species that tolerate lower P availability without reduction of biomass yield or quality. Informed approaches to develop such plants build upon a system-wide understanding of how plants react, respond, and adjust to P limitation. P limitation responses in plants are drastic and occur at different levels (Fredeen et al., 1989;Robinson, 1994;Lynch, 1995;Darcy et al., 2018;Hou et al., 2020). Among these, the transcriptional response may be the most amenable. Using hybridization-based technologies for transcript profiling, such as Affymetrix ATH1 and tiling gene chips encompassing probe sets for over 22,000 genes, steady state levels of many gene transcripts were reported to change significantly during P-limitation in Arabidopsis thaliana (Misson et al., 2005;Morcuende et al., 2007;Bustos et al., 2010;Woo et al., 2012). Similar transcriptome studies were also done in rice (Li et al., 2010;Secco et al., 2013) and soybean (Wang et al., 2016). The three ATH1 gene chip studies by These studies revealed major insights in a plant species' response to P limitation/deprivation, as the functional classification of P-status responsive gene transcripts (PRGTs) indicated adaptations of (i) uptake and transport of Pi and other inorganic ions (sulfate, iron), (ii) Pi salvage systems, (iii) lipid biosynthetic pathways, (iv) primary and secondary metabolism, (v) phytohormone synthesis/response pathways, (vi) cell wall metabolism and (vii) transcriptional and posttranslational signaling.
Despite such major insight, the readout from ATH1 gene chips remained incomplete due to partial and biased gene representation, and limited sensitivity. For comparison, The Arabidopsis Information Resource (TAIR10) genome annotation contains 33,602 genes i.e. roughly 50% more than what is represented in ATH1 Gene chips. Dedicated and highly sensitive approaches such as quantitative reverse-transcription PCR (RT-qPCR) revealed additional information about the transcriptional P-starvation response in Arabidopsis. For example, Morcuende et al. (2007) profiled ~2200 transcription factor genes and other gene families (e.g., phosphate transporters, purple acid phosphatases, PHO1-like and SPX-domain proteins, phospholipases and GDPDs) for their response to P-limitation, revealing P-responsive genes not found or represented on ATH1 arrays. Similarly, Pant et al. (2009) analyzed the expression of ~200 microRNA primary transcripts, none of which were found on the ATH1 array. This study confirmed the strong induction of several miR399 species (Bari et al., 2006;Pant et al., 2008), and also revealed strong induction of miR778, miR827 or miR2111 during P-limitation.
The advent of high-throughput RNA sequencing (RNA-seq) technology provided the opportunity to reanalyze the P-starvation response of Arabidopsis and other plants at truly genome-wide depth, without bias and with sensitivity limited only by the number of RNA-seq reads produced. As yet, no dedicated attempt has been made to comprehensively describe the transcriptional response to P-starvation response in the model plant Arabidopsis with this technology. A major obstacle for doing so may be related to bioinformatics challenges, including the availability of software for read mapping and proper transcript assembly (Pertea et al., 2015) to ensure accurate reconstruction of a transcriptome from the RNA-seq reads. Genome-wide DNA methylation patterns in response to phosphate starvation were determined and correlated with the RNA-seqbased expression levels of P-starvation responsive genes (Yong-Villalobos et al., 2015). However, this study lacked information about the unknown transcriptional P-starvation response. RNA-seq was also performed to compare the response of P-limitation in roots of Arabidopsis wild-type and PHOSPHATE RESPONSE1 (PHR1)-LIKE 2 (phl2) mutant (Sun et al., 2016). This analysis provided a broad comparison of 581 Pstarvation-upregulated (≥2-fold) gene transcripts in the two genotypes at the gene ontology level. The result showed a >30% reduction in expression of a majority (88%) of these phosphate starvation-inducible (PSI) genes in the phl2 mutant. A third study (Yuan et al., 2016) mined RNA-seq data from P-starved Arabidopsis thaliana for long non-coding RNAs with altered abundance but information about other gene classes or novel protein-coding PRGTs was not reported.
To help fill this void, we analyzed RNA-seq data from roots and shoots of P-sufficient and P-limited Arabidopsis plants, grown in sand culture, to identify previously unknown P-starvation responsive gene transcripts in this model plant species. We also compared the obtained expression data with the ATH1 results reported by Bustos et al. (2010), because this study, provided data for both shoots and roots and for all genes represented on the gene chip. We analyzed novel PSI genes for their reversible induction to P-availability by RT-qPCR, their regulation by PHR1/PHL1 and PHOSPHATE2 (PHO2), and screened for P-responsive alternative splicing forms. We analyzed the novel PRGTs for their involvement in different biological processes and cross-talk with other biotic and abiotic stresses and we analyzed P-limited Arabidopsis plants for resistance to a bacterial pathogen and drought stress.

Plant materials and nutrient media
Arabidopsis thaliana Columbia-0 wild-type (WT) seeds were surface sterilized and grown in half strength Murashige and Skoog (MS) nutrient agar plates for 10 days under short-day conditions with light/dark regime of 8/16 hours in a Percival growth chamber (www.percival-scientific.com). Seedlings were transferred to acidwashed sand and perlite mixture and grown in a growth chamber with a light/dark regime of 12/12 hours at 23/22 ºC, respectively, with a light intensity of ~125 µmol m -2 s -1 and a relative humidity of 50%. The plants were divided into two groups and were watered with the nutrient solution with full nutrition (P-sufficient, +P) and without phosphate (-P) for the first three days and then the nutrient solution was diluted 5X and supplied to the plants every other day. Plants were grown under -P condition for nine days until the typical P-limitation phenotype (dark-green leaves and smaller plants) appeared. The location of the plants within a growth chamber were randomized every two days to minimize position effects. Pi-limitation in plants was confirmed before harvesting by measuring inorganic phosphate in the leaves. The roots and shoots were washed with deionized water, harvested, snap frozen in liquid nitrogen and stored at -80 ºC until further analysis. The nutrient media contained 4 mM KNO 3 , 2 mM NH 4 NO 3 , 3 mM KH 2 PO 4 /K 2 HPO 4 , pH 5.7, 4 mM CaCl 2 , 1 mM were grown in axenic liquid culture as described previously (Pant et al., 2015).

Phosphate (Pi) measurement
Measurement of the inorganic phosphate was done using the colorimetric micro method (Itaya and Ui, 1966) as described in Pant et al. (2015). In summary, frozen plant material was ground in liquid nitrogen, weighed and soluble inorganic phosphate was extracted in Milli-Q water (Millipore) by repeatedly (3X) freezing and heating the samples at 95 °C for 5 minutes. 10μL of the sample was mixed with 100μL HCl and 100μL malachite green color reagent in a 96 well plate. The plate was incubated at room temperature for 15 min. and the absorbance was measured at 630 nm. Phosphate concentration in the samples was determined based on a calibration curve.

RNA isolation, cDNA synthesis and quantitative reverse transcription PCR
RNA was isolated using TRIzol reagent (Invitrogen) according to the manufacturer's standard protocol. For RT-qPCR and RNA-seq experiments, residual DNA was removed from RNA samples by treating with TURBO DNA-free™ DNase (Ambion) according to manufacturer's instructions. The complete removal of genomic DNA (gDNA) from RNA was confirmed by RT-qPCR using the intron or exon primers, where gDNA was taken as a positive control and water sample as a negative control. NanoDrop spectrophotometer (ND-8000) was used to check the concentration and purity of RNA, and Agilent Bioanalyzer 2100 was used to check the RNA integrity.
RNA samples with RNA Integrity Number (RIN) values >8 were taken for RNA-seq. Selection of primers, RT-qPCR conditions and procedures for cDNA synthesis were as described previously (Pant et al., 2008;Pant et al., 2015). DNA oligos used in the study are given in supplemental table S8.

RNA-seq library preparation and deep sequencing
Three independently prepared batches of plants grown under P-sufficient and P-limited condition were used for

RNA-seq read mapping and data processing.
Paired-end reads generated from the RNA-seq were used for each sample. All RNA-seq reads were demultiplexed using a known list of barcodes (Illumina), allowing zero mismatches. After demultiplexing, adapter sequences were removed from the sequenced libraries using the FastQ/A clipper found in the FastX Toolkit (http://hannonlab.cshl.edu/fastx_toolkit/). Reads were quality trimmed by removing low-quality bases until two consecutive bases with quality scores of 30 or above were found. Only the quality trimmed reads longer than 50 nucleotides were used for alignment with the Arabidopsis reference genome version 10 performing RT-qPCR expression analysis for some selected genes/loci.
To obtain more inclusive list of differentially regulated genes in response to P-limitation, the RNA-seq reads were also analyzed using an independent method. Illumina adapters and low-quality bases were trimmed using  , 2015). SAM files were converted into compressed and sorted BAM files using Samtoolsview, -bSh, and -sort-sort commands, respectively (Li et al., 2009). Reads mapping to annotated and unannotated regions of the Arabidopsis genome were extracted using StringTie (Pertea et al., 2015). Genes with less than one read count per million across at least three samples were discarded from the analysis. For the differential expression analyses, significant transcripts were converted to the corresponding gene for performance evaluation, such that if a single transcript was called as differentially expressed, the corresponding gene was also called differentially expressed. We note that because of this unavoidable difference between gene-level and transcript-level comparisons, quantitative comparisons of recall and/or precision between a gene-level and a transcript-level workflow should be avoided. Rather, we recommend evaluating the relative performance of a given workflow as compared with other workflows with matched genelevel or transcript-level estimation.

Water Loss Assay
For drought tolerance analysis, water loss assay was performed as described previously (Pant et al., 2022). pv. tomato DC3000 at concentration of 1 X 10 8 CFU/ml in water containing 0.01% Silwet L-77 and the plant phenotype was recorded after 5 days. Pathogen assay for Erwinia carotovora ssp. carotovora and P. syringae pv. maculicola was performed as described previously (Ishiga et al., 2011;Pant et al., 2020). Briefly, 3-weekold Arabidopsis plants grown on ½ MS media plates (+P and -P condition) were flood inoculated at concentration of 1.6 X 10 5 CFU/ml in water containing 0.01% Silwet L-77. Quantification of bacterial growth was done at 3 dpi for E. carotovora and 5 dpi for P. syringae pv. maculicola by counting colony forming units (CFUs).

RNA-seq study identified several known and novel P-limitation responsive gene transcripts (PRGTs)
Plants grown in P-limitation (-P) were distinctly smaller than plants grown in P-sufficient (+P) conditions and showed anthocyanin accumulation (supplemental Figure S1 A, B). Inorganic phosphate content in the shoots of -P grown seedlings was about 10-fold lower than in +P grown seedlings (supplemental Figure S1 A Figure S1C). RNA-Seq library production and analysis resulted in 10 RNA-seq files (supplemental data files 1-10). Analysis of the RNA-seq data was performed using HiSat StringTie and the TopHat Cufflinks pipelines using the stringent approach.
The output file from the HiSat StringTie pipeline (supplemental Table S1 . In the shoot, 50 PRGTs were determined to be induced by >100-fold ( Figure 1B). Among those 50 PRGTs, 28 are represented on ATH1 gene chips and are known PRGTs. The remaining 22 include some additional known PRGTs identified by other approaches (e.g. IPS1, miRNA399d), along with several annotated as "unknown protein" or "other RNA" and six that are not annotated in TAIR10 ( Figure 1B).

The filtered
≥ 5-fold changed PRGTs were manually inspected for their response in the Integrative Genomics Viewer resulting in 909 and 861 PRGTs in the HISAT-StringTie and TopHat Cufflinks output files, respectively, with an overlap of 696 and a non-redundant total of 1,074 PRGTs (supplemental Table S1, supplemental Figure S2). Reliance on only one of the two RNA-seq analysis pipelines would have precluded the identification of some known PRGTs already identified in ATH1 studies. For example, SPX DOMAIN GENE were identified with HiSat StringTie, whereas AT2G17280, GALACTINOL SYNTHASE 1 were extracted with the TopHat Cufflinks pipeline only (supplemental Figure S2A). Both analysis pipelines detected 405, 116 and 138 GTs with strong P-response in shoots, in roots or in both shoots and roots, respectively.
Of the 1,074 non-redundant PRGTs, 802 (74.7%) were induced and 272 (25.3%) were repressed during P-limitation ( Figure 1C and S1B, supplemental  (Figure 2A, green sector). In addition, 259 of the 1074 >5-fold changed PRGTs (Figure 2A, pink sector) detected by RNA-seq correspond to TAIR10 annotated genes for which no probe sets are present on ATH1 arrays, and yet another 106 PRGTs (yellow sector) were found for which no gene models exist in TAIR10. Altogether, this RNA-Seq study reports 648 (= 283+259+106; red, pink and yellow sectors respectively) ( Figure 2A) strongly differentially expressed PRGTs for which no evidence was previously obtained in Bustos et al. (2010) and the P-responsive status of many was previously unknown. Figure 2B shows Sashimi plots for three PRGTs currently not annotated in TAIR10 (yellow sector), tentatively named AT5G27993, AT1G53605 and AT1G53615.
RT-qPCR was used to verify the P-limitation response of thirty randomly selected gene transcripts from the green (4), pink (19), yellow (5) blue (1) and red (1) sectors ( Figure 2A, supplemental Figure S5). Shoot and root samples from an independent set of sand culture-grown plants were investigated for this purpose. The plants were either supplied with phosphate (+P) or were P-limited. P-limited plants were also re-supplied with phosphate for 30 min or 3 hours before harvesting. RT-qPCR demonstrated that virtually all thirty tested gene transcripts again displayed a strong P-response (supplemental Figure S5A, B). Re-addition of phosphate to Plimited plants only slightly affected the abundance of these transcripts within 30 min but led to a marked reversion for most of the gene transcripts within three hours, indicating that these PRGTs probably directly respond to a change in P status rather than an indirect and delayed change in metabolism or development.
More detailed comparison of the RNA-Seq and RT-qPCR data of these 30 gene transcripts also revealed good quantitative correlation (supplemental Figure

Sand-grown, P-limited Arabidopsis plants display many of the familiar transcriptional P-limitation responses known from plate-or liquid culture-grown younger seedlings
In the present study, we investigated 19-day old Arabidopsis thaliana Col-0 plants grown in a sand/perlite mix in phytotrons (Figure 1 cf. materials and methods). In terms of growth conditions and age, there is a considerable difference between the plant materials used in this study and the plate-or liquid culture-grown younger seedling materials used in previous ATH1 P-limitation studies (Misson et al., 2005;Morcuende et al., 2007;Bustos et al., 2010). Nonetheless, many previously identified PRGTs were also identified here. For example, of the 100 most strongly P-limitation induced gene transcripts detected in shoots of plate-grown seedlings by Bustos et al. (2010), 62, 24 and 12 were present in the green, orange and brown sectors, respectively ( Figure 2A). Consequently, the major (hallmark) transcriptional responses to P-limitation were also identified in our study (supplemental Table S1). These hallmark transcriptional responses include induction of 1 gene transcripts for (i) uptake and transport of Pi and other inorganic ions (Pht1.3/1.4/1.  Table S1).

Arabidopsis plants
To obtain a broad functional overview of the novel PRGTs in the red and pink sectors (Figure 2A), their Arabidopsis Gene Identifiers (AGIs) were subjected to a Gene Ontology (GO) analysis in TAIR ( Figure 2C).
Compared to the 229 GTs that were found to be strongly P-limitation induced in the present RNA-seq study and in the Bustos et al., 2010 ATH1 study (green sector), the 283 GTs in the red sector (identified from this study, Figure 2A) displayed a notably increased frequency of GO terms "transcription factor" (7.0%), "DNA/RNA binding" (10.2%) and "unknown molecular function" (12.8%), while the terms "hydrolase activity", "transferase activity" and "transporter activity" were reduced. The abundance of GO term ""unknown molecular function" further increased (to 33.2%) in the set of 259 PRGTs in the pink sector ( Figure 2A, novel transcript identified from this study without ATH1 probe set on ATH1 array), while the terms "hydrolase activity" (8.2%) and "transferase activity" (6.8%), further declined.
Next, individual inspection of the PRGTs in the pink and red sectors (

Novel PRGTs encoding transporters
Based on ATH1 Affymetrix gene chip studies, P-limitation is known to affect expression of ion transporters, most notably to induce PHT1 phosphate transporters, but also potassium and sulfur transporters, indicating the profound effect of P-limitation on ion homeostasis. PRGTs in this study with GO annotation "transporter activity" in the blue, red and pink sectors ( (AT5G50200) were strongly repressed in roots, suggesting that P-limitation leads to a decrease of nitrate uptake. NRT1.1 and NRT2.1 also showed repression in shoots, whereas NRT1.7 (AT1G69870) and NRT1.8 (AT4G21680) were induced. Repression in P-limited roots was also found for AT4G08620, AT1G77990 and AT1G78000 encoding sulfate transporters 1.1, 2.2 and 1.2, respectively, whereas induction of sulfate transporters 1.3 (AT1G22150) and 3.4 (AT3G15990) is known from ATH1 Affymetrix gene chip studies and was confirmed here (see above). In addition to these changes, the RNA-seq analyses also revealed novel PRGTs encoding other cation/anion, peptide, sugar, ABC, MtN21/UMAMIT and MATE transporters (Table 2).

Novel PRGTs that encode transcription factors
PRGTs annotated as "encoding transcription factors" were analyzed in more detail, due to their importance and because of their abundance (33, i.e. 11.7%) in the red sector ( Figure 2C and supplemental Table S2), suggesting that the sequencing depth and sensitivity of RNA-seq was high enough to detect gene transcripts with low average expression (supplemental Figure S3), such as those encoding transcription factors previously missed using ATH1 arrays. In the green, blue and pink sectors (figure 2A) 11 (4.8%), 19 (9.6%) and 11 (4.2%), respectively, PRGTs have the GO annotation "transcription factor" yielding a total of 74 P-responsive transcription factor gene transcripts (TFGTs) identified by RNA-Seq. Sashimi plots of eight novel P-limitation induced and two repressed TFGTs, four (AT1G18860, AT1G49900, AT1G60280, AT2G46790) from the red and four (AT1G27045, AT3G18010, AT3G28857, AT5G55690) from the pink sector ( Figure 2A) are shown in Figure 3. While the RNA-seq results confirmed the existing annotations of genes like WUSCHEL-related homeobox AT3G18010 and AT2G46790/PRR9 ( Figure 3B, G), the Sashimi plots for AT5G55690 ( Figure 3A), the highly induced MADS box TF gene transcript AGL47, previously identified as P-responsive by RT-qPCR (Morcuende et al., 2007), and AT1G18860 ( Figure 3E), a WRKY TF, both revealed an additional unannotated exon. Sashimi plot for AT1G49900 (C2H2 zinc finger TF; Figure 3D) suggests that the current annotation encompasses two distinct PRGTs. AT1G27045/HB54 ( Figure 3F) and AT1G60280/NAC23 ( Figure 3H) are two examples of TFGTs with very low expression, yet their P-response was clearly identified using RNA-seq (supplemental Table S1).
Besides two TFGTs encoding nuclear factor Y subunits (AT1G30500, AT3G05690) from the green sector, four additional ones (AT5G43250, AT1G54160, AT1G72830 and AT2G34720) were present in the blue and pink sectors. Furthermore, four TFGTs (AT2G46790/PRR9, involved in clock and circadian function were considerably induced during Plimitation ( Figure 3G). In addition, AT1G09530/PIF3 which interacts with promoter elements of LHY and CCA1 and with photoreceptors PhyA and PhyB, and AT3G62090/PIF3-like2 which physically associates with TOC1/PRR1, a pseudo response regulator involved in the generation of circadian rhythms, were induced. It thus appears that P-limitation can affect clock and circadian functions through several of the known, central regulators.

Novel PRGTs involved in signaling and regulatory processes
In addition to the above mentioned PRGTs, we sought for other novel potential regulatory gene transcripts that respond to P availability. These included gene transcripts encoding (i) regulatory RNAs (microRNAs, longnoncoding RNAs), (ii) proteins with functions in signaling, including receptor kinases, small secreted peptides, and MAP kinases, and (iii) proteins with predicted functions in post-translational protein modification (kinases and phosphatases) and protein degradation for example F-box proteins and E3 ligases (Table 1).
Primary transcripts of sixteen microRNAs were detected among the 5-fold changed gene transcripts (supplemental Table S1, Table 1). While induction of miR399b-f, miR827, miR2111a transcripts was very strong and expected (Pant et al., 2009), the induction in roots of four primary transcripts encoding miR156a, c, e and miR157a (all targeting SPL transcription factors) and repression of two transcripts encoding miR169a/b (targeting NF-YA HAP-type transcription factors) in shoots and roots, respectively, was less pronounced and unreported. Repression was also observed for the primary transcripts encoding miR161, miR823a and miR775a, while the one for miR472a was induced. Besides microRNAs, 11 long non-coding RNAs (lncRNAs) were identified among the genes in the pink sector (supplemental Table S1, Table 1), including the known IPS1, AT4/IPS2 and TAS4, four genes annotated as "potential natural antisense RNA" (AT3G12502, AT5G01732, AT5G44562, AT5G53048) and another four annotated as "other RNA" (AT1G21529,

AT3G57157, AT2G14878, AT1G43765).
Gene transcripts encoding small secreted peptides (SSPs), known to play a role in many signaling pathways, were (and continue to be) poorly annotated and thus hardly represented on ATH1 gene chips.
During recent years the importance of SSPs emerged especially with regard to their crucial developmental functions, and more recently also with regard to their role in plant nutrition (Tabata et al., 2014;de Bang et al., 2017;Ohkubo et al., 2017). The roles of SSPs during P-limitation are yet uncharacterized and P-status responsive SSP gene transcripts are thus far almost undescribed. We detected fifteen SSP transcripts among the 1,074 PRGTs (supplemental Table S1; Table 1), including those encoding root meristem growth factor RGF7/GLV4, rapid alkalinization factor RALFL19, EPF2, PIP1 and PIP3, CLE19, DVL1 and three plant defensins. Moreover at least two PRGTs (AT1G08165, AT3G13433) encoding potential novel peptides were identified among the 1,074 PRGTs. Inspection of PRGTs with responses <5-fold ( Figure 1) yielded a total of 38 annotated SSP transcripts with at least 2-fold changed abundance in shoots and roots during P-limitation. It will be important to investigate how P-status responsive SSPs, either produced in planta or synthetically, affect the plant P-limitation response including root development, and what their cognate peptide receptors (usually LRR-RLKs) are.
Over 20 additional PGRTs with potential regulatory functions were identified in this study as shown in the blue, red and pink sectors (Figure 2, supplemental Table S1; Table 1). For example, AT4G26890, encoding MAP3K16, was ~9-fold induced in shoots, and, in addition to MAP3K18 and 19 (AT1G05100, AT5G67080) which are present in the green sector, represents the third MAP3K responsive to P-status. The transcript encoding calcineurin B-like-interacting serine/threonine protein kinase 14 (CIPK14, AT5G01820) was known from ATH1 studies to be P-limitation induced in shoots, and this induction was confirmed here. More interestingly, five additional CIPK transcripts (4, 12, 18, 19 and 20; Table 1) were identified in this study to be induced by P-limitation in shoot (CIPK12,18,19,20) or root (CIPK4), thus suggesting that several signaling pathways involving CIPKs calcium sensors are modulated by P-status. Furthermore, at least six novel genes encoding P-responsive F-box proteins, a Skp1-like ubiquitin protein ligase, a RING-finger E3 ubiquitin ligase, two protein phosphatases 2C, as well as PYR1-like 6, a regulatory component of the abscisic acid receptor, were found to respond to P-status at the transcript level (supplemental Table S1; Table 1).

Novel PRGTs encoding cell wall, peptidases and domain of unknown function proteins
Over thirteen novel PRGTs related to cell wall, peptidases, domain of unknown function (DUF506) and unknown proteins were identified in this study (supplemental Table S1; Table 1). This includes four gene transcripts in red sector such as AT2G20750/Beta expansin and AT2G20870/cell wall protein that were repressed in shoots and AT4G30280/xyloglucan endotransglucosylase and AT2G26620/pectate lyase that were induced in shoots during P-limitation. Differential regulation of these cell wall related genes fits well with fine-tuning of cell wall formation during P-limitation. Five additional PRGTs encoding peptidases (AT5G36180/SCPL1, AT2G24010/SCPL23, AT4G15100/SCPL30, AT2G27920/SCPL51, AT5G22860/Serine carboxypeptidase S28 family) were identified in this study and included in blue, red and pink sectors (supplemental Table S1; Table 1). P-limitation induction of previously known AT2G24000/SCPL22 was confirmed here. Other novel PRGTs identified by this study include DUF506s or unknown protein encoding gene transcripts such as AT3G25240, AT4G32480, AT2G20670, and AT1G03106 as shown in blue and pink sectors (supplemental Table S1; Table 1, Figure 2A). AT3G07350/DUF506 was known from ATH1 studies to be P-limitation induced in shoots and roots, and this induction was confirmed here.

Unannotated transcriptional units, long non-coding RNAs, and a novel P-responsive microRNA
Unannotated transcripts (106 in total) with >5-fold change in abundance during P-limitation were identified ( Table S1). The unannotated transcripts have low or no coding potential, i.e. only small or no open reading frames are discernible, thus classifying them, a priori, as long non-coding RNAs (lncRNAs). A comparison with the set of ~1200 lncRNAs suggested by Yuan et al.
LncRNAs have various biochemical and molecular functions as target mimics (e.g. IPS1) (Franco-Zorrilla et al., 2007), small RNA precursors, antisense RNAs (blocking recognition of exons by the spliceosome; generating endo siRNAs) or binding partners of proteins (thus affecting protein activity, alter protein localization or providing a structural/organizational role) (Wilusz et al., 2009). Although in-depth structural analyses of the lncRNAs or their functional assignments are not within the scope of this study, the highly P-limitation induced and expressed MSTRG.19189 (AT5G27993) transcript was found to be produced from a genome region that can form an extended stem-loop structure ( Figure  genes were up-regulated in pho2 mutants in response to P. Since PHR1 is known to regulate phosphate starvation-inducible (PSI) genes by binding to their PHR1 binding site (P1BS) element that is present in the promoters of most PSI genes (Rubio et al., 2001), these 30 genes were also inspected for the presence of P1BS element in their 3000 bp upstream region. Nearly 25 out of 30 genes (83%) had PHR1 binding site in their upstream sequence and 87% of these genes were found to be down-regulated in phr1phl1 mutants ( Figure 6, supplemental table S5). These results indicate that the majority of these PRGTs act downstream of PHR1/PHL1 and PHO2 signaling pathway.

Alternative splicing of PRGTs and identification of miR399 resistant PHO2.2
Inspection of the genome using splice junction mapper TopHat2 v2.0.9 (Kim et al., 2013) in 100 nt windows and looking for significant changes in coverage within genes or hotspot list identified almost 100 additional PRGTs, many significantly improved annotations and some alternative splicing events (supplemental Table S6 and S7). GO analysis of the 44 alternately spliced genes revealed chromatin-binding or -regulatory protein as the most enriched GO category followed by gene-specific transcriptional regulator, membrane traffic protein, metabolite interconversion enzyme, nucleic acid metabolism protein, protein modifying enzyme, transmembrane signal receptor, and transporter (supplemental Figure S9). Among these alternately spliced PRGTs was PHO2 (PHOSPHATE 2), which has a splice form that encodes only about half the transcript by splicing out of the miR399 binding sites, 1 st , 2 nd and 3 rd introns as well as 2 nd exon of 5' UTR and 1 st exon of the coding sequence ( Figure 7A). As PHO2 is known to play a critical role in Pi-signaling pathway and homeostasis (Bari et al., 2006;Pant et al., 2008), its alternative splicing was further investigated in detail.
Interesting splice variants of PHO2 (PHO2.1 and PHO2.2) were detected using IGV, therefore six pairs of splice form specific primers were designed to verify its existence by PCR and RT-qPCR. Data from the PCR based assay confirmed the existence of PHO2.1 and PHO2.2, a larger and smaller splice form of PHO2, respectively ( Figure 7B). Furthermore, transcript abundance of PHO2 was monitored with splice form specific RT-qPCR primers ( Figure 7C, D). PHO2.2 (shorter splice form of PHO2) has relatively higher abundance in Plimited shoot, whereas it is almost undetectable in shoots of FN grown plants ( Figure 7C) but in roots that form is equally abundant in both conditions. This shows that PHO2.2 level is not affected by miR399 abundance, which is very high in P-limited roots. Splice form PHO2.1 is more abundant in P-sufficient root compared to Plimited root, this fits with the miR399 mediated negative regulation of PHO2 during P-limitation. But in shoot, miR399 mediated negative regulation of PHO2 does not take place (Bari et al., 2006;Pant et al., 2008). The splice form PHO2.2 is equally abundant in P-sufficient and P-limited root, fitting with the model that PHO2.2 is resistant to miR399-mediated transcript cleavage, as it does not have a miR399 binding site. Intriguingly, PHO2.2 is undetectable in P-sufficient shoot but more abundant in P-limited shoot ( Figure 7A1, D). This suggests a novel mechanism involving negative regulation of PHO2.2 in P-sufficient shoot.

Phosphate limitation modulates plant defense via PHR1-PHL1
Phosphate deficiency was previously found to change in the expression of jasmonic acid (JA) biosynthetic (e.g. AT3G25760) or signaling genes (e.g. AT3G12500) (Morcuende et al., 2007;Bustos et al., 2010) and to induce the jasmonate pathway enhancing resistance to insect herbivory (Khan et al., 2016). JA induction and P starvation also share some common phenotypes, such as growth reduction and anthocyanin accumulation.
Therefore, the response of gene transcripts involved in defense responses was analyzed in more detail. carotovora. The results showed that P-limited plants grown either in soil or on agar were resistant to all the pathogens tested as evidenced by reduced necrotic/chlorotic spots when compared to P-sufficient plants ( Figure 8B-G). In addition to the disease symptoms, the bacterial growth curve assay also showed significantly less bacterial accumulation in P-limited Arabidopsis plants when compared to P-sufficient plants. Pathogenicity assay also showed that phr1 -P and phr1-phl1 -P were more susceptible to these pathogens when compared to wild-type -P plants ( Figure 8B-G). These results indicate that P-limitation enhances disease resistance in Arabidopsis and is regulated via PHR1-PHL1.

P-limited plants transpire less water compared to P-sufficient plants and is mediated by PHR1-PHL1
We found that AT4G27920, which encodes a member of the PYRABACTIN RESISTANCE (PYR)/PYR1-like (PYL)/REGULATORY COMPONENTS OF ABA RECEPTOR (RCAR) family proteins was induced in both shoots and roots during P-limitation. Plant induction of ABA during drought stress is sensed by PYR/PYL/RCAR receptors and these in turn regulate the activity of protein phosphatases 2Cs ABI1/ABI2 and downstream signaling pathway during drought stress (Santiago et al., 2009;Cutler et al., 2010;Yu et al., 2016). AT1G69260, which encodes ABI FIVE BINDING PROTEIN (AFP1) was also induced in shoot of Plimited plants.
To further determine the cross talk between P-limitation and drought, we identified the genes confirmed to play role in drought based on literature search in TAIR10 database (www.arabidopsis.org) and found that 11 of these genes were induced by P-limitation by more than 5-fold ( Figure 9A). GO term enrichment for these genes revealed that majority of these genes are transcriptional regulators and transporters. Among these, AT1G21529, which encodes a DROUGHT INDUCED lncRNA (DRIR), was ~300 fold induced by P-limitation.
Overexpressing DRIR in Arabidopsis was shown to increase drought tolerance and ABA sensitivity (Qin et al., 2017). MiR399f, a highly P-limitation induced gene, also modulates plant responses to ABA and drought (Baek et al., 2016), in addition to its role in phosphate signaling. AT5G09570 (AT12CYS-2), which is highly induced by P-limitation ( Figure 9B), encodes a twin CX9C domain protein and loss of both AT12CYS-1 and AT12CYS-2 was shown to enhance drought tolerance (Wang et al., 2016). We performed RT-qPCR expression analysis of eight genes common between drought and P-limitation and found that seven out of eight genes were regulated via PHR1-PHL1 transcription factors ( Figure 9B). However, AT12CYS-2, which is a negative regulator of drought tolerance, was also induced in phr1 or phr1-phl1 mutants during P-limitation.
Based on all these results, we hypothesized that P-limitation can enhance drought tolerance in plants. Water loss assay in the detached leaves of Arabidopsis WT, phr1 and phr1-phl1 mutants was performed and it confirmed that P-limited plants lost the least amount of water compared to P-sufficient plants, while P-limited phr1 and phr1-phl1 mutants lost more than WT-P ( Figure 9C). This suggests that P-limitation may enhance plant drought tolerance and is mediated via PHR1-PHL1.

DISCUSSION
Owing to the unbiased nature and relatively high sensitivity of the sequencing approach, our RNA-Seq analysis revealed a larger number of novel, previously undescribed P-limitation responsive gene transcripts (PRGTs).
These include novel annotated and unannotated genes and alternative splice forms. The majority of Plimitation induced genes showed reversion after 30 minutes and 3 hours of P addition to P-limited plants, and were regulated by PHR1, PHL1 and PHO2. Our results also showed the involvement of many novel PRGTs in the cross-talk with other stresses including enhanced plant defense and drought tolerance. We have also shown that PHR1-PHL1 mediates the P-limitation induced plant defense and drought tolerance.
We identified more than 1,000 PRGTs that are ≥ 5-fold induced or repressed in shoot and/or root during P-limitation. Of these, ~75% PRGTs were up-regulated and ~25% were repressed, which is similar to 3:1 ratio of induced vs. repressed PRGTs in the ATH1 study of Bustos et al., (2010). Importantly, ~60% of the PRGTs detected in this study were not previously reported in ATH1 experiments, although 26% have an associated ATH1 gene identifier. Some of the differences in the number of PRGTs in different studies could be due to biological variation between experiments (plant, age etc.) or technical differences. For example, difference in P-response to Bustos et al. 2010 experiment could be due to lack of specificity and cross-hybridization of ATH1 probe sets ( Figure S3 and S4) (e.g., 259221_s_at for AT3G03540) with RNA from a highly homologous locus (i.e., AT3G03530). The probe set (259221_s_at) perfectly matches AT3G03530, and hence picks up the signal of AT3G03530, and attributes it to AT3G03540. RNA-Seq however indicates that only AT3G03530 is expressed ( Figure S10), as the two homologous gene transcripts are still sufficiently different (86% identity) for unambiguous read mapping. Experimental variations are also commonly observed among studies (Lee et al., 2020). There might also be genes that were simply not detected by StringTie for example miR399a (AT1G29265), At2g34210, AT4G38080 or AT2G41070 because their transcripts were fused with those from neighboring genes (supplemental Figure S11). . We found that none of these transcript assemblers are perfect, and both miss some transcripts; however, StringTie was better at assembling more transcripts and therefore led to the identification of more PRGTs in this study. AT1G76430/Pht1.9 was not identified by StringTie and is an example showing that also HiSat StringTie doesn't work perfectly. Therefore, one should not completely rely on assembler output and visual inspection in IGV is important for robustness.
Our data showed a link between P-limitation and ion homeostasis. AT2G43500/NLP8, a member of the  et al., 2007;Pant et al., 2015), induction of these UMAMIT genes, mainly in P-limited shoots, suggest their involvement in importing amino acids from xylem which could contribute to higher amino acid accumulation in the leaves during P-limitation.
We identified the genes linking P-limitation to hormone signaling. AT1G68320/MYB62 was found among the P-responsive TFGTs. We identified more than 100 P-limitation inducible unannotated transcripts which were not reported previously and the majority of those were lncNAs. LncRNAs that regulate gene expression by chromatin remodeling or transcriptional interference, by acting as microRNA target mimics s, antisense RNAs, alter the stability and translation of mRNAs (Franco-Zorrilla et al., 2004;Statello et al., 2021). Some of these transcripts can also encode small functional peptides (Choi et al., 2019) and it will be interesting to investigate the function and small peptide coding potential of these lncRNAs. We identified and confirmed that one of these lncRNAs was a microRNA that targets a protein kinase involved in signal transduction (Shiu and Bleecker, 2001) and an We found that the genes related to fatty acid (FAs) metabolism such as AT4G01950 (glycerol-3- We found that PYR/PYL/RCAR was induced by P-limitation suggesting it as a common component for cross-talk between drought and P-limitation. PYR/PYL/RCAR family proteins mediate ABA-dependent regulation of protein phosphatase 2Cs (Fujii et al., 2009;Ma et al., 2009;Park et al., 2009)  Identification of unannotated and annotated novel P-responsive genes involved in regulatory and signaling function, noncoding RNAs, small peptides, novel P-responsive splice forms including miR399 resistant PHO2.2, transporters, genes involved in root hair formation and defense add to our understanding of Psignaling ( Figure 10). Furthermore, the involvement of PRGTs in hormone signaling, drought and defense pathway entails cross-talk between P-limitation with other stress pathways including plant defense. Regulation of the majority of PRGTs reversibly by P-availability and their regulation by PHR1, PHL1 and PHO2 add new potential players to the P-signaling network. Novel insights provided by this study can be useful for translational research to develop genetically enhanced stress resilient crop varieties that can grow in Plimitation conditions. In addition, these results could inform more efficient management strategies.

Supplementary data
Supplemental Figure S1.  Figure S6. Sashimi plots of unannotated highly P-status responsive transcripts.
Supplemental Figure S7. Additional information on annotation and sequence information for the novel PRGTs that are currently unannotated in TAIR10 database.
Supplemental Figure S8. Inter-study comparison of PRGTs in shoots.
Supplemental Figure S9. GO analysis of 44 alternately spliced genes.
Supplemental Figure S10. Sashimi plots for selected strongly P-status responsive gene transcripts.
Supplemental Figure S11. Sashimi plot for miR399a (AT1G29265) as an example of PRGTs missed by RNAseq analysis pipeline.
Supplemental Table S2. Transcription factor genes with P-status responsive transcript abundance as detected in this study using RNA-seq and comparison with previous studies Supplemental   (  2  0  0  8  )  M  i  c  r  o  R  N  A  3  9  9  i  s  a  l  o  n  g  -d  i  s  t  a  n  c  e  s  i  g  n  a  l  f  o  r  t  h  e  r  e  g  u  l  a  t  i  o  n  o  f  p  l  a  n  t  p  h  o  s  p  h  a  t  e  h  o  m  e  o  s  t  a  s  i  s  .  T  h  e  P  l  a  n  t  J  o  u  r  n  a  l   5  3  :   7  3  1  -7  3  8   P  a  n  t  B  D  ,  B  u  r  g  o  s  A  ,  P  a  n  t  P  ,  C  u  a  d  r  o  s  -I  n  o  s  t  r  o  z  a  A  ,  W  i  l  l  m  i  t  z  e  r  L  ,  S  c  h  e  i  b  l  e  W  R   (  2  0  1  5  )  T  h  e  t  r  a  n  s  c  r  i  p  t  i  o  n  f  a  c  t  o  r  P  H  R  1  r  e  g  u  l  a  t  e  s  l  i  p  i  d  r  e  m  o  d  e  l  i  n  g  a  n  d  t  r  i  a  c  y  l  g  l  y  c  e  r  o  l  a  c  c  u  m  u  l  a  t  i  o  n  i  n  A  r  a  b  i  d  o  p  s  i  s  t  h  a  l  i  a  n  a  d  u  r  i  n  g  p  h  o  s  p  h  o  r  u  s  s  t  a  r  v  a  t  i  o (  2  0  0  6  )  P  H  O  2  ,  m  i  c  r  o  R  N  A  3  9  9  ,  a  n  d  P  H  R  1  d  e  f  i  n  e  a  p  h  o  s  p  h  a  t  e  -s  i  g  n  a  l  i  n  g  p  a  t  h  w  a  y  i  n  p  l  a  n  t  s  .  P  l  a  n  t  P  h  y  s  i  o  l  o  g  y   1  4  1  :   9  8  8  -9  9  9   B  u  s  t  o  s  R  ,  C  a  s  t  r  i  l  l  o  G  ,  L  i  n  h  a  r  e  s  F  ,  P  u  g  a  M  I  ,  R  u  b  i  o  V  ,  P  e  r  e  z  -P  e  r  e  z  J  ,  S  o  l  a  n  o  R  ,  L  e  y  v  a  A  ,  P  a  z  -A  r  e Table S1 for more details on these altogether 1,074 gene transcripts.   (A) Transcript length (nucleotides, nt) of the unannotated novel PRGTs identified in present study. Length of these transcripts range from ~1700 nt for the longest and ~150 nt for the shortest with median length 520 nt. Insert shows the number of transcript count (y-axis) vs. number of introns (x-axis) revealing that most of the genes have no introns. (B), (C) and (D) shows Sashimi plots for the novel P-limitation inducible gene transcripts MSTRG.8269, MSTRG.13996, and MSTRG.11040, respectively. Plots were generated using IGV and BAM-files originating from RNA sequencing of the P-limited (-P, red) and P-sufficient (+P, black) shoot samples. Scaling of FPKM read numbers for each transcript is shown in parentheses (upper left corner of each panel) in each panel (1) P-limited shoot; (2) P-sufficient shoot; (3) P-limited root and (4) P-sufficient root.  (A) Sashimi plots and splice junction tracks are shown for PHO2 (At2g33770) in P-sufficient (black) and Plimited (red) shoots and roots. The alternative splicing event that leads to elimination of a large portion of PHO2 including 1 st , 2 nd and 3 rd introns, 2 nd exon of 5' UTR and 1 st exon of coding sequence is shown in pink in the 3 rd and 4 th Sashimi plot. MicroRNA399 binding sites (miR399 BS; green arrow heads) and their consensus sequence are shown in green in the third plot, as well as the presumable alternative translational start in exon 4 (MERSEY…, blue). (B) Agarose gel detection of PHO2 splice variants amplified with primers F1+R or F2+R (indicated in blue in the 4 th plot) and cDNA produced from roots of P-limited (-P) or P-sufficient (FN, +P) Arabidopsis plants. PCR product sizes amplified with F1+R or F2+R are 2460 / 315, or 2286 nucleotides. (C) Schematic representation of PHO2 splice forms (PHO2.1 and PHO2.2) showing RT-qPCR primer locations with arrows and dotted arrows for exon-exon junctions. Thick boxes represent exons, thin boxes represent UTRs and lines represent introns.. (D) Bar diagram showing RT-qPCR expression levels for PHO2.1 and PHO2.2 transcripts in shoot and root in +P and -P conditions. The expression levels are shown on a log 2 scale as 40-ΔCT. The data are average of three replicates ± standard error (SE

Figure 8. Phosphate limitation modulates plant defense via PHR1-PHL1
(A) RT-qPCR expression analysis of plant defense related genes in wild-type (WT), phr1 and phr1phl1 mutants grown under P-limited (-P) condition. Gene expression levels are shown on a Log 2 scale along the y-axis. Gene expression in WT plants grown under P-sufficient condition (WT+P) was taken as a control. The data are average of three replicates ± standard error (SE). (B) Phenotypic response of 4-week-old Arabidopsis thaliana Col-0 grown under short day condition (8:16 light-dark cycle and light intensity 140 µmol m -2 s -1 ) at 5 dpi. Plants were spray inoculated with bacterial suspension of P. syringae pv. tomato DC3000 at concentration of 1X10 8 CFU/ml in water containing 0.01% Silwet L-77. Control plants were inoculated with sterile water containing 0.01% Silwet L-77. (C) Quantification of P. syringae pv. tomato DC3000. (D) Arabidopsis seedlings were flood inoculated with Erwinia carotovora ssp. Carotovora at concentration of 1.6X10 5 CFU/ml in water containing 0.01% Silwet L-77 at 3 dpi. (E) Quantification of Erwinia carotovora ssp. carotovora. (F) Arabidopsis seedlings were flood inoculated with P. syringae pv. maculicola at concentration of 1.6X10 5 CFU/ml in water containing 0.01% Silwet L-77 at 5 dpi (G) Quantification of P. syringae pv. maculicola. Standard error from six biological replicates is shown as error bars in all graphs and experiments were repeated with similar results. Asterisks is shown for statistically significant difference with Student's t test (p value  As described in Figure 2, PRGTs are shown in pink (≥5-fold changed in RNA-seq and without probe set on ATH1 study), blue (≥5-fold changed in RNA-seq and 2-5-fold changed in ATH1 study), red (≥5-fold changed in RNA-seq without support from ATH1 study), yellow (PRGTs of yet unannotated genes), and green (PRGTs ≥ 5fold changed in RNA-seq and ATH1 study. Transcripts labeled with superscripts were previously reported as P-status responsive in the respective reference: (1) Bari et al. (2006), (2)   PRGTs are shown in green (≥5-fold changed in present RNA-seq and ATH1 study), blue (≥5-fold changed in RNA-seq and 2-5-fold in ATH1 study), red (≥5-fold changed in RNA-seq but not changed in ATH1 study) and pink (≥5-fold changed in RNA-seq and no probes set present on ATH1 gene chips).