Genome accessibility dynamics in response to phosphate limitation is controlled by the PHR1 family of transcription factors in Arabidopsis.

As phosphorus is one of the most limiting nutrients in many natural and agricultural ecosystems, plants have evolved strategies that cope with its scarcity. Genetic approaches have facilitated the identification of several molecular elements that regulate the phosphate (Pi) starvation response (PSR) of plants, including the master regulator of the transcriptional response to phosphate starvation PHOSPHATE STARVATION RESPONSE1 (PHR1). However, the chromatin modifications underlying the plant transcriptional response to phosphate scarcity remain largely unknown. Here, we present a detailed analysis of changes in chromatin accessibility during phosphate starvation in Arabidopsis thaliana root cells. Root cells undergo a genome-wide remodeling of chromatin accessibility in response to Pi starvation that is often associated with changes in the transcription of neighboring genes. Analysis of chromatin accessibility in the phr1 phl2 double mutant revealed that the transcription factors PHR1 and PHL2 play a key role in remodeling chromatin accessibility in response to Pi limitation. We also discovered that PHR1 and PHL2 play an important role in determining chromatin accessibility and the associated transcription of many genes under optimal Pi conditions, including genes involved in the PSR. We propose that a set of transcription factors directly activated by PHR1 in Pi-starved root cells trigger a second wave of epigenetic changes required for the transcriptional activation of the complete set of low-Pi-responsive genes.

Phosphorus (P) is an essential nutrient for life (1, 2), playing a fundamental role as a backbone of nucleic acids, in membranes as a component of phospholipids, and by participating in countless energy-dependent metabolic processes (3, 4). Orthophosphate (Pi) availability is a key factor limiting plant growth and yield in many natural and agricultural ecosystems. Although Pi can be present in substantial amounts in the soil, its bioavailability is often severely reduced because of its rapid incorporation into insoluble soil particles due to its high reactivity and its conversion by microbial activity into organic forms not readily taken up by the plant (4, 5). To mitigate Pi limitation, plants have acquired strategies during evolution to better cope with low Pi availability, collectively referred to as the phosphate starvation response (PSR). PSR includes changes in biochemical pathways that reduce Pi requirements, the expression of high-affinity Pi transporters that enhance uptake capacity in the root, the expression of genes encoding enzymes that facilitate uptake of organic sources of Pi, and changes in the root system architecture that increase soil exploration capacity (4). Therefore, the PSR is a complementary set of strategies that allows plants to enhance their capacity to survive and reproduce in soils with low Pi availability (6).

The biochemical and molecular responses to Pi deprivation are relatively well characterized, and several critical components in the underlying signaling transduction pathways have been identified (7–9). The PSR master regulator PHR1 (PHOSPHATE STARVATION RESPONSE1) controls the transcriptional activation of a large set of low-Pi–responsive genes (8, 10). PHR1 regulates transcription via a Pi-dependent interaction with proteins containing SPX domains (7, 11, 12). The Arabidopsis thaliana (Arabidopsis hereafter) genome encodes four highly similar SPX proteins that differ in their subcellular localization and the expression pattern of their encoding genes, suggesting that these proteins have a variety of regulatory roles. The nucleus-localized proteins SPX1 (SPX DOMAIN GENE1) and SPX2 bind to PHR1 in a Pi-dependent manner and prevent its binding to PHR1 binding sites (P1BS motifs) in the promoters of many low-Pi–responsive genes (8). In Pi-limited seedlings, SPX1, SPX2, and SPX3 are degraded via the 26S proteosome pathway, thus allowing PHR1 to bind to P1BS motifs and promote transcription of target genes in response to low Pi, including SPX1. Therefore, the PHR1–SPX1 module regulates Pi starvation responses and establishes a negative-feedback loop for Pi sensing (7, 13, 14). PHR1 is a member of a small subfamily of MYB domain transcription factors that consists of PHR1, PHR1-like 1 (PHL1), PHL2, PHL3, and PHL4. Although PHR1 is a major player in low Pi responses, it is partially redundant with PHL1 and PHL2 (8, 10), as phr1 phl1 and phr1 phl2 double mutant seedlings are less responsive than phr1 seedlings to low Pi (8, 10, 15, 16).

In addition to transcriptional activation by PHR1, global methylation changes have been reported in Arabidopsis and rice (Oryza sativa) in response to low Pi availability (17–19). Changes in methylation patterns were detected in genomic regions harboring phosphate starvation-responsive genes and were often associated with P1BS motifs, suggesting that methylation patterns are related to transcriptional gene regulation in response to Pi limitation (17, 18). The coordination of Pi-induced transcription and epigenetic changes is supported by the important role played in this process by the MEDIATOR complex, which is involved in transcriptional activation and chromatin remodeling (20). Pi limitation also induces differential deposition of the methylation histone marks H3K27me3 and H3K9me2 (21, 22), as well as the deposition of the H2A.Z histone variant (23), suggesting that chromatin remodeling is an integral component of the plant response to low Pi. However, whether the PSR involves changes in chromatin accessibility and whether changes in chromatin accessibility affect the amplitude of the transcriptional response to low Pi remain unknown.

Epigenomic approaches now allow the determination of chromatin alterations associated with developmental disorders and diseases in animals and humans (24–28). Numerous techniques that document genome-wide chromatin states have been developed to investigate the function of chromatin regulators (29–32). One such technique has gained widespread use due to its simplicity and low requirement for input of biological material: assay for transposase-accessible chromatin, followed by sequencing (ATAC-seq) (33, 34). This method has since been employed to characterize and compare the dynamics of chromatin accessibility between cell types and different growth conditions (35–42). In this study, we applied ATAC-seq, together with transcriptome deep sequencing (RNA-seq) to the roots of Arabidopsis wild-type (WT) and phr1 phl2 seedlings to explore the relationship between chromatin accessibility and differential gene expression during Pi limitation, as well as the potential role of the transcriptional regulators PHR1 and PHL2. Our data show that low-Pi conditions induce drastic changes in chromatin accessibility and reveal that PHR1 and PHL2 play key roles in the ensuing chromatin remodeling.

Results

