All Arabidopsis (Arabidopsis thaliana) plants used in this study were in the Colombia-0 (Col-0) ecotype background. The T-DNA insertional lines SALK_067629 (phr1) and SAIL_731_B09 (phl1) were obtained from the Arabidopsis Biological Resource Center (ABRC). These two mutants have been proved to be the null mutants (Nilsson et al., 2007; Bustos et al., 2010). The null allele of PHL4 gene, phl4-2 (referred as phl4 in this work), was generated by CRISPR/Cas9-based genome editing system (Wang et al., 2018). phr1phl1, phr1phl4, phl1phl4, and phr1phl1phl4 were produced through genetic crossing using above mentioned single mutants. PHL4 OX-1 (referred as PHL4 OX in this work) was generated as described in Wang et al., 2018.

Surface-sterilized seeds were sown on Pi-sufficient (+Pi) or Pi-deficient (-Pi) medium. The +Pi medium contained half-strength Murashige & Skoog basal salts (Caisson Labs, MSP01-01190008), half-strength Murashige & Skoog vitamin powder (1000×) (Phyto Technology Laboratories, M533), 1.0% (w/v) sucrose, 0.5% MES, and 1.2% (w/v) agar (Sigma-Aldrich, catalog no. A1296) or 0.8% (w/v) agarose (Biowest, catalog no.111860). The pH of the medium was adjusted to 5.8 with NaOH. In the -Pi medium, half-strength Murashige & Skoog without Pi (Caisson Labs, MSP11-05160009) was used to replace Murashige & Skoog basal salts. After seeds were stratified for 2 days at 4°C, the agar plates were placed vertically in a phytotron with a photoperiod of 16 h light and 8 h dark at 22-24°C. The light intensity was 100 μmol m^-2 s^-1. Nicotiana benthamiana plants were grown in soil under the same lighting conditions.

Cellular Pi content and total P content were quantified according to Wang et al. (2011). Specifically, the weighed fresh shoot and root tissues were submerged in 1 ml of 1% glacial acetate and were alternatively freeze-thawed 10 times in liquid nitrogen and in a water bath maintained at 65°C. A 100-μl volume of the extract was mixed with 200 μl of ddH₂O and 700 μl of Pi reaction buffer containing a mixture of [0.48% NH₄MoO₄, 2.85% (v/v) H₂SO₄] and [10% (w/v) ascorbic acid] in a ratio of 6:1. The reaction was allowed to proceed at 37°C for 1 h. The cellular Pi content was determined at A₈₂₀ according to a prepared standard curve and was expressed as μmol/g FW.

To determine the total P content, about 50 mg of fresh tissues was oven-dried at 500°C for 3 h and flamed to ash. The ashes were dissolved in 100 μl of 30% (v/v) HCl and 10% (v/v) HNO3. About 10 μl of the dissolved sample was mixed with 290 μl of ddH₂O and 700 µl of Pi reaction buffer, and the Pi was quantified by above method. The total P contents of plant tissues were determined and expressed as Pi contents extracted from flamed ashes.

Anthocyanins in shoots were extracted with propanol: HCl: H₂O (18: 1: 81, v/v/v) in dark at room temperature for 24 h. Absorbance was measured at 530 and 650 nm. Anthocyanin content was expressed as (A₅₃₀ – A₆₅₀)/g FW.

Total RNAs of 8-day-old seedlings were extracted using the Magen HiPure Plant RNA Mini Kit. A 2 μg quantity of the RNAs were reversely transcribed to cDNA using M-MLV reverse transcriptase (Takara). qRT-PCR analyses were carried out using 2×RealStar Green Fast Mixture (GenStar, catalog no. A301-10) on a Bio-Rad CFX96 real-time PCR detection system. The ACTIN 2 gene (At3g18780) was used as an internal control, and the relative expression level of each gene was calculated by the 2^−ΔΔCt method (Livak and Schmittgen, 2001). The primers used for RT-qPCR analyses are listed in Supplementary Table 1.

