Low-phosphate chromatin dynamics predict a cell wall remodeling network in rice shoots

One-sentence summary: Combined chromatin structural data reveals specific chromatin-state transitions that correlate with subsets of functionally distinct rice genes differentially expressed under phosphate starvation. hydrolase; Integ, cell wall integrity protein; IQ CaM, IQ calmodulin-binding motif protein; JAZ, 777 ZIM domain-containing JAZ protein; kCoA, 3-ketoacyl-CoA synthase; kinase; KNOT, knotted- 778 1-like homeobox protein; LIG, lignin dirigent; lipase; LRX, leucine-rich repeat extensin; LTPL, 779 Protease inhibitor/seed storage/LTP protein; LYM, lysM domain-containing GPI-anchored 780 protein; MEE, maternal effect embryo arrest; MYB; N-Gly, shiga/ricin-like N-glycosidase; NPC, 781 non-specific phospholipase; OXI, oxidoreductase; P450, cytochrome P450; PDD, PD-(D/E)XK 782 nuclease superfamily protein; PEC, pectinase; PLA, phospholipase A; PLATZ; PNM, 783 phosphoethanolamine N-methyltransferase; POX, peroxidase; PPR, pentatricopeptide repeat 784 protein; RALF, Rapid ALkalinization Factor; RLK, receptor-like kinase; S1P, Subtilisin Site-1 785 Protease; SAM, S-adenosyl-L-methionine-dependent methyltransferases; SRO, OsSRO1c; 786 SULF, sulfotransferase; UDP, UDP-glucuronosyl/UDP-glucosyltransferase; UVB, ultraviolet-B- 787 repressible protein; VQ, VQ domain containing protein; WAK, wall-associated kinase; WAX, WRKY; ZGT, ZGT circadian clock coupling factor; ZR1, FYVE zinc finger 789 domain protein. synthase; kinase; KNOT, knotted-1-like homeobox protein; LIG, lignin dirigent; lipase; LRX, leucine-rich repeat extensin; LTPL, Protease inhibitor/seed storage/LTP protein; LYM, lysM domain-containing GPI-anchored protein; MEE, effect arrest; MYB; N-Gly, shiga/ricin-like N-glycosidase; NPC, non-specific phospholipase; OXI, oxidoreductase; P450, cytochrome P450; PDD, PD-(D/E)XK nuclease superfamily protein; PEC, pectinase; PLA, phospholipase A; PLATZ; PNM, phosphoethanolamine N-methyltransferase; POX, peroxidase; PPR, pentatricopeptide repeat protein; RALF, Rapid ALkalinization Factor; RLK, receptor-like kinase; S1P, Subtilisin Site-1 Protease; SAM, S-adenosyl-L-methionine-dependent methyltransferases; OsSRO1c; SULF, sulfotransferase; UDP, UDP-glucuronosyl/UDP-glucosyltransferase; UVB, ultraviolet-B-repressible wall-associated WAX,