Phosphate Limitation Triggers Changes in Chromatin Accessibility in Arabidopsis Root Cells.

To characterize the chromatin accessibility landscape in response to Pi starvation, we performed ATAC-seq on root cells of Arabidopsis Col-0 (WT) seedlings grown in a hydroponics system (43) under sufficient (1,000 µM, +P) and limiting (10 µM, –P) Pi conditions (). We first validated our experimental conditions by looking for known responses to Pi limitation such as anthocyanin accumulation, reduced shoot growth, increased production of secondary roots, and increased density and length of root hairs (). We collected 10-d-old seedlings; at this time point, we observed a robust establishment of PSR under low-Pi conditions, and it is the same time point at which the role of PHR1 as master regulator was characterized, and PHR1 direct targets were identified (8). We then collected the entire root system and isolated nuclei from two independent sets of seedlings, growing in +P or –P conditions, to prepare ATAC-seq libraries for sequencing (). We obtained an average of 132 million high-quality reads per sample mapping to the Arabidopsis nuclear genome (). Each ATAC-seq replicate showed high correlation, as determined by principal-component analysis (PCA). Importantly, PCA revealed the clear effect of Pi limitation on global chromatin accessibility in Arabidopsis ().

Next, we identified open chromatin regions (peaks) using HOMER (44), resulting in 48,477 high-confidence peaks across all WT samples (+P and –P), of which 14,446 and 14,388 were specific to +P and –P samples, respectively, with the remaining 19,643 peaks being shared (). These results suggested that Pi limitation markedly affected chromatin accessibility in Arabidopsis root cells. Peak calling with HOMER, however, does not provide a sufficiently robust statistical basis to determine differential signals between treatments. We therefore turned to csaw, a sliding-window–based method more effective in determining differential chromatin accessibility regions (DARs) (45–47) (see ). We observed statistically significant differences in chromatin accessibility between +P and –P conditions across all five Arabidopsis chromosomes (Fig. 1 ). An analysis of DARs revealed that of the 6,886 DARs detected, Pi limitation increases chromatin accessibility at 5,712 genomic regions (upDARs) (83%) and decreases chromatin accessibility at 1,174 genomic regions (downDARs) (17%) (Fig. 2 ). The higher proportion of upDARs indicated that chromatin is more accessible in root cells of Pi-limited seedlings compared to Pi-sufficient seedlings. PDLZ2 (PHOSPHOLIPASE D ZETA2) illustrates the typical pattern seen with increased accessibility in response to low Pi over its promoter region (Fig. 1). PDLZ2 encodes phospholipase D Z2 and is up-regulated during PSR (9).

Fig. 1.

Differential accessibility signal in response to phosphate starvation. (A) Heatmap representation of greater and lower ATAC-seq differential chromatin accessibility regions (DARs) in Arabidopsis Col-0 root cells. Each row represents one DAR. The color represents the intensity of chromatin accessibility, from gain (yellow) to loss (dark blue). DARs are grouped based on K-means clustering and aligned to the center of genomic regions. (B) General overview of DARs within a 11-Mb window from a lower arm on chromosome 3 in WT and phr1 phl2 mutant seedlings. Genes are represented by black vertical lines (Top). Red, greater accessibility; blue, lower accessibility. The graphical summary for accessibility was extracted from the Integrated Genome Viewer (IGV) and group-scaled. (C) Example of the PDLZ2 locus, whose chromatin accessibility increases in response to phosphate limitation. Top track: locus organization, with the arrow indicating direction of transcription; red differential signal, greater accessibility; blue differential signal, lower accessibility; red lines, upDARs; blue lines, downDARs; green, +P accessibility profile; and purple, –P accessibility profile. Signal was group-scaled making comparable for the same set of data (+P and –P profiles).

Fig. 2.

Differential accessibility regions (DARs) in response to phosphate starvation. (A) Pie chart representing the total number of upDARs and downDARs. (B) Percentage of upDARs and downDARs as a function of gene region. (C) IPS1 exhibits two upDARs in response to phosphate limitation. (D) PHR1 shows no changes in its chromatin accessibility in response to low Pi. (E) The chromatin at the NR1 locus becomes less accessible in response to Pi limitation. From the Top to Bottom, tracks represent locus organization, with the arrow indicating direction of transcription; red differential signal, greater accessibility; blue differential signal, lower accessibility; red lines, upDARs; blue lines, downDARs; green, +P accessibility profile; and purple, –P accessibility profile. Signal was group-scaled. (F) Gene Ontology (GO) enrichment analysis of genes associated with upDARs in response to phosphate limitation. Values represent adjusted P values.

Since changes in chromatin accessibility are often associated with transcriptional activation or repression of nearby genes (48), we selected the gene closest to each DAR to explore the possible association of DARs with changes in gene expression in response to low Pi. While 62.5% of all upDARs mapped to the promoter regions (1 kb upstream of the transcription start site [TSS]), another 16.9% were located in distal intergenic regions, and 3.8% in 3′-untranslated regions (3′-UTRs). Among downDARs, 62% mapped to promoter regions, 19% to 3′-UTRs, and 5.8% to distal intergenic regions (Fig. 2). We selected representative genes belonging to each category: the chromatin at IPS1 (INDUCED BY PHOSPHATE STARVATION1), a gene whose transcription is known to be highly induced by low Pi (49), became more accessible over its promoter region (Fig. 2); PHR1, whose transcription does not respond to Pi limitation, lacked a significant DAR over the length of its locus (50) (Fig. 2); and NR1 (NITRATE REDUCTASE1), known to be repressed in low-Pi conditions, was associated with a downDAR over a large portion of its coding region (51) (Fig. 2).

To understand the potential influence of changes in chromatin accessibility on transcriptional responses to Pi nutrition, we conducted Gene Ontology (GO) classification to determine enriched categories of genes with upDARs. Notably, we identified no enriched categories related to low-Pi responses for genes associated with upDARs but did observe enrichment for categories related to gene expression, cellular metabolic processes, and nitrogen metabolism (Fig. 2). These results are in line with the known effects of Pi limitation on nitrogen uptake and assimilation (51).

Differential Chromatin Accessibility Is Associated with Changes in Gene Expression.