For LCI assays, the CDSs of PHR1, PHL1, and PHL4 were inserted into the vectors pCAMBIA-nLUC and pCAMBIA-cLUC (Chen et al., 2008), respectively. For BiFC assay, the CDSs of PHR1, PHL1, and PHL4 were individually cloned into the vector nYFP or cYFP using a one-step isothermal in vitro recombination procedure (Gibson et al., 2009). The resultant constructs were mobilized into the Agrobacterium tumefaciens strain GV3101 and used to transform the leaves of Nicotiana benthamiana. The LCI and BiFC assays were performed as described by Sun et al. (2016). The primers used for vector construction are listed in Supplementary Table 1.

The yeast two-hybrid assays were performed using the Matchmaker GAL4 Two-Hybrid System (Clontech) according to the manufacturer’s instructions. The CDSs of PHR1 and PHL1 were cloned into pGBKT7 vector, while the CDS of PHL4 was cloned into pGADT7 vector. The various combinations of the two constructs were co-transformed into the yeast strain AH109. Transformants were selected on SD/-Trp/-Leu/-His medium and SD/-Trp/-Leu/-His/-Ade medium.

The RNAs of the root of 8-day-old Col, phr1, phl1, phl4, phr1phl1, phr1phl4, phr1phl1phl4, and PHL4 OX grown on +Pi and –Pi medium were extracted and sent to Anoroad genome company for RNA-sequencing analysis. Three biological replicates were used in this experiment. Sequencing was performed on an Illumina platform generating 150bp pair-end reads. The expression levels of genes were normalized by TPM (transcripts per million). Genes with a fold change of |Log₂FC|≥1 and FDR<0.05 were considered as differentially expressed. Comparisons were made between -Pi and +Pi wild-type plants to identify PSR genes (Supplementary Table 2) as well as Pi-deficient mutant and wild-type plants to identify mutant-affected genes (Supplementary Table 4). Heatmap was drawn using Morpheus (https://software.broadinstitute.org/morpheus/). GO enrichment was performed using the Agrigo v2 (http://systemsbiology.cau.edu.cn/agriGOv2/).

To determine the functional relationship among PHR1, PHL1, and PHL4 in regulating plant development and responses to Pi starvation, we further generated Arabidopsis phl1phl4 double mutant and phr1phl1phl4 triple mutant through genetic cross (for details, see Materials and Methods). We first compared their growth phenotypes with phr1, phl1, phl4, phr1phl1, and phr1phl4 in soils. The seeds of the WT and all the mutants were germinated on 1/2 MS medium and grown for seven days before they were transferred to soils. Previously, we reported that none of the single mutant of phr1, phl1, or phl4 showed obvious developmental defects in soils, but phr1phl1 showed early senescence phenotype (Wang et al., 2018, 2023). This work also showed that the one-month-old soil-grown plants of phr1phl1 and phr1phl1phl4, but not of phr1phl4, displayed early senescence phenotype (yellowish leaves) (Figure 1B). When the plants became 45-days old, phr1phl1phl4 had more senescent leaves and higher degree of senescence than those of phr1phl1 (Figure 1C). phr1phl1phl4 also flowered later than the WT and phr1phl1. These results suggested that PHL1 and PHL4 not only act redundantly with PHR1 in regulating leaf senescence, but also in regulating flowering time.

Next, we analyzed developmental and physiological responses of the mutants to Pi starvation. The seeds of the WT and all the mutants were directly germinated on 1/2 MS media with different amount of Pi supplementations and primary root (PR) growth was evaluated eight days after germination (8 DAG). When the Pi concentration was reduced from 625 μM Pi (normal 1/2 MS medium, or +Pi medium) to 300 μM, the seedlings of all the mutants displayed similar PR length with that of the WT (Figures 2A, G). When Pi concentration was reduced to 50 μM, the PR length of all the mutants still did not differ from that of the WT, except phr1phl1phl4 whose PR length was only about 70% of the WT (Figures 2C, I). When grown on the medium with 25 μM Pi, the PR length of phr1, but not phl1 and phl4, was shorter than the WT (Figures 2D, J). On this medium, the PR growth of two double mutants and the triple mutant also became progressively more inhibited than phr1. On the medium with 10 μM Pi or without Pi supplementation (-Pi medium), the PR lengths of phr1, phr1phl1, and phr1phl4 were similar, but that of phr1phl1phl4 was further reduced (Figures 2E, F, K, L).

Induction of anthocyanin production in shoots is a hallmark response of plants to Pi starvation. We then analyzed the anthocyanin contents in the shoots of these mutants. The seeds of plants with different genotypes were sown on +Pi and -Pi media and grown for 12 days with strong light intensity (200 μmol m^-2 s^-1) before the analyses were performed. On +Pi medium, all types of seedlings accumulated very low levels of anthocyanins in their shoots (Supplementary Figure 1). Under -Pi condition, WT seedlings accumulated a large quantity of anthocyanins in their shoots. The anthocyanin content of phr1, but not of phl1 and phl4, was reduced to about 20% of the WT. Both phr1phl1 and phr1phl4 displayed further reduction of anthocyanin content with similar extents compared with phr1. The triple mutant phr1phl1phl4 did not show further decrease of anthocyanin content compared with phr1phl1 (Supplementary Figure 1).

The above results together indicated that PHR1 plays a major role, and PHL1 and PHL4 play equally minor role in regulating Pi starvation-induced inhibition of PR growth and accumulation of anthocyanins.

Finally, we investigated the role of PHR1, PHL1, and PHL4 in regulating Pi homeostasis by analyzing the Pi contents in all mutants at 10 DAG (Figure 3A). Under +Pi condition, the Pi content of phr1 decreased to about 60% of the WT in shoots. The Pi content of phr1phl1, but not phr1phl4 displayed further decrease compared to that of phr1. phl1 and phl1phl4 also had a small but significant reduction of Pi content compared to that of the WT. In +Pi shoots, phr1phl1phl4 did not show further decrease of Pi content compared to phr1phl1. In Pi sufficient-roots and Pi deficient-shoots, Pi contents of all mutants were similar to that of the WT. However, in Pi-deficient roots, the Pi content of phr1phl1 was almost doubled of that of the WT, consistent with what we reported before (Wang et al., 2018, 2023). The Pi content of phr1phl1phl4 were further increased compared to that of phr1phl1 (Figure 3A). We also measured the total P contents in all mutants. The total P contents in all samples were closely correlated with their cellular Pi levels (Figure 3B). These results suggested that PHR1, PHL1, and PHL4 redundantly regulated Pi homeostasis, but both PHL1 and PHL4 only played a minor role. And, among them, PHL4 had the least contribution.

Previously, PHR1 has been shown to directly interact with PHL1 and PHL4 (Bustos et al., 2010; Wang et al., 2023). To determine whether PHL4 could also directly interact with PHL1, we first performed luciferase (LUC) complementation imaging assays. The coding sequences (CDS) of PHL1 and PHL4 were fused to the C-terminal half of the luciferase (cLUC) and the CDS of PHR1 was fused to the N-terminal half of the LUC (nLUC). PHL1-cLUC and PHL4-cLUC were co-transformed with PHR1-nLUC, respectively, into the leaves of N. benthamiana with proper controls (nLUC or cLUC). All the gene constructs mentioned in this study were under the control of 35S CaMV promoter. Like PHL1-cLUC and PHR1-nLUC, the co-expression of PHL4-cLUC and PHR1-nLUC resulted in a strong fluorescence signal from reconstituted LUC activity (Figure 4A), consistent with the previous reports (Bustos et al., 2010; Wang et al., 2023). A PHL4-nLUC construct was also made. When this PHL4-nLUC fusion gene was co-expressed with PHL1-cLUC or PHL4-cLUC in the leaves of N. benthamiana, they also generated strong reconstituted LUC signals. These results demonstrated that PHL4 not only could interact with PHL1 but also with itself. The bimolecular fluorescence complementation assays performed in the leaves of N. benthamiana supported above conclusions (Figure 4B).

We further carried out the yeast two-hybrid experiments to determine the relative strength of the interactions among these transcription factors. The CDS of PHR1 and PHL1 were individually fused to the DNA-binding domain of the transcription factor GAL4 (PHR1-BD and PHL1-BD) while the CDS of PHL4 was fused to its transcriptional activation domain (PHL4-AD). The yeast cells carrying PHR1-BD and PHL4-AD or PHL1-BD and PHL4-BD constructs were grown on the selection media. The results not only confirmed the interactions between PHL4 and PHL1, and PHL4 and PHR1, but also indicated that the interaction between PHL4 and PHL1 was stronger than that between PHL4 and PHR1 (Figure 4C).

To determine the individual roles of PHR1, PHL1, and PHL4 in regulating plant transcriptional responses to Pi starvation at genomic levels, we conducted transcriptomic analyses of the roots of WT, phr1, phl1, phl4, phr1phl1, phr1phl4, phr1phl1phl4 and PHL4-overexpressing lines. Total RNAs were extracted from the roots of 8-day-old seedlings grown under +Pi and -Pi conditions. We first performed RT-qPCR analyses of six hallmark PSI genes for all batches of RNAs. These six PSI genes included a non-coding transcript, IPS1 (Burleigh and Harrison, 1999); a microRNA, miR399D (Fujii et al., 2005); two high-affinity phosphate transporters, AtPT1 (Pht1;1) and AtPT2 (Phtl;4) (Muchhal et al., 1996); a ribonuclease, RNS1 (Bariola et al., 1994); and an acid phosphatase, ACP5 (AtPAP17) (del Pozo et al., 1999). The results showed that the expression of all six genes were significantly induced in the WT (Supplementary Figure 2), indicating that the experimental condition was set up properly. These batches of RNAs were then subjected to sequencing using next generation sequencing technology. Three biological replicates were used for all genotypes. The RT-qPCR analyses of the same six PSI genes for all the mutants were also closely correlated with the results from RNA-seq data, validating the approaches of RNA-seq experiments (Supplementary Figure 2).

Using |Log₂FC| ≥ 1 (differential expression level greater than 2 folds) and FDR (false discovery rate) < 0.05 as cutoff, we identified 2579 PSR genes which included 1339 PSI genes and 1240 Pi starvation-suppressed (PSS) genes for the WT (Supplementary Table 2). In this work, we focused our analyses only on the PSI genes for all genotypes. GO analysis showed that most PSI genes identified in the WT were involved in glucosinolate biosynthetic process (GO:0019761), Pi transport (GO:0006817), cellular response to Pi starvation (GO:0016036) (Supplementary Table 2). The remaining PSI genes were involved in responses to abiotic and biotic stresses such as hypoxia, wounding, salt, water deprivation, pathogens as well as in some hormone signaling pathways, including jasmonic acid, karrikin, abscisic acid. These results were largely consistent with what we reported recently (Wang et al., 2023) and with those reported by other research groups in the past (Wu et al., 2003; Misson et al., 2005; Morcuende et al., 2007; Bustos et al., 2010; Osorio et al., 2019; Hani et al., 2021).

To obtain a global view of the Pi starvation-induced changes of transcriptomes for the WT and all the mutants, we drew a heatmap of PSI genes using the Log₂FC value of Pi-deficient WT and Pi-deficient mutants versus Pi-sufficient WT. Our visual examination revealed that under Pi deficiency, the transcriptome of phr1 displayed significant changes compared to that of the WT, whereas those of phl1and phl4 only showed a slight change (Figure 5A; Supplementary Table 3). Both phr1phl1 and phr1phl4 double mutants showed further alterations of the transcriptomes compared to that of phr1, but the alterations in phr1phl4 was lower than those in phr1phl1. phr1phl1phl4 showed the further changes of the expression profile compared to those of two double mutants. The heatmap therefore indicates that PHR1, PHL1, and PHL4 all were involved in regulating PSI gene expression, but PHL1 and PHL4 alone only played minor roles, and between these two, PHL4 contributed the least.

Next, we investigated the individual roles of PHR1, PHL1, and PHL4 in regulating the expression of PSI genes by quantitative analyses. To determine how the mutation of these three PHR family members affected gene expression under Pi starvation, we defined those genes in the mutants whose expression was significantly higher than those in the WT as “PSI-up genes” and those genes whose expression was significantly lower than those in the WT as “PSI-down genes”. We also arbitrarily used a 2-folds difference of expression levels between the WT and a mutant with FDR < 0.05 to define a given PSI gene as “significantly mutation-affected”. With these criteria, we found that the expression of the majority of mutation-affected genes was downregulated compared to those in the WT (Figure 5B; Supplementary Table 4). This result indicated that the major functions of PHR1, PHL1, and PHL4 were the promotion of gene transcription under Pi starvation. In phr1, the expression of total 393 PSI genes was affected; among them, 343 were downregulated and 50 were upregulated; whereas in phl1 and phl4, only 39 and 31 PSI genes whose expression were affected, respectively, i.e. 35 downregulated and 4 upregulated in phl1 and 29 downregulated and 2 upregulated in phl4 (Figure 5B). The number of mutation-affected genes was increased in two double mutants, but the increase in phr1phl1 (620) was larger than the increase in phr1phl4 (462). The number of mutation-affected genes in phr1phl1phl4 was further increased to 686. These results demonstrated there were synergistic interactions among PHR1, PHL1, and PHL4 in regulating the transcription of PSI genes.

All the PSI genes were further divided into three categories according to their responsiveness to Pi starvation (High: Log₂FC > 5; Medium:2 < Log₂FC ≤ 5; Low: 1 ≤ Log₂FC ≤ 2). Our analyses indicated that 79.4% of the PSI genes of high Pi responsiveness (Log₂FC > 5) were downregulated in phr1 and this percentage decreased to 13.7% in phl1 and 6.9% in phl4, respectively (Figures 5B, C). In the double and triple mutants, this percentage did not significantly change compared to that of phr1 (83.3% in phr1phl1, 79.4% in phr1phl4 and 83.3% in phr1phl1phl4). In terms of medium-responsive PSI genes, the percentages of PSI-down genes were decreased to 33% in phr1, 3.8% in phl1, and 3.6% in phl4. This percentage was significantly increased to 49.0% in phr1phl1, 43.6% in phr1phl4, and 53.3% in phr1phl1phl4. For the low-responsive PSI genes, the percentages for the PSI-down was further decreased to 15.6% in phr1, 0.7% in phl1, and 0.9% in phl4, respectively. Compared to the single mutants, this percentage was increased to 28.3 in phr1phl1, 19.0% in phr1phl4, and 33.4% in phr1phl1phl4. Taken together, these results suggested that when combined with PHR1, PHL1 and PHL4 mutations did not significantly increased the number of high-responsive PSI-down genes, but increased the number of medium- and low-responsive PSI-down genes. Moreover, compared to PHL1, PHL4 has weaker enhancing effect when it acted with PHR1.

Considering that the mutation of PHL1 and PHL4 in phr1 background did not significantly increase the number of mutation-affected high-responsive PSI-down genes, we wondered whether the combined mutations might more affect the magnitude of PSI gene transcription. We then quantified the difference of average expression level (DAE) of PSI-down genes between the WT and each mutant. For a given mutant, its DAE was determined by calculating the difference of the expression level of each differentially expressed gene between the -Pi WT and -Pi mutant, then taking the average value of these differences (expressed as Log₂FC) (Figure 6A). The DAE of PSI-down genes for phr1 was highest among all single mutants in all type of PSI-down genes, especially for high-responsive ones. The DAE was increased in both phr1phl1 and phr1phl4, but the increase in phr1phl1 was more significant than that in phr1phl4. Compared to phr1phl1, the DAE was further significantly increased in phr1phl1phl4 for high-responsive PSI-down genes (6.00 vs 4.82), but less significantly increased for medium-responsive PSI-down gens (2.58 vs 2.38), and not significantly increased for low-responsive PSI-down genes (1.86 vs 1.85). When such analyses were performed using all PSI genes, including both PSI-up and PSI-down genes, the results displayed the similar trends (Supplementary Figure 3). This was probably because the PSI-up genes only account for a small portion in all the mutants.

Furthermore, we selected a group of PSI genes that are known to be involved in various Pi starvation responses, such as Pi signaling (SPX proteins, IPS1, miRNA399s), Pi transport (Pi transporters), Pi remobilization (purple acid phosphatases, PAP), and lipid remodeling (the enzymes involved in the generation of glycolipids and sulfate lipids to replace phospholipids in biomembranes). The transcriptomic changes of these genes were compared among the WT and each mutant under Pi starvation. For the nine members of Pht1 family (high affinity Pi transporters, hereafter, were collectively referred as Pht1 genes), we analyzed all of them except for Pht1:6 that is specifically expressed in young floral bud and pollen (Nussaume et al., 2011). Mutation of PHR1 resulted in more than 50% decrease in all Pht1 genes except for Pht1;1 (18%) and Pht1;4 (45%), whereas mutation of PHL1 and PHL4 only reduced the expression of Pht1;5 to less than 50% of the WT (Figure 6B). For the double mutants, both phr1phl1 and phr1phl4 had further reduction of the expression of Pht1 genes compared with phr1. Overall, the extent of the reduction of Pht1 gene expression in phr1phl1 was greater than that in phr1phl4. For the triple mutant, the expression levels of the Pht1 genes were further decreased; in which the expression of six of Pht1 genes was reduced to less than 5% of that of the WT. The regulatory effects of PHR1, PHL1, and PHL4 on the expression of the PSI genes that were involved in Pi signaling, Pi remobilization, and lipid remodeling were similar with that on Pht1 genes (Figure 6B).

We have previously reported that overexpression of PHL4 increased plant Pi starvation responses including enhanced root-associated APase activity, elevated cellular Pi and total P contents, and increased expression of typical PSI genes (Wang et al., 2018). To further understand the working mechanism of PHL4, we compared the expression profile between the WT and PHL4 OX line by RNA-seq analysis under normal growth condition. Compared to the WT, there were 869 genes whose expression was upregulated and 1213 genes whose expression was downregulated in the PHL4 OX line (Supplementary Table 5). We then examined whether PHL4 functions by affecting the expression of other PHR genes. The results showed that in the PHL4-overexpressing (PHL4 OX) line, the expression of PHL4 was increased around 500 folds while the expression of other PHR genes were not obviously altered (Figure 7A).

Finally, we conducted GO analysis of differentially expressed genes at genomic scale in the PHL4 OX line under normal growth condition. For those PHL4 OX-upregulated genes, they were enriched in cellular responses to Pi starvation (GO:0016036), Pi transport (GO:0006817), including all Pht1 genes, except for Pht1:6, and glucosinolate biosynthetic process (0019761) (Figure 7B). For the PHL4 OX-downregulated genes, they were enriched in response to abiotic stress (water deprivation, salt stress), hormone-related process (jasmonic acid, abscisic acid, salicylic acid), cell wall-related process (cell wall modification, lignin metabolic process, pectin catabolic process), and Fe-related process (intracellular sequestering of iron ion, cellular response to iron ion) (Figure 7C). These PHL4 OX-affected genes were largely consistent with the PSR genes observed in Pi-deficient plants (Wang et al., 2023). Therefore, these results provided further evidence for the function of PHL4 in regulating transcriptional responses to Pi starvation.