Introduction
6 TEG/TE is larger than that of H3K4me3. To further examine the apparent association between 139 H3K4me3 and H2A.Z, we first computed a correlation coefficient using deepTools (Ramírez et 140 al., 2016), which showed that both chromatin marks were correlated across the rice genome (r = 141 0.77; Pearson correlation coefficient; Figure 1C). Next we identified and compared distinct 142 H3K4me3 and H2A.Z peaks using SICER (Zang et al., 2009) and BEDTools (Quinlan and Hall,143 2010). This identified 32,886 H3K4me3 peaks and 44,804 H2A.Z peaks, of which 30,813 (93% 144 of H3K4me3 peaks) overlapped ( Figure 1D). Finally, we plotted the average abundance of both 145 marks centered at peak summits of H3K4me3 (Supplemental Figure S2A) or H2A.Z 146 (Supplemental Figure S2B). Together, these analyses demonstrate substantial co-occurrence of 147 H3K4me3 and H2A.Z in the rice genome.  Table S1). PCG were ranked according to FPKM and divided into five expression 152 quintiles, as well as a sixth group of genes that were not expressed (i.e. FPKM = 0). We found a 153 clear, positive correlation between transcript abundance and H3K4me3 localization around the 154 TSS ( Figure 2A,B), consistent with studies from a variety of species (Bernstein et al., 2002, 155 Santos- Rosa et al., 2002, Barski et al., 2007, Zhang et al., 2009, Van Dijk et al., 2010. In 156 contrast, transcript abundance exhibited a general negative correlation with TES-and gene body-7 To evaluate a potential role for H3K4me3 in modulating Pi-deficiency responses, we 168 carried out H3K4me3 ChIP-seq on shoots from plants subjected to a 24-hour Pi-deficiency 169 treatment (Supplemental Table S1). As shown in Supplemental Figure S3A,  H3K4me3 distribution at PCG, such that the prominent 5′ peak was reduced. Separating PCG 171 into two broad functional categories showed that house-keeping genes exhibited the 5′ peak 172 reduction (Supplemental Figure S3B), whereas stress-responsive genes exhibited an increase in 173 H3K4me3 at the TSS (Supplemental Figure S3C). These data along with our prior studies 174 (Zahraeifard et al., 2018 indicate that nucleosome occupancy, H2A.Z, and 175 H3K4me3 each exhibit distinct changes in response to Pi starvation.  Figure 3A). Because the enrichment values are relative, a low score for a particular mark does 188 not indicate a lack of the mark. CS1 and CS2 were each deficient in both H2A. Z and H3K4me3, 189 CS3 was enriched in only H2A.Z, CS4 was enriched in both H2A. Z and H3K4me3,and CS5 was 190 enriched in only H3K4me3. Regarding nucleosome density, CS2 and CS3 had moderately higher 191 nucleosome enrichment compared to the other three states. To support the presence of these five 192 particular chromatin states in the rice genome, we repeated the ChromHMM analysis by 193 combining our three datasets from control samples with publically available datasets for two 194 marks recognized as largely repressive, DNA methylation (Secco et al., 2015) and/or H3K27me3 195 (Zhang et al., 2012) from control, shoot samples. As shown in Supplemental Figure S4, the 196 addition of these marks did not change the combinatorial effect of the three marks we used to 8 designate the five chromatin states, verifying the relevance of these states in our original 198 ChromHMM analysis. The inclusion of DNA methylation or H3K27me3 alone each yielded a 199 sixth chromatin state, whereas including both repressive marks yielded a sixth and seventh state 200 (Supplemental Figure S4A). In addition, a genome comparison of the five and seven chromatin-201 state analyses showed similar percentages of the genome represented by each CS except CS1. 202 When the repressive marks were added, some regions of the genome defined as CS1 were 203 changed to CS6 or CS7 (Supplemental Figure S4B,C). 204 Next we mapped the distribution of the five chromatin states across the genome (divided 205 into 200-bp bins), which revealed several biases with genomic features ( Figure 3B,D). CS1 was 206 the major chromatin state, accounting for 63% of the rice genome, and was enriched at TE and 207 TEG. It should be noted that highly repetitive regions of the genome were likely designated CS1 208 due to low numbers of mappable reads rather than bona fide depletion of the chromatin marks 209 examined. TE and TEG were also enriched in CS2 and CS5. This means that the transposable 210 element-related loci were either deficient in both H2A.Z and H3K4me3 or enriched in H3K4me3 211 only. In contrast, PSG were enriched in CS2 and CS3, consistent with depletion of both H2A.Z 212 and H3K4me3 or enrichment of only H2A.Z. Finally, PCG were enriched in CS4, consistent 213 with enrichment of both H3K4me3 and H2A.Z. To more specifically characterize PCG, we 214 calculated enrichments at the TSS and TES separately ( Figure 3B). Compared to all bins within 215 PCG, the TSS was more enriched in CS4, CS5, and CS3, whereas the TES was more enriched in 216 CS3 and less enriched in CS4. These results indicate for PCGs generally an overall high 217 occupancy of H2A.Z and/or H3K4me3 at the TSS, but an enrichment of only H2A.Z at the TES.

218
Pi starvation has a dramatic impact on chromatin signatures 219 To characterize the impact of Pi starvation on chromatin signatures we compared the 220 distribution of chromatin states between control and Pi-deficiency conditions. First we analyzed 221 the enrichment of each chromatin state within the four gene types ( Figure 3C,D). In response to 222 Pi starvation, the genomic bins within TE, TEG, and PSG increased in CS1, but decreased in 223 CS2, CS4, and CS5, consistent with a loss of H3K4me3. On the other hand, PCG bins decreased 224 in CS1, but increased in CS2, CS3, and CS5. At the TSS of PCG, CS2 and CS4 decreased and 225 CS1 increased, whereas at the TES, CS1 and CS4 decreased and CS2 increased (p-value < 0.01) 226 ( Figure 3B,C, Supplemental Figure S5). Because the TSS is an important regulatory region for 227 9 gene expression, we investigated the effect of Pi starvation at the TSS of PCG in more detail.  Figure 4, the enriched GO terms for CS4-CS3 genes fell 240 into five GOMCL clusters, including transcription factor activity, response to endogenous 241 stimulus, cell wall, oxygen binding, and response to extracellular stimulus. In contrast, CS4-CS5 242 genes were enriched in GO terms defined by nine GOMCL clusters, which among other 243 functional categories, were related to translation and gene expression, nuclear functions, plastid 244 functions, nucleic acid metabolism, development, and RNA binding. Interestingly, CS5-CS1 245 genes shared essentially the same enriched GO terms (Supplemental Dataset 3) and GOMCL 246 clusters ( Figure 4) as CS4-CS5 genes. One explanation for this is that the CS5-CS1 and CS4-247 CS5 transitions are frequently found together at the TSS. Indeed, examination of the bins that 248 flank the TSS (Supplemental Figure S8) showed that CS5-CS1 genes were approximately four 249 times more likely than random to exhibit a CS4-CS5 transition in the bin downstream of the TSS 250 (p-value < 0.001). Similarly, the CS4-CS5 genes were 3.5-fold more likely to contain a CS5-CS1 251 transition upstream of the TSS (p-value < 0.001). In contrast, CS5-CS1 genes with a CS4-CS5 252 upstream bin, and CS4-CS5 genes with a CS5-CS1 downstream bin were similar to random or 253 under-represented, respectively. Thus, the identification of subgroups of functionally similar 254 genes with CS5-CS1 and CS4-CS5 transitions at their TSS is reflective of these genes containing 255 a specific pair of transitions (CS5-CS1 + CS4-CS5, 5′-3′) at the TSS. We analyzed our recent RNA-seq experiments (Zahraeifard et al., 2018;Supplemental 259 Table S1) to investigate the relationship between gene expression and chromatin state transitions 260 in response to Pi starvation. Differential expression analysis with DESeq2 identified 1385 261 differentially-expressed genes (DEGs) in response to Pi starvation, 694 up-regulated and 691 262 down-regulated (adjusted P-value < 0.001; Supplemental Figure S9, Supplemental Dataset 4).

263
GO terms enriched for up-regulated genes included response to stress, lipid metabolic process, 264 and signal transduction, whereas down-regulated genes were enriched in growth, cell-cell 265 signaling, and lipid, carbohydrate, and secondary metabolic processes (Supplemental Table S2, 266 Supplemental Figure S9B,C). Although lipid metabolism was overrepresented in both groups of 267 DEGs, genes linked to carotenoid biosynthesis and alpha-Linolenic acid metabolism were among 268 the up-regulated DEGs, whereas cutin, suberin, and wax biosynthesis were among the down-269 regulated DEGs. Overall, the functional categories of these DEGs were similar to those from 270 previous transcriptome studies of Pi-deficient plants (Thibaud et al., 2010, Cai et al., 2013, Secco 271 et al., 2013, Zahraeifard et al., 2018.

272
Next we investigated whether differential expression in response to Pi starvation 273 correlated with distinct CS transitions. We quantified the overlap between the up-or down-  Figure S8), at the TSS are unlikely to be differentially 294 expressed after 24-hour Pi deficiency. Because a number of translation-related genes were 295 previously shown to be down-regulated by long-term (21-day) Pi deficiency in rice shoots 296 (Secco et al., 2013), we carried out a bootstrapping analysis to test whether our CS5-CS1 and 297 CS4-CS5 genes were enriched among those DEGs. Indeed, both CS5-CS1 and CS4-CS5 groups 298 were enriched (p-value < 0.01) among genes down-regulated by long-term Pi deficiency 299 (Supplemental Figure S10). This suggests that the chromatin dynamics observed at these genes 300 after 24 hours of Pi starvation is a prelude to decreased transcript abundance not observable until 301 after a longer duration of Pi deficiency.

302
In addition to the correlations between DEGs and distinct chromatin state transitions, 303 there were also correlations between DEGs and groups of genes that did not transition ( Figure 5).

429
Our data support this by revealing that more than 40% of all rice PCG in shoots exhibit a 430 chromatin state transition at their TSS in response to a 24-hour Pi-deficiency treatment, and that 16 several specific transitions correlate with subgroups of genes differentially expressed by Pi 432 starvation. It is noteworthy that although we found significant correlations between chromatin 433 dynamics and subsets of differentially expressed genes, there were many genes exhibiting a 434 chromatin state transition that were not differentially expressed (Figure 4). This lack of a strong 435 global correlation between chromatin dynamics and changes in gene expression has been 436 reported previously (Zong et al., 2013, Fiziev et al., 2017, and likely reflects regulation of also showed high similarity to our CS4-CS3 DEG profile, suggesting the apparent apoplastic 524 signaling network overlaps with multiple stressors. Our CS4-CS3 DEG list contains many 525 orthologs of Arabidopsis components involved in salinity stress responses, including FER, LRX, 526 and RALF peptides (Zhao et al., 2018). It is of interest to evaluate whether the rice orthologs 527 exhibit similar functions in response to stressors including salinity and Pi starvation.  Also, the Spp1-independent genes tend to be more highly expressed during control conditions,

20
and repressed during environmental stress, whereas the Spp1-dependent genes generally exhibit 554 low expression during control conditions and induced expression during stress. Finally, in 555 response to diamide stress, many yeast ribosomal protein genes are down-regulated and exhibit a 556 decrease in H3K4me3 (Weiner et al., 2015). Our data suggest that rice employs different 557 mechanisms to modulate H3K4me3 levels at distinct gene groups, similar to yeast. This is Sterilization and pre-germination (1 day at 37°C followed by 2 days at 28°C) were 567 carried out on rice cultivar Nipponbare (Oryza sativa ssp. japonica) seeds. Seeds were 568 transferred to 12-h light/12-h dark at 30°C/22°C conditions to germinate for 14 days. Seedlings 569 were grown hydroponically in modified Yoshida Rice culture media as described (Yoshida et al., 570 1971, Secco et al., 2013. The solution was replaced every 7 days. After 21 days, seedlings were 571 used for a 24-hour Pi-deficiency treatment (modified Yoshida Rice solution without NaH 2 PO 4 ).      Table S1. Summary of ChIP-seq and RNA-seq libraries (short reads).