We next performed RNA-seq experiments to compare transcript levels in WT roots grown in Pi-limited and Pi-sufficient conditions (). We detected 1,012 up-regulated differentially expressed genes (upDEGs) and 1,273 down-regulated genes (downDEGs) in response to Pi limitation (Fig. 3). UpDEGs included most of the classical genes induced by low Pi in Arabidopsis, including PDLZ2, PHT1; 4 (PHOSPHATE TRANSPORTER 1; 4), IPS1 (INDUCED BY PI STARVATION1), AT4/IPS2, SPX1, SPX2, PHF1 (PHOSPHATE TRANSPORTER TRAFFIC FACILITATOR1), and PAP10 (PURPLE ACID PHOSPHATASE10), among others. The most conspicuous gene among downDEGs was PHO2 (PHOSPHATE2), a gene known to be repressed by Pi limitation (52) (Fig. 3).

Fig. 3.

RNA-seq analysis of Col-0 root cells in response to low Pi. (A) MA plot illustrating the number of differentially expressed genes (upDEGs and downDEGs) in the WT. Red dots: up-regulated genes; blue dots: down-regulated genes, as determined by log2(fold change) between Pi-limited and Pi-sufficient conditions. (B) Heatmap representation of normalized gene expression (as Z score) of 28 canonical low-Pi–responsive genes.

We then determined the extent of overlap between DEGs and DARs. Of the 2,285 DEGs, 537 were associated with at least one DAR (Fig. 4). A GO enrichment analysis of these 537 genes yielded categories related to abiotic stress, including “response to phosphate starvation” (Table 1). The set of 537 genes was obtained by comparing all DEGs and all DARs irrespective of their direction of regulation, prompting us to divide this initial gene list into the four possible regulatory outcomes. Out of the 537 genes above, 247 were both up-regulated for their transcript levels (upDEGs) and associated with increased chromatin accessibility (upDARs) in Pi-limited seedlings (Fig. 4). These genes were enriched for GO categories related to cellular response to stress, including “phosphate starvation” and “response to starvation” (Table 1), demonstrating a clear relationship between increased chromatin accessibility and higher transcriptional activity. A second subset derived from the 537 genes consisted of 194 genes (or 15% from 1,273 downDEGs) that exhibit reduced transcript levels and are associated with upDARs in Pi-limited seedlings (Fig. 4); they were enriched in GO categories related to photosynthesis, organonitrogen compound biosynthetic processes, and translation, among others (Table 1). A third subset comprised 63 down-regulated genes (or 4.9% of all downDEGs) associated with downDARs in Pi-limited seedlings (Fig. 4) enriched in GO categories related to plastid functions (Table 1). The final subset representing the overlap between downDARs and upDEGs contained 42 genes (Fig. 4). Association between DARs and DEGs, upDARs and upDEGs, and downDARs and downDEGs was highly significant with P values lower than 1.8 × 10–6. Perhaps as expected the association between upDARs with downDEGs and downDARS with upDEGs had a higher P value; nevertheless, was statistically significant (P > 0.001) (). However, GO analysis of the subset representing the overlap between downDARs and upDEGs did not retrieve any enriched genetic category. Together, these results therefore suggest that the transcriptional response to low Pi is accompanied by changes in chromatin accessibility. We illustrated the behavior of genes from each category with the representative genes PEPC1 (PHOSPHOETHANOLAMINE/PHOSPHOCHOLINE PHOSPHATASE1) (Fig. 4), NADP-MDH (NADP-DEPENDENT MALATE DEHYDROGENASE) (Fig. 4), PHO2 (Fig. 4), and ZAT6 (ZINC FINGER OF ARABIDOPSIS THALIANA6) (Fig. 4).

Fig. 4.

Relationship between DAR-associated genes and DEGs. (A) Venn diagram showing the overlap between DEGs and DARs in WT root cells. (B–E) Venn diagrams of the overlap between upDARs and upDEGs in WT root cells (B), upDARs and downDEGs in WT root cells (C), downDARs and downDEGs in WT root cells (D), and downDARs and upDEGs (E). (F) PEPC1 transcription is strongly increased in response to low Pi and is associated with an upDAR. (G) NADP-MDH transcription is turned off in response to low Pi and is associated with an upDAR. (H) PHO2 transcription is repressed by low Pi and is associated with a downDAR. (I) ZAT6 chromatin accessibility increases specifically in the phr1 phl2 double mutant in low-Pi conditions and is associated with a downDAR. From the Top to Bottom, tracks represent the following: locus organization, with the arrow indicating direction of transcription; red lines, upDARs; blue lines, downDARs; green, +P accessibility profile; purple, –P accessibility; orange, +P RNA-seq; blue, –P RNA-seq. Signal was group-scaled making comparable for the same set of data.

Table 1.

GO enrichment of DEGs associated with DARs

Comparison (gene number)	GO category enriched	P adjusted
DARs and DEGs (537)	Response to abiotic stimulus	5.428 × 10−18
	Cellular response to hypoxia	2.041 × 10−18
	Cellular response to oxygen levels	2.616 × 10−18
	Response to chemical	2.796 × 10−17
	Cellular response to phosphate starvation	4.191 × 10−6
	Cellular response to starvation	8.387 × 10−6
	Cellular response to nutrient levels	1.746 × 10−5
upDARs and upDEGs (247)	Cellular response to stress	4.157 × 10−11
	Cellular response to starvation	8.424 × 10−7
	Cellular response to nutrient levels	2.297 × 10−6
	Response to starvation	1.265 × 10−5
	Cellular response to phosphate starvation	1.464 × 10−5
	Cellular response to extracellular stimulus	3.138 × 10−5
	Response to hormones	2.804 × 10−5
upDARs and downDEGs (194)	Photosynthesis	2.716 × 10−12
	Organonitrogen biosynthetic process	3.942 × 10−6
	Translation	3.012 × 10−5
	Peptide biosynthetic process	3.301 × 10−5
	Chlorophyll metabolic process	3.399 × 10−5
	Chlorophyll biosynthetic process	2.379 × 10−5
	Peptide metabolic process	1.272 × 10−4
downDARs and downDEGs (63)	Chloroplast thylakoid membrane	4.068 × 10−5
	Plastid thylakoid membrane	4.146 × 10−5
	Thylakoid membrane	5.896 × 10−5
	Photosynthetic membrane	6.004 × 10−5
	Chloroplast	8.026 × 10−4
	Chloroplast thylakoid	1.934 × 10−4
	Plastid thylakoid	1.964 × 10−4

GO enrichment of shared genes between DARs and DEGs.

Changes in Chromatin Accessibility in Response to Phosphate Limitation Are Dependent on PHR1 and PHL2.

Based on the observed relationship between DARs and DEGs with cellular responses to Pi limitation, we hypothesized that PHR1 might play an important role in shaping chromatin accessibility in response to this nutritional stress. Accordingly, we first performed RNA-seq on root cells of the Arabidopsis phr1 phl2 double mutant grown in +P or –P conditions. The phr1 phl2 double mutant failed to accumulate anthocyanins when grown in low-Pi conditions, consistent with previously reported observations for phr1 seedlings and further validating our growth conditions (8) (). An analysis of this new RNA-seq dataset identified 269 upDEGs and 137 downDEGs in phr1 phl2 in response to Pi limitation (Fig. 5), which represented a reduction of 74% and 90%, respectively, in DEGs relative to those detected in WT.

Fig. 5.

RNA-seq analysis of phr1 phl2 root cells in response to low Pi. (A) MA plot illustrating the number and the log2(fold change) of upDEGs and downDEGs (differentially expressed genes) in root cells of phr1 phl2. Red dots: up-regulated genes; blue dots: down-regulated genes in response to phosphate limitation. (B) Heatmap representation of the log2(fold change) (logFC) for 1,012 upDEGs in WT and their corresponding logFC in phr1 phl2. Each row represents one gene; red represents higher transcript levels, and blue, lower transcript levels in response to phosphate limitation. (C) Venn diagram showing specific and shared genes between direct PHR1 targets genes, upDEGs detected in WT, and upDEGs detected in phr1 phl2 from our RNA-seq data. (D) Heatmap representation of normalized gene expression (as Z score) of 159 direct PHR1 targets whose transcription was induced in the WT in low-Pi conditions and their corresponding Z score. (E) Heatmap representation of normalized gene expression (as Z score) of 70 direct PHR1 targets whose transcription was induced in the WT and phr1 phl2 in low-Pi conditions.

To independently confirm the weaker transcriptional response of the phr1 phl2 double mutant in low-Pi conditions, we plotted the fold change in estimated transcript levels for the 1,012 genes up-regulated by low Pi in the WT in the form of a heatmap (Fig. 5). Indeed, many genes showed little or no induction by low Pi in the phr1 phl2 double-mutant background. We then examined the expression of 28 markers genes induced by low Pi in different plant species: Their transcriptional activation was also drastically affected in phr1 phl2 specifically during Pi limitation (Fig. 3). Our data corroborate that the phr1 phl2 double mutant is drastically affected in this transcriptional response to Pi limitation.

A set of 2,364 genes was previously identified by chromatin immunoprecipitation sequencing (ChiP-seq) as direct targets of PHR1 in Arabidopsis (53). We thus characterized the changes in gene expression and chromatin accessibility for these PHR1 direct targets. We first determined that of 2,364 genes, only 1,011 showed detectable expression in WT roots under our +P and −P experimental conditions. In addition, 229 of these 1,011 genes were up-regulated in WT roots in response to low Pi, while only 35 were up-regulated specifically in phr1 phl2 but not in the WT (Fig. 5). Of the 229 upDEGs detected in WT roots in low-Pi conditions, 159 appeared to be fully dependent on PHR1 and PHL2, as their transcript levels failed to increase upon exposure to Pi limitation in the phr1 phl2 double mutant. The remaining 70 genes retained responsiveness to Pi limitation in phr1 phl2 (Fig. 5). Moreover, a heatmap representation of the expression changes of these 229 direct PHR1 targets in the WT and phr1 phl2 highlighted the dependency of 159 genes on PHR1–PHL2 (Fig. 5), while most of the 70 PHR1 targets responding to low Pi in the phr1 phl2 double mutant showed a weaker transcriptional induction relative to the WT (Fig. 5). We also discovered a set of 733 genes that are up-regulated in the WT but not the phr1 phl2 double mutant in response to low Pi and that were not previously identified as PHR1 direct targets according to ChIP-seq data (8), indicative of an indirect activation by PHR1 and PHL2. Furthermore, we identified 50 genes induced by Pi limitation that are equally activated in the WT and the phr1 phl2 double mutant and were not previously identified as PHR1 direct targets and therefore may be activated by a PHR1-independent signaling pathway. Finally, 114 genes were up-regulated specifically in the phr1 phl2 double mutant in response to low Pi, suggesting that PHR1 and/or PHL2 repress their transcription (Fig. 5). Interestingly, PHO2 displayed higher transcript levels in +P conditions and an attenuated repression in –P conditions in the phr1 phl2 double mutant compared to the WT (Fig. 3). The severely reduced response exhibited by the phr1 phl2 double mutant confirms that PHR1 and PHL2 play a major role in the direct and indirect regulation of low-Pi–responsive genes in Arabidopsis.

An effect on transcription does not necessarily imply that chromatin accessibility is modulated. To explore the involvement of chromatin accessibility, we performed ATAC-seq on phr1 phl2 root cells () and compared DARs obtained from the WT and the double mutant. This analysis revealed the relative insensitivity to changes in chromatin accessibility in response to Pi limitation in the phr1 phl2 double mutant (Fig. 6), as evidenced by our identification of 1,942 upDARs (34% of WT) and 155 downDARs (13% of WT) in phr1 phl2 (Fig. 6). In addition, only 824 DARs were located in the promoter region of the closest associated gene in phr1 phl2, in contrast to the 3,573 DARs detected for the same region in the WT. The difference between DARs located in distal intergenic regions was less pronounced but nevertheless similar, with 587 DARs for the phr1 phl2 double mutant and 965 for the WT. PHR1 thus appeared to affect chromatin accessibility at the promoter of genes whose transcription is affected by Pi limitation more substantially than that of DARs located at more distal regions.

Fig. 6.

Relationship between DARs and DEGs in phr1 phl2. (A) Heatmap representation of greater and lower ATAC-seq differential chromatin accessibility regions (DARs) in phr1 phl2 root cells. Each row represents one DAR. The color represents the intensity of chromatin accessibility, from gain (yellow) to loss (dark blue). DARs are grouped based on K-means clustering and aligned to the center of genomic regions. (B) Annotation of upDARs as a function of their genomic context. (C) Venn diagram showing the overlap between upDAR-associated genes and upDEGs in phr1 phl2. (D) IPS2 transcription is induced by phosphate limitation in a PHR1- and PHL2-dependent manner. Note the complete loss of DAR in the phr1 phl2 double mutant. From Top to Bottom, tracks represent the following: locus organization, with the arrow indicating direction of transcription; red lines, upDARs; orange lines, downDARs in phr1 phl2; dark green, +P accessibility profile in WT; purple, –P accessibility in WT; light green, +P accessibility profile in phr1 phl2; pink, –P accessibility profile in phr1 phl2; orange, +P RNA-seq in WT; blue, –P RNA-seq in WT; yellow, +P RNA-seq in phr1 phl2; light blue, –P RNA-seq in phr1 phl2. Signal was group-scaled making comparable for the same set of data.

In contrast to the WT, for which upDEGs associated with upDARs were enriched for GO categories related to Pi limitation (Table 1), those for upDEGs associated with upDARs in the phr1 phl2 double mutant were related to responses to chemicals, chemical stimuli, and phytohormones rather than to PSR itself (). This observation suggested that PHR1 and PHL2 play a specific role in PSR but have a limited influence on chromatin remodeling induced by other types of stress indirectly caused by low Pi. We also observed only 52 upDEGs associated with upDARs in the phr1 phl2 double mutant compared to the 247 upDEGs associated with upDARs in the WT (Fig. 6), consistent with the notion that PHR1 and related transcription factors are required for the expression of genes that exhibit greater chromatin accessibility and higher transcript levels. We selected IPS2, one of the most highly responsive genes to low-Pi conditions, to illustrate the effects of the phr1 phl2 double mutant on chromatin accessibility. Chromatin accessibility was not largely affected in the phr1 phl2 double mutant regardless of Pi status and was accompanied by a marked reduction in its transcriptional activation compared to the WT (Fig. 6). These results strongly suggest that PHR1 and PHL2 are not only transcriptional regulators of the low Pi response but also modulate chromatin remodeling in response to Pi limitation.

PHR1-Like Proteins Indirectly Regulate the Chromatin Accessibility Response to Phosphate Limitation in Arabidopsis.

To better understand the functional role of PHR1 and PHL2 in chromatin remodeling in response to Pi limitation, we took a closer look at the 229 direct PHR1 target genes that are up-regulated in low-Pi conditions in the WT and the upDEGs associated with upDARs (hereafter designated upDARDEGs). Of these 229 direct PHR1 targets, 74 were classified as upDARDEGs (Fig. 7 ) in the WT, with another 16 genes behaving as upDARDEGs in both the WT and the phr1 phl2 double mutant (Fig. 7), although their transcriptional activation was lower in the double mutant than in the WT (Fig. 7). Interestingly, we also identified a group of six upDARDEGs in phr1 phl2 but not in the WT, including PPCK1 (PHOSPHOENOLPYRUVATE CARBOXYLASE KINASE1), GDPD1 (GLYCEROPHOSPHODIESTER PHOSPHODIESTERASE1), At1g05000 (PHOSPHOTYROSINE PROTEIN PHOSPHATASE), PPCK2 (PHOSPHOENOLPYRUVATE CARBOXYLASE KINASE2), At3g25240 (SULFATE/THIOSULFATE IMPORT ATP-BINDING PROTEIN), and At3G25795 (Fig. 7), suggesting that PHR1 and PHL2 may reduce chromatin accessibility for these genes in response to Pi limitation. We also identified 156 genes, not previously reported as PHR1 target genes, with higher transcript levels and greater chromatin accessibility in the WT but not in phr1 phl2, effects that we attribute to be under the control of transcription factors and/or chromatin remodelers that are directly activated by PHR1 (Fig. 7).

Fig. 7.

DARs related to the PHR1 and PHL2 transcription factors. (A) Venn diagrams showing the overlap between upDEGs associated with upDARs (upDARDEGs) from WT and phr1 phl2 root cells and direct PHR1 targets. (B) Heatmap representation of the Z score of 16 upDARDEGs showing partial dependency on PHR1 and PHL2, as evidenced by the lower expression in the phr1 phl2 double mutant in low-Pi conditions. (C) Heatmap representation of Z score of 74 upDARDEGs that are fully PHR1- and PHL2-dependent. Each row represents a gene; red represents positive change, and blue, negative change. (D) Average accessibility profile of 74 upDARDEGs that are fully dependent on PHR1 and PHL2. (E) Accessibility profile of 16 upDARDEGs that are partially dependent on PHR1 and PHL2. For D and E, analysis was conducted over a genomic region per gene starting 3 kb upstream of the transcription start site (TSS) and ending 3 kb downstream from the transcription end site (TES). (F) Set of genes responsive to phosphate limitation among the direct PHR1 targets.

The increase in chromatin accessibility seen for the 74 upDARDEGs and PHR1 target genes specific to the WT was completely abolished in the phr1 phl2 double mutant, as evidenced by genome accessibility profiles (Fig. 7). By contrast, the 16 upDARDEGs and PHR1 target genes shared between the WT and phr1 phl2 still exhibited increased chromatin accessibility in phr1 phl2 seedlings, although to a lower extent than the WT (Fig. 7). We noticed a number of genes encoding transcription factors from various families among the 74 upDARDEGs and PHR1 target genes, including NAC048, NAC047, MYB107, MYB34, and WRKY18 (Fig. 7). We hypothesize that these transcription factors initiate a second wave of transcriptional activation in response to Pi limitation and mediate the changes in chromatin accessibility at loci that are not direct PHR1 or PHL2 targets. When subjected to GO enrichment analysis, the set of 74 upDARDEGs, which failed to respond to low Pi in phr1 phl2, returned categories related to stress responses, phosphate-containing metabolic process, and cellular response to phosphate starvation (P value, 1.54 × 10−4) (). When combining these 74 upDARDEGs with the other 16 upDARDEGs shared with the phr1 phl2 double mutant, the GO category related to cellular responses to phosphate starvation became the most enriched with a P value of 6.3 × 10−7 (). Finally, we extended the GO enrichment analysis to the 74 upDARDEGs and direct PHR1 targets, the shared 16 upDARDEGs between the WT and phr1 phl2, and the 156 upDARDEGs only in the WT that are not direct PHR1 targets (Fig. 7). The most enriched category was cellular response to phosphate starvation, with a P value of 6.3 × 10−14, suggesting that the chromatin of PSR-associated genes that are direct PHR1 targets as well as upDARDEGs indirectly activated by PHR1 and PHL2 becomes more accessible, leading to their transcriptional activation upon Pi limitation. We conclude that PHR1 and PHL2 are key players in regulating differential chromatin accessibility and gene expression during the Pi starvation response.

Discussion

Mounting evidence suggests that epigenetic processes, such as changes in DNA methylation patterns and chromatin accessibility, play critical roles in evoking tailored transcriptional programs during development and in response to environmental factors (54). Here, we report that chromatin accessibility substantially contributes to the response of Arabidopsis roots to Pi limitation. Using HOMER, we detected 34,031 and 34,089 peaks of open chromatin in roots exposed to Pi-sufficient and Pi-limited conditions, respectively, numbers that are broadly in line with the 41,419 open chromatin peaks from a previous report assessing genome accessibility in Arabidopsis WT roots using the same bioinformatic approach (38). The differences between our study and the previous study may stem from the experimental design: a hydroponics system (our study) or Petri plates with solid medium (38). Of the 48,477 high-confidence peaks detected in the WT in either growth condition, 14,388 euchromatin regions appeared specifically in Pi-limited conditions.

To perform a comprehensive analysis with statistical support, we turned to csaw, a bioinformatics suite that determines statistically significant changes in differential chromatin accessibility (47), leading to the identification of 5,712 upDARs and 1,174 downDARs in response to low-Pi conditions, with 2,692,770 bp of chromatin becoming accessible and another 1,653,100 bp losing at least partial accessibility. Most DARs (62.5% in +P and 62% in –P) were located within promoter regions, which is comparable to the 58% of DARs located in promoter regions previously reported (38).

We established that 24% of DARs are associated with upDEGs or downDEGs, suggesting that a significant proportion of regions that gain or lose chromatin accessibility experience changes in transcription rates. Although not many studies have been published in plants directly comparing changes in chromatin accessibility with changes in gene expression, ATAC-seq and RNA-seq assays in the parasitic fungus Sparassis latifolia showed that 23% of DARs were associated with upDEGs or downDEGs between control cells and light-induced cells (42). In addition to this first set of DARs, we identified 1,748 DEGs not associated with detectable changes in chromatin accessibility and another 5,509 DARs for which the closest genes did not show differential gene expression, suggesting that changes in gene expression do not necessarily require a change in chromatin accessibility and that changes in chromatin accessibility do not necessarily affect the expression of the closest transcriptional unit. With the former, transcription factors may already be bound to DNA but do not initiate transcription until their cognate signal activates them; with the latter, changes in chromatin accessibility may be related to long-distance transcriptional activation or facilitate other process such as DNA recombination or repair (54, 55). The lack of direct correlation between DEGs and DARs has been previously noted in comparable analyses (56, 57). We acknowledge that associating DARs with the closest transcriptional unit may to some extent bias the analysis, since genes can be regulated by distal regulatory elements (24, 28, 58), but to date there are no methods to directly associate changes in chromatin accessibility and changes in gene expression. Analysis of GO categories of genes associated with low-Pi–induced DARs did not show a clear correlation with PSR genes, confirming that nutritional stress activates a very small set of stress-specific genes and a larger set of general stress-related genes (59). However, GO analysis of DEGs associated with DARs showed a high enrichment for “response to phosphate starvation,” which was more pronounced when running the analysis for upDEGs associated with upDARs that are also direct PHR1 targets.

Only 11.8% of the genes induced by Pi limitation in the WT similarly displayed a transcriptional activation in Pi-limited phr1 phl2 roots (Figs. 3 and 5), which is similar to previous results using the phr1 phl1 double mutant under the similar conditions (8). Among the genes induced by low Pi in the WT, 229 were previously identified as PHR1 targets by ChIP-seq; in agreement, they showed strong enrichment for the GO category “cellular response to phosphate starvation.” Of these 229 genes, the expression of 70 of them was still responsive to Pi limitation in the phr1 phl2 double mutant, albeit to a reduced extent relative to the WT (Figs. 3 and 5). Notably, these 70 genes showed a strong enrichment for the GO category “cellular response to phosphate starvation,” as might be expected since this list of genes included classical Pi limitation genes that are conserved among multiple plant species: genes encoding SPX proteins, purple acid phosphatases, high-affinity Pi transporters, and galactolipid biosynthesis proteins. The fact that these prototypical genes are still partially induced by low Pi in the phr1 phl2 double-mutant background points to a role for additional members of the PHR1 family in the response to Pi limitation, which is only revealed in the absence of the central regulator PHR1.

Several reports have described the relationship between phosphate limitation and other abiotic responses such as nitrogen starvation (1, 51, 59) and drought stress (60, 61). Among the enriched GO categories for the 159 PHR1-dependent DEGs that are direct PHR1 targets (Fig. 5), several were related to drought responses and responses to abscisic acid (ABA) (), confirming that PHR1/PHL2 integrate various stress signaling cascades. Notably, this enrichment was only seen with PHR1-dependent DEGs that are direct PHR1 targets ().

upDARs located in promoter regions were the most affected in phr1 phl2 seedlings in response to low Pi (3,573 in the WT versus 824 in phr1 phl2), whereas upDARs located in distal intergenic regions were less affected (965 in the WT versus 587 in phr1 phl2). Changes in chromatin accessibility can facilitate the binding of transcription factors to the promoters of their target genes, for instance, the binding of PHR1 to the P1BS motif frequently located within 1,000 bp of the transcription start site of PSR genes (8, 62). Notably, none of the direct PHR1 targets whose transcription is induced in response to low Pi encode chromatin remodelers, suggesting that PHR1 and PHL2 recruit the existing chromatin remodeling machinery under Pi-limited conditions to increase chromatin accessibility. This notion is supported by the observation that 90 direct PHR1 genes were up-regulated and gained chromatin accessibility in response to Pi limitation in WT roots. Strikingly, 16 genes out of these 90 were still partially up-regulated and gained chromatin accessibility in the phr1 phl2 double mutant in low-Pi conditions, suggesting that additional members of the PHR1 family can regulate their chromatin accessibility and transcriptional activation. We conclude that the transcription factors PHR1 and PHL2 are required for the changes in chromatin accessibility in response to Pi limitation in Arabidopsis roots. Effects on chromatin remodeling take place largely at promoter regions and to a much lesser extent at distal intergenic (but potentially regulatory) regions (8).

While open chromatin regions provide unobstructed access of transcription factors to their cognate motifs to trigger transcriptional activation or repression, it is more difficult to envision how transcription factors might do so in silent chromatin with the packing of DNA around nucleosomes. However, several transcription factors have been shown to reach their cognate binding sites even when wrapped around a nucleosome by recruiting the chromatin remodeling machinery (63–65). We identified several instances of transcriptional activation in response to low Pi for genes exhibiting a closed chromatin state in Pi-sufficient conditions but gaining chromatin accessibility in response to Pi limitation in the WT (). In the phr1 phl2 double mutant, both chromatin accessibility and transcriptional activation of these genes in response to Pi limitation were lost. These observations suggest that PHR1 may bind to PSB1 sites even within closed chromatin to then recruit the chromatin remodeling machinery, thus acting in a similar fashion as pioneer transcription factors with important roles in cell fate reprogramming.

We also identified two sets of genes relevant to chromatin remodeling that appeared to be down-regulated in the phr1 phl2 double mutant relative to the WT in Pi-sufficient conditions: 1) genes playing an important role in the low-Pi response, such as IPS2, several PAP genes, and three genes encoding high-affinity phosphate transporters; and 2) genes involved in gene expression and chromatin remodeling, such as nine genes encoding Mediator subunits and genes encoding histone acetylases (). That phosphate transporter genes and other phosphate assimilation genes are down-regulated in the phr1 phl2 double mutant when grown in Pi-optimal conditions may provide an explanation for the lower Pi levels measured in phr1 and phr1 phl1 seedlings grown in Pi-sufficient conditions compared to the WT (8). The reduced transcript levels for a set of genes involved in chromosome and chromatin organization in phr1 phl2 seedlings in Pi optimal conditions would, at least in part, explain the drastic loss of the chromatin remodeling response to low Pi in phr1 phl2 roots ().

Several reports have suggested that the response to Pi limitation does not only involve transcriptional activation of genes involved in Pi uptake, transport, and remobilization, but also relies on the remodeling of the photosynthetic apparatus by reducing gene expression and the redistribution of photosynthates, in turn influencing the shoot-to-root biomass ratio (16, 59). In accordance with this, we identified 194 down-regulated genes in Pi-limited WT roots displaying a gain in genome accessibility (Fig. 4) that were mainly related to the regulation of photosynthesis, nitrogen metabolism, and translation (Table 1), suggesting the binding of a transcriptional repressor to their promoters. Another set of 63 down-regulated genes was also associated with downDARs, and were mainly related to mitochondrion and chloroplast function (Table 1). However, their physiological relevance needs to be experimentally validated to determine whether their down-regulation is due to closed chromatin mediated by chromatin remodelers. Although PHR1 is considered a bona-fide transcriptional activator when Pi supply is low (7, 8), we identified six direct PHR1 target genes that are only up-regulated and have increased chromatin accessibility in the phr1 phl2 double mutant, suggesting that PHR1 and PHL2 may directly or indirectly act as transcriptional repressors. Several transcription factors have been described as having dual functions as activators or repressors depending on the cognate DNA motif to which they bind (66). However, it is also possible that PHR1 acts as a repressor by interacting with a negative transcriptional regulator to repress the expression of these genes. Less likely but also possible is that SPX1 may perform an as-yet-undescribed role as a transcriptional activator in the absence of PHR1 and PHL2, although there is no evidence suggesting that SPX proteins bind to DNA. It will be interesting to explore by which mechanisms P-limited plants repress photosynthesis and nitrogen assimilation to develop strategies for breeding plant varieties that can better withstand the PSR.

We showed here that both direct and indirect PHR1 targets gain chromatin accessibility in response to Pi limitation, suggesting that PHR1 and PHL2 are indirectly responsible for driving chromatin accessibility changes triggered by low Pi. In agreement with this possibility, we identified several genes encoding transcription factors including MAF5 (MADS AFFECTING FLOWERING5), NAC048, NAC047, MYB107, ERF1B (ETHYLENE RESPONSE FACTOR 1B), WRKY18, MYB34, and ZAT9, which are all direct PHR1 targets and may in turn activate epigenetic regulators in a second wave of PHR1-independent chromatin remodeling to activate long-term transcriptional responses to Pi limitation. It is also possible that PHR1 modulates chromatin accessibility by altering DNA methylation levels, as it has been previously reported that the expression of several DNA methylates and dimethylases are regulated by PHR1 and harbor PSB1 binding sites in their promoter sequence (17).

Based on our results, we propose a two-wave model to explain the changes in chromatin accessibility observed in Arabidopsis roots subjected to Pi limitation (). In Pi-sufficient conditions, promoter regions are accessible to PHR1, but their local chromatin state does not change in response to Pi limitation; other loci exhibit low chromatin accessibility but gain chromatin accessibility upon binding by PHR1 (). SPX1 is an example of such a direct PHR1 target whose chromatin accessibility does not change in response to Pi limitation. PHR1 triggers transcription of PSR genes with low chromatin accessibility in low-Pi conditions, which gain chromatin accessibility once PHR1 binds to P1BS sites, possibly by recruiting the chromatin remodeling machinery. A typical example is PDLZ2 (). We also propose that, in response to Pi limitation, PHR1 indirectly triggers a second wave of changes in chromatin accessibility by activating the transcription of other transcription factor–encoding genes or genes encoding enzyme that modify DNA methylation levels, whose concerted action may enhance gene methylation, histone modification, and chromatin remodeling to regulate the expression of genes within closed chromatin regions that are not direct PHR1 targets. This class of genes represents the largest subset of genes whose expression is activated by low Pi and gains chromatin accessibility but is not among direct PHR1 targets (). Further research to decipher the exact cascade of regulation is necessary and may lead to the design of new breeding strategies to enhance Pi nutrition in crops.

Methods

Plant Materials and Growth Conditions.

WT Arabidopsis (Arabidopsis thaliana) accession Columbia-0 (Col-0) was used for all experiments; the phr1 phl2 double mutant in the Col-0 background was obtained from Dong Liu (Tsinghua University, Beijing, China). Seeds were surface sterilized with 90% (vol/vol) ethanol for 5 min and 50% (vol/vol) bleach solution for 5 min before four washes with sterile distilled water. Seedlings were grown using a hydroponic system with 0.1× Murashige and Skoog (MS) medium (43) supplemented with one of two phosphate concentrations: +P (1,150 µM) and –P (10 µM) using KH2PO4. All genotypes were analyzed 10 DAG (days after germination). Seedlings were grown at 20 °C in an 18-h light/6-h-dark photoperiod. To test the phenotypes of the double mutant, seedlings were grown for 20 DAG on MS medium containing 5 µM phosphate before they were photographed ().

Isolation of Nuclei.

For nuclei isolation, we followed the method of Bajic et al. (36) with some modifications: The extraction buffer consisted of 15 mM Tris⋅HCl, pH 7.5, 20 mM NaCl, 80 mM KCl, 0.2% Triton X-100, and 5 mM β-mercaptoethanol. Nuclei suspensions were filtered using a 30-µm mesh; to eliminate mitochondria and chloroplasts, we used sucrose sedimentation buffer and centrifuged at 600 rpm for 20 min: 20 mM Tris⋅HCl, pH 8, 0.2 mM MgCl2, 2 mM EDTA, 0.2% Triton X-100, 15 mM β-mercaptoethanol, and 1.7 M sucrose.

ATAC Library Preparation and Sequencing.

For ATAC-seq assays, two replicates per sample were processed using between 80,000 and 120,000 nuclei, as determined by flow cytometry. Chromatin was digested by Tn5-mediated tagmentation and adapter incorporation, according to the manufacturer’s protocol (Nextera DNA sample preparation kit; Illumina) at 37 °C for 30 min; each library was amplified for 12–15 cycles according to the published protocol (37). The quality of the libraries was assessed by a DNA-based fluorometric assay and by electrophoresis. Samples were sequenced on a HiSeq2500 Illumina sequencer system as paired-end reads of 2 × 50 bp.

RNA Extraction.

Harvested root tips were frozen in liquid nitrogen and ground to a fine powder. Total RNA was isolated using TRIzol reagent (Invitrogen) according to the manufacturer’s instructions. mRNA-seq libraries were generated by Novogene.

RNA-Seq Analysis.

Adapter sequences and low-quality reads were removed from raw reads with trimGalore v0.6.4 (https://github.com/FelixKrueger/TrimGalore). Mapping of reads to the genome and gene counts were performed using RNA-STAR v2.7.5b (67) and Galaxy (68) through the usegalaxy.eu server, and read counts over genes were obtained using htseq-count v0.9.1+galaxy1 (69). Differential gene expression analysis was performed using the edgeR package in R (70). Analysis of GO enrichment and cluster analysis by biological process were performed using g:profiler (71). Heatmaps of differentially expressed genes were constructed following published bioinformatics methods (72).

ATAC-Seq Bioinformatic Analysis.

Trimming of adapter sequences and removal of low-quality reads from raw reads were performed using trimGalore v0.6.4 (https://github.com/FelixKrueger/TrimGalore). Clean reads were then aligned to the Arabidopsis TAIR10 release 43 reference genome using Bowtie2 v2.3.5.1 (73) with options -k 10 –very-sensitive. PCR duplicates were marked with sambamba-markdup v0.7.0 (74); all steps up to this point in the analysis were automated using snakePipes (75). PCR duplicates and reads mapping to the organellar genomes were removed with samtools v1.10 (76). Quality control of filtered mapped data were performed using ATACseqQC v1.10.4 (77).

Peaks were called for each replicate using the find Peaks function within HOMER suite v4.11 (44) with the following parameters: -style histone -size 75 -minDist 75 and -gsize 1.2e8. The resulting peak files were merged by experimental condition with HOMER merge Peaks function. Regions with differential accessibility were estimated using csaw (46). Mapped reads were counted genome-wide in 75-bp windows, and only windows with a Log2(signal enrichment) > 1 relative to background were considered for further steps. For differential accessibility estimation, replicates were normalized using the trimmed mean of M values (TMM) method, and adjacent windows (up to 150 bp apart) with differential signal between conditions were merged up to a maximum size of 5,000 bp. The resulting regions were filtered for significance using absolute log2(fold change) ≥ 0.8 relative to the control condition and with a false-discovery rate < 0.05. Peaks and DARs datasets were annotated to the TSS of the nearest gene using ChipSeeker v1.22.1 with org.At.tair.db and TxDb.Athaliana.BioMart.plantsmart28 Bioconductor packages (78–80). Promoters were defined as spanning 1,000 bp of sequence upstream of the TSS and 400 bp of sequence downstream of the TSS. Analysis of GO categories and cluster analysis by biological process were performed using g:profiler (71). Signal visualization files and images were generated using deepTools v.3.5.0 (81). MultiBamSummary scaling factors were used to generate bigwig files with bamCoverage. Overlap between genomic regions was determined using Intervene 0.6.4 (82).