Deposition of Histone Variant H2A.Z within Gene Bodies Regulates Responsive Genes

Summary

The regulation of eukaryotic chromatin relies on interactions between many epigenetic factors, including histone modifications, DNA methylation, and the incorporation of histone variants. H2A.Z, one of the most conserved but enigmatic histone variants that is enriched at the transcriptional start sites of genes, has been implicated in a variety of chromosomal processes. Recently, we reported a genome-wide anticorrelation between H2A.Z and DNA methylation, an epigenetic hallmark of heterochromatin that has also been found in the bodies of active genes in plants and animals. Here, we investigate the basis of this anticorrelation using a novel h2a.z loss-of-function line in Arabidopsis thaliana. Through genome-wide bisulfite sequencing, we demonstrate that loss of H2A.Z in Arabidopsis has only a minor effect on the level or profile of DNA methylation in genes, and we propose that the global anticorrelation between DNA methylation and H2A.Z is primarily caused by the exclusion of H2A.Z from methylated DNA. RNA sequencing and genomic mapping of H2A.Z show that H2A.Z enrichment across gene bodies, rather than at the TSS, is correlated with lower transcription levels and higher measures of gene responsiveness. Loss of H2A.Z causes misregulation of many genes that are disproportionately associated with response to environmental and developmental stimuli. We propose that H2A.Z deposition in gene bodies promotes variability in levels and patterns of gene expression, and that a major function of genic DNA methylation is to exclude H2A.Z from constitutively expressed genes.

Introduction

In addition to packaging the DNA to fit within the cell, histones function to control the structure and accessibility of the chromatin environment by altering the biochemical properties of the nucleosome or through the recruitment of distinct binding partners. These actions promote changes in transcription that regulate the proper timing of developmental decisions and appropriate responses to the external environment. One such method of histone-mediated control comes from the exchange of the canonical histones with non-allelic histone variants, which alter the fundamental structure and stability of the nucleosome [1]–[4].

H2A.Z is one of the most enigmatic of these histone variants, as well as the most well-conserved, with a single origin at the root of eukaryotic evolution [1]. H2A.Z has been implicated in a number of apparently disparate and even contrary chromosomal processes, including heterochromatic silencing, gene activation, transcriptional memory, cell-cycle progression and thermal-sensory response [5]–[9]. A common aspect of H2A.Z biology is its enrichment within the few nucleosomes surrounding transcription start sites (TSS), which has been demonstrated by genome-wide localization experiments in protozoa, fungi, animals, and plants [10]–[19].

The conserved H2A.Z distribution pattern at the TSS in many species has lead to considerable effort to understand the effect of H2A.Z on transcription. H2A.Z enrichment at promoters in yeast is simultaneously required for transcription and inversely correlated with transcription level [11], [20], [21]. Studies in animals have reported that H2A.Z exhibits a positive correlation with transcription [12], [22], [23], although some have found that this relationship is only true up to a point, after which the association becomes negative [15], [24]. In plants, the relationship between H2A.Z at the TSS and transcription appears to be roughly parabolic, with the highest and lowest expressed genes having the least H2A.Z enrichment [16].

In the yeasts Saccharomyces cerevisiae and Schizosaccharomyces pombe, H2A.Z regulates genes that respond to changes in the environment [21], [25], [26], and loss-of-function mutants fail to react appropriately to external cues [27], [28]. Arabidopsis thaliana plants lacking PIE1 (AT3G12810), the plant homolog of the SWR1 catalytic subunit of protein complexes responsible for the deposition of H2A.Z in yeast and mammals [29]–[36], exhibit misregulation of many genes involved in the innate immune response [37]. Recent work has shown that Arabidopsis plants with a mutated ARP6, which encodes a component of the PIE1 complex, inappropriately express temperature response genes, leading to the proposal that H2A.Z may act specifically as a thermosensor in plants [8].

The genomic distribution and biological functions of DNA methylation, another well-conserved feature of chromatin, are in many aspects strikingly different from those of H2A.Z. DNA methylation in the form of 5-methylcytosine is present in all vertebrates examined to date, as well as in many invertebrates, fungi, and plants [24], [38], [39]. The primary function of eukaryotic DNA methylation has long been considered to be the silencing of the sequences it decorates, particularly transposable elements [40], although the recent discovery of gene body methylation in plants and animals, the functional significance of which is still unknown, has complicated this view [24], [38], [41]–[45]. Whereas H2A.Z is enriched near the TSS of most genes, TSS-proximate DNA methylation is strongly associated with transcriptional repression in plants and vertebrates [46].

Recently, we reported a strong, genome-wide anticorrelation between H2A.Z and DNA methylation in Arabidopsis, including in bodies of active genes [16]. Results from similar studies in vertebrates suggest that this anticorrelation is a conserved feature of eukaryotes [24], [47], [48]. In Arabidopsis, we showed that changes in DNA methylation caused by a mutation in the DNA methyltransferase MET1 induced reciprocal alterations in H2A.Z deposition, demonstrating that DNA methylation antagonizes H2A.Z recruitment [16]. We also used a null mutation in PIE1 (pie1–5) to examine the effect of disrupted H2A.Z function on DNA methylation. By coupling methylated DNA immunoprecipitation (MeDIP) to microarray analysis, we found a low magnitude but genome-wide DNA methylation increase in genes that suggested a mutual antagonism between H2A.Z and DNA methylation [16].

There is now considerable evidence that the PIE1 complex deposits H2A.Z into chromatin in Arabidopsis, though whether it has H2A.Z-independent functions, as has been shown for other eukaryotic SWR1 homologs, remains unclear [29], [49], [50]. Other Arabidopsis chromatin remodelers are probably also able to deposit H2A.Z, as does the yeast INO80 complex [51], because H2A.Z is incorporated into nucleosomes at low levels in pie1 and swr1 mutants [29], [34], [52]. Given that the sets of genes that are misregulated in H2A.Z and SWR1-related mutants only partially overlap in both S. cerevisiae and Arabidopsis [29], [37], we sought to use an H2A.Z-deficient plant line, as opposed to SWR1-related mutants, for further analysis of H2A.Z function.

Here, we describe the characterization of an H2A.Z loss-of-function line in Arabidopsis thaliana. We find that loss of H2A.Z in Arabidopsis has little effect on the level or profile of DNA methylation in genes, and propose that the global anticorrelation between DNA methylation and H2A.Z is primarily caused by the exclusion of H2A.Z from methylated DNA. We show that the level of H2A.Z enrichment in gene bodies is generally correlated with gene responsiveness and that lack of H2A.Z causes misregulation of genes that respond to a variety of stimuli. We propose that H2A.Z deposition in gene bodies promotes gene responsiveness, but may destabilize constitutive expression, and that a major function of gene body DNA methylation is to exclude H2A.Z from constitutively expressed genes.

Results

Construction of a near-null Arabidopsis h2a.z mutant line

Three of the thirteen Arabidopsis H2A genes, HTA8 (AT2G38810), HTA9 (AT1G52740), and HTA11 (AT3G54560), have been classified as encoding H2A.Z based on phylogenetic analyses [53], [54], and distribution patterns and genetic studies suggests that these proteins are largely functionally redundant [16], [33], [34]. Recently published work has demonstrated that a double mutant of hta9-1 and hta11-1 produced plants with phenotypes similar to those found in null pie1–5 mutants [37]. To generate a line devoid of H2A.Z, we crossed hta9-1 and hta11-1 plants with a line bearing an insertion in HTA8, hta8-1 (Figure 1A). Contrasting with recent evidence that individual knockouts of the two vertebrate H2A.Z isoforms exhibit different phenotypes [55], we did not observe morphological abnormalities in any of the three single mutant lines. The resulting triple mutant line, which we will refer to as h2a.z, is both viable and phenotypically distinguishable from WT (Figure 1B). Transcripts of HTA8 and HTA11 were not detectable in the h2a.z mutant by RT-PCR, but low levels of HTA9 RNA were present (∼26% of wild-type; Figure 1C and 1D) in h2a.z plants but not in hta9-1 single mutants, suggesting that the intronic T-DNA insertion in HTA9 is spliced out in a fraction of transcripts, as confirmed by sequencing of the cDNA (Figure S1). To test whether this low level of expression was the result of a genetic rearrangement at the HTA9 locus that occurred in our crosses, we recreated the h2a.z line using hta9-1 plants lacking HTA9 transcript (Figure 1D). The h2a.z progeny from the independent cross produced similar phenotypes to the original h2a.z line and similar RT-PCR results for HTA9, suggesting upregulation of HTA9 in the triple mutant.

A fourth gene, HTA4 (AT4G13570), is the closest H2A family member to the three H2A.Z genes and has been categorized as H2A.Z-like [54], but all publically available data indicate that HTA4 is not expressed at significant levels in any WT tissue. To ensure that HTA4 is not upregulated as a result of the drop in H2A.Z levels in our h2a.z line, we tested the expression of HTA4 by RT-PCR (Figure 1E), and did not detect HTA4 RNA in h2a.z or in WT. Taken together, our data indicate that the h2a.z line has less than ten percent of wild-type H2A.Z transcript levels. Despite reduced fertility (Figure 1F), h2a.z plants are viable and produce offspring, differing markedly from the lethality of strong H2A.Z mutations in other multicellular organisms [15], [56], [57], [58], [59].

The h2a.z mutant phenotype is distinct from that caused by lack of PIE1

We measured the number of leaves present when the plant produced its first flower buds in h2a.z and WT (Figure 2A). In short days (SD), the h2a.z line flowered significantly earlier than WT, with 23.2+/−1.1 leaves vs. 49.7+/−1.5 leaves (P-value<0.0001, two sample T-test). In long days (LD), the difference in flowering time between h2a.z and WT was less pronounced, with 8.3+/−0.2 leaves and 10.6+/−0.2 leaves, respectively (P-value<0.0001), but the difference in rosette size and plant stature was greater in LD than SD (Figure S2). Of the first ten flowers, 22+/−3.1% in LD and 76+/−4.6% in SD exhibited extra petals (between 5 and 8) in the h2a.z mutant line, compared to 1.5+/−0.6% (LD) and 2+/−0.8% (SD) in WT (Figure 2B and 2C). The h2a.z mutant also exhibited short, thickened siliques, a phenomenon potentially related to decreased fertility. The h2a.z siliques averaged 4.8+/−0.1 and 5.6+/−0.1 mm in length under LD and SD conditions, compared to 10.6+/−0.1 and 11+/−0.2 mm for WT (Figure 2D–2F). The h2a.z phenotypes described above, as well as increased leaf serration and petiole length in SD (Figure 2G), are similar to those previously published for hta9-1; hta11-1 and pie1–5 mutants [32], [37], [60].

The h2a.z line exhibited several phenotypes not previously reported for pie1–5 or hta9-1;hta11-1. First, while both pie1–5 and h2a.z have reduced stature, pie1–5 plants are more severely dwarfed and have a bushy, extensively branched architecture, whereas h2a.z plants are spindly and have trouble remaining upright (Figure 2H). This aspect of the h2a.z phenotype might be caused by contributions from the WS ecotype of hta8-1 (all other lines are in the Col ecotype), but this is unlikely because the WT siblings from the same cross do not show these traits. Additionally, many of the siliques in the h2a.z mutant exhibited a strong asymmetric curvature, most likely due to the improper fusion of its carpels (Figure 2F). Other novel phenotypes occurred only rarely, affecting multiple aerial plant tissues including leaf and stem structures, but were most prevalent among floral organs (Figure S3). The most striking examples were the inappropriate emergence of petals and stamens directly from the stem, and flowers with improperly fused carpels, leading to severely compromised reproductive structures.

A mutation in yeast swr1 (pie1) ameliorates many of the phenotypes observed with the htz1 (h2a.z) single mutant, as well as the severe phenotypes of the double mutant between htz1 and set3 [61], [62]. The cause of the htz1 phenotypes was proposed to be chromatin disruption by the SWR1 complex in the absence of its proper substrate, a hypothesis supported by SWR1-dependent accumulation of DNA damage in the absence of htz1. To test whether simultaneous removal of the PIE1 chromatin remodeler and H2A.Z would reduce the severity of phenotypes seen in h2a.z plants, we crossed the h2a.z mutant line to pie1–5. Contrary to the results from yeast, the phenotype of the Arabidopsis double mutant is more severe than that of either parent – progeny exhibit early developmental arrest, dying shortly after germination (Figure 2I). Taken together with the phenotypic disparity, our results suggest that H2A.Z and PIE1 have non-redundant functions in Arabidopsis. Because h2a.z is not a complete loss of function line, the stronger phenotype of h2a.z; pie1 plants might be caused by a further reduction of H2A.Z incorporation into chromatin, but nevertheless demonstrates that pie1–5 does not entirely abolish H2A.Z function. While we cannot rule out the possibility that all pie1 phenotypes are associated with H2A.Z, we consider this unlikely because of the stronger effect of pie1 on plant architecture compared to h2a.z (Figure 2H).

Lack of H2A.Z does not substantially perturb genic DNA methylation

To test our hypothesis that H2A.Z protects genes from DNA methylation, we generated genome-wide methylation profiles for the h2a.z mutant and WT using shotgun bisulfite sequencing. Because plants have DNA methylation in three different sequence contexts, CG, CHG, and CHH (H = A, T or C), which are largely controlled by distinct families of methyltransferases and have different genome-wide distributions [24], [38], it is advantageous to use an assay that has single base-resolution to distinguish between these contexts. Two biological replicates each of h2a.z and WT were generated for each of three different tissue types that represent different stages along a developmental continuum: 14 day-old whole seedlings, 6 week-old rosette leaves, and 6 week-old cauline leaves. One biological replicate was taken from the original h2a.z mutant line, and the second from the additional h2a.z line generated from independent crosses with the same T-DNA insertional alleles. Analysis of the average methylation levels across all genes revealed that a loss of H2A.Z in Arabidopsis has little effect on the global patterns of DNA methylation in CG, CHG or CHH contexts (Figure 3A and Figure S4). For comparison, we generated bisulfite sequencing data for two biological replicates each of pie1 and sibling WT seedlings, and one replicate of h2a.z;pie1 seedlings. As with the results for the h2a.z mutant, the pie1 and h2a.z;pie1 mutants showed only subtle changes compared with WT in the global patterns of genic DNA methylation (Figure S5). Previously, we used locus-specific bisulfite sequencing to validate our microarray results at five candidate genes scored as hypermethylated in the pie1 mutant; all five showed modest but consistent gains in CG methylation [16]. Similar analyses performed here on both our h2a.z and pie1 mutants demonstrate an overall consistency between current and previous pie1 data. They also show that the hypermethylation at these loci in pie1 is less consistently present in h2a.z (Figure S6). Statistical analyses of the methylation differences between the h2a.z, pie1, and h2a.z;pie1 mutants and their respective WT controls suggest that there are subtle increases in genic methylation as a result of H2A.Z loss (Figure S7). These results are consistent with our published data, showing a small but statistically significant increase of genic methylation in the pie1 mutant [16]. However, the quantitative data generated here demonstrate that this increase is of a very low magnitude and is unlikely to substantially contribute to the global anticorrelation between H2A.Z and DNA methylation.

Unexpectedly, the h2a.z mutant exhibited tissue-specific DNA methylation changes in transposable elements (TEs; Figure 3B and Figure S8). CG methylation was marginally increased over wild-type in four of the six replicates, with the most consistent change in seedlings, whereas CHG methylation decreased more heavily in the older tissues, though there is considerable variation between replicates (Figure 3B and Figure S8). CHH methylation was substantially reduced specifically in cauline leaves (Figure 3B and Figure S8). Kernel density estimations of these changes demonstrate that the majority of transposons show a modest change in methylation, rather than a larger effect in a small subset of TEs (Figure S9). Analyses of DNA methylation in pie1 and h2a.z;pie1 seedlings show that, like h2a.z seedlings, these lines exhibit increased CG methylation in TEs (Figure S10). Curiously, the h2a.z;pie1 seedlings exhibit decreases in CHG and CHH TE methylation that are not seen in seedlings of pie1 or h2a.z, but which are reminiscent of the decreases in h2a.z plants later during development (cauline and rosette leaves; Figure 3B–3C and Figure S10). Our data indicate that whereas a loss of H2A.Z does not substantially change DNA methylation within genes, lack of H2A.Z affects TE methylation in all three sequence contexts in a development-specific manner.

We hypothesized that a larger effect of H2A.Z on DNA methylation may be undetectable in our h2a.z mutant due to the carefully targeted nature of routine DNA methylation maintenance. By this logic, H2A.Z's role may be to protect against random and spurious accumulation of DNA methylation at the TSS over evolutionary timescales, rather than to act as a barrier to regular DNA methylation processes. To test this hypothesis, we performed crosses of h2a.z and pie1 to two methylation mutants, ibm1–6 and met1–6, in which normal methylation targeting to genes is perturbed. We expected that if H2A.Z were acting to prevent methylation from accumulating at the TSS, there would be greater increases in methylation in these double mutants than in the parental lines. IBM1 (AT3G07610) encodes a H3 lysine 9 demethylase, MET1 (AT5G49160) encodes the primary CG DNA methyltransferase, and both ibm1 and met1 mutations cause increased CHG methylation in gene bodies [43], [44], [63], [64], [65], [66], [67]. Single mutant plants are viable and fertile (Figure 3D), but h2a.z;ibm1, h2a.z;met1, pie1;ibm1, and pie1;met1 double mutants die shortly after germination and exhibit severe developmental abnormalities, including the production of undifferentiated callus-like material, under-sized root systems, and premature flowering (Figure 3E).

Bisulfite sequencing of h2a.z;ibm1, h2a.z;met1, and pie1;ibm1 seedlings revealed that a loss of H2A.Z does not strongly alter the genic methylation profile in any context from that seen in the parental backgrounds (Figures S11, S12, S13). Statistical analyses of the CG methylation differences between the h2a.z;ibm1, pie1;ibm1, and the ibm1 control suggest that the loss of H2A.Z in these double mutant lines leads to a subtle increase in genic methylation, as we found in the h2a.z and pie1 single mutants (Figure S14). Once again, however, the magnitude of this change is extremely small.

The defining characteristic of ibm1 mutants is a major increase in genic CHG methylation [63], [64], [68]. The h2a.z;ibm1 and pie1;ibm1 double mutant lines were generated such that h2a.z;ibm1 seedlings were newly homozygous for ibm1 (first generation), whereas pie1;ibm1 seedlings came from first generation ibm1 homozygous parents (second generation). The h2a.z;ibm1 seedlings in their first generation of ibm1 homozygosity have higher levels of CHG methylation than first generation ibm1 seedlings, and pie1;ibm1 seedlings in their second generation of ibm1 homozygosity have lower levels of CHG methylation than second generation ibm1 seedlings (Figure S12). Both first generation datasets, h2a.z;ibm1 and ibm1, show similar levels of CHH hypermethylation to one another; likewise, the second generation pie1;ibm1 and ibm1 data exhibit similar CHH hypermethylation levels (Figure S13). Importantly, the control data show that genic CHG and CHH methylation is unstable in ibm1, increasing greatly in the second generation (Figure S12), making interpretation of changes in h2a.z;ibm1 and pie1;ibm1 CHG methylation difficult.

Whereas there is little difference between the double mutant lines and their parental lines in TE CG methylation (Figure S15), we found CHG hypomethylation in the double mutants as compared to their respective parental lines (Figure 3F and Figure S16). Additionally, while CHH methylation is unaltered in h2a.z and pie1 seedlings, there is a significant reduction of TE CHH methylation in h2a.z;ibm1, h2a.z;met1, and pie1;ibm1 seedlings compared to the ibm1 and met1 single mutants, which is similar to the reduction seen in h2a.z;pie1 seedlings (Figure 3C and Figure S17). Taken together, our results suggest that while H2A.Z may play a modest role in the regulation of DNA methylation in TEs, the genome-wide anticorrelation between H2A.Z and DNA methylation is due to DNA methylation preventing the incorporation of H2A.Z.

H2A.Z is enriched in responsive genes

Given the published work linking H2A.Z with regulation of several types of genes that respond to the environment [8], [21], [26], [27], [28], [37], [69], we sought to examine H2A.Z enrichment with respect to gene responsiveness. To do so, we generated a genome-wide map of H2A.Z using our published tagged H2A.Z Arabidopsis line [16] by coupling affinity purification of H2A.Z-bound DNA with high-throughput sequencing. Metaanalyses of the new dataset demonstrate a strong peak of H2A.Z at the 5′ end and a smaller peak at the 3′ end of most genes, with varying levels of H2A.Z distributed within gene bodies (Figure 4A–4C). Our new data are consistent with our published microarray results, though the resolution provided by high-throughput sequencing is significantly better (Figure S18). To determine if a loss of H2A.Z had a preferential effect on genic methylation in the subset of genes that were enriched for H2A.Z in WT, as suggested by our published results in pie1, we compared the h2a.z and WT bisulfite sequencing datasets for those genes with the most and least H2A.Z across the TSS and the gene-body (Figure S19). Although there are major methylation differences between these four groups of genes, the profiles of h2a.z and WT within each group were virtually indistinguishable from one another. However, the subtle increase in genic methylation we detected in h2a.z, pie1, and h2a.z;pie1 plants was stronger in H2A.Z enriched genes (Figure S7), consistent with our published results [16].

In support of earlier data [16], [24], we found a negative correlation between H2A.Z enrichment in gene bodies and WT transcript levels (Spearman's rho = −0.4039, P-value<0.0001). Genes with the most gene body H2A.Z (n = 4,081, Table S1) have median WT expression more than six-fold lower than that of genes with the lowest H2A.Z within their bodies (n = 3,920, Table S2) (Figure 4D). In comparison, levels of H2A.Z enrichment near the TSS showed a different trend: genes with the most and least H2A.Z at the TSS had lower levels of expression than those with intermediate levels of H2A.Z (Figure 4D), as we showed earlier for both Arabidopsis and pufferfish [16], [24].

We also discovered a positive correlation between enrichment of H2A.Z across gene bodies and gene responsiveness – the degree to which a gene is differentially expressed among different tissue types or experimental conditions (including hormone, nutrient, and chemical treatments, as well as biotic or abiotic stimulus), with higher response scores associated with greater differential expression [70]. H2A.Z body-enriched genes (n = 4,081) have a six-fold higher median gene responsiveness score than that of genes with the lowest H2A.Z levels across their bodies (n = 3,920) (Figure 4E). Levels of H2A.Z at the TSS are considerably less correlated with response score than levels of H2A.Z in the body (Spearman's rho = 0.0748 and 0.3325, P-values<0.0001, respectively), and highly responsive genes have more body H2A.Z than genes with low responsiveness (Figure S20). The least and most responsive genes [70], defined as housekeeping genes (n = 371) and hypervariable genes (n = 117) [70], are depleted for and enriched in H2A.Z across the gene body, respectively (Figure 4A and Figure S20). These results suggest that H2A.Z deposition in the gene body may facilitate rapid activation or inactivation of genes.

H2A.Z regulates responsive genes

To determine which genes are misregulated upon loss of H2A.Z, we profiled the transcriptomes of the h2a.z mutant and WT in 4-week old rosette leaves with three replicates each of RNA sequencing. 1,800 genes were upregulated and 544 genes were downregulated in h2a.z with a P-value cut-off of 0.001. This is consistent with transcriptome analyses of hta9;hta11 and pie1, which showed three-fold and two-fold more genes upregulated than downregulated, respectively [37]. The genes exhibiting up and downregulation in the h2a.z mutant show statistically significant overlap with lists of up and downregulated genes in pie1 and hta9;hta11 [37], despite differences in tissues type, developmental stage, growth conditions, and transcriptional profiling platform used to generate these data (Figure S21). Gene Ontology analysis of the misregulated genes in h2a.z revealed enrichment of categories related to immune response (P-value = 8.6×10⁻⁹) and temperature response (P-value = 4.8×10⁻⁸), consistent with previous studies of pie1 and arp4 mutants [8], [37] (Table S3). Strikingly, all of the most-enriched categories (P-value<1×10⁻⁵) are specifically response-related, and include many previously unreported GO-terms involved in the perception of a variety of external cues (Figure 5A). Many of these GO terms are also overrepresented in the smaller subset of genes upregulated in at least two of the h2a.z, pie1, and hta9;hta11 mutant datasets (Figure 5A).

Consistent with our Gene Ontology analysis (Tables S3, S4), we discovered a relationship between the degree of misregulation in the h2a.z mutant and the responsiveness score of a gene (Figure 5B). Genes exhibiting greater than 4-fold upregulation (n = 938) had a 2.5-fold higher median responsiveness score than that of the least upregulated genes (less than 1.4-fold up or downregulated, n = 9,300). The relationship between downregulation and response score, on the other hand, was roughly parabolic, with the most downregulated and least downregulated genes showing the lowest levels of responsiveness, and genes with intermediate levels of downregulation (2 to 4-fold) showing the greatest responsiveness (Figure 5B). Hypervariable genes are generally strongly upregulated in h2a.z plants, despite a lack of change in DNA methylation (Figure S20), whereas the expression of housekeeping genes is largely unchanged (Figure 5C).

Because H2A.Z is enriched in bodies of response genes, we investigated whether changes in transcriptional regulation in the h2a.z mutant correlated with specific H2A.Z enrichment patterns in WT. As expected, we found a positive relationship between misregulation in the h2a.z line and H2A.Z gene body enrichment (Figure 5D) (Spearman's rho = 0.2634 for downregulated genes and 0.2540 for upregulated genes, P-value<0.0001). Genes with the greatest misregulation (greater than four-fold up or downregulated, n = 1,258) have more than a 36-fold higher median H2A.Z-body enrichment score than that of genes with the lowest levels of change in transcription between h2a.z and WT (less than 1.4-fold up or downregulated, n = 9,300). Taken together, our data demonstrate that loss of H2A.Z leads to a general transcriptional misregulation of response genes that are enriched for H2A.Z within the gene body in wild type, including genes that respond to developmental, biotic, and abiotic stimuli (Figure 5E). Our results also suggest that one function of gene body methylation, which is strongly anticorrelated with gene responsiveness in plants and animals [70], [71], [72], is the exclusion of H2A.Z from the bodies of constitutively expressed genes.

Discussion

We have generated a viable H2A.Z-deficient mutant line in Arabidopsis thaliana that shares many, but not all of the phenotypic characteristics of pie1 mutants. We show that unlike in yeast, combining Arabidopsis h2a.z and pie1 mutations exacerbates the phenotype. Loss of H2A.Z does not substantially affect the level or profile of DNA methylation in genes, even when combined with mutations that alter the normal genic methylation landscape, whereas DNA methylation in transposons is perturbed in a tissue-dependent manner. We show that differences in gene body H2A.Z levels are correlated with gene expression and gene responsiveness. Finally, we show that a loss of H2A.Z causes misregulation of many genes involved in the response to developmental and environmental cues, and that these genes tend to have high levels of gene-body H2A.Z.

Residual H2A.Z function remains in pie1 mutant plants

Whereas the fungi S. pombe and S. cerevisiae can tolerate mutations in H2A.Z [27], [73], H2A.Z is essential in many species, including Tetrahymena thermophila, Drosophila melanogastor, Xenopus laevis, Caenorhabditis elegans and mice [15], [56], [58], [74], [75]. Consequently, many studies of H2A.Z function outside of yeast have utilized mutants in components of the chromatin remodelers that deposit H2A.Z to emulate H2A.Z loss-of-function [8], [34], [37], [69]. The substantial overlap between the phenotypes of Arabidopsis pie1 and h2a.z mutants suggests that PIE1 is the primary remodeler responsible for H2A.Z deposition. However, h2a.z;pie1 double mutants exhibit early developmental arrest not seen in either of the single mutant lines, indicating that H2A.Z may be deposited in the absence of the PIE1 complex, potentially by the Arabidopsis homolog of INO80 [76], which can deposit H2A.Z in yeast [51]. The PIE1 complex might also have H2A.Z-indpendent roles, as has been hypothesized for the PIE1/SWR1 orthologs in animals [49], [50]. Indeed, a recent study showed that H2A.Z deposition by p400 and SRCAP, the human orthologs of SWR1, could not account for all the regulatory roles of these complexes [50]. These results emphasize that phenotypes caused by mutations in chromatin remodeling complexes must be interpreted with caution.

DNA methylation excludes H2A.Z from chromatin

DNA methylation and H2A.Z are tightly anticorrelated in plants and animals [16], [24], [47], [48], and we have shown that DNA methylation quantitatively excludes H2A.Z from chromatin [16]. Here, we demonstrate that H2A.Z does not have a large influence on DNA methylation in genes, even when genic DNA methylation is in flux, but that loss of H2A.Z does cause a small increase in genic methylation, particularly in H2A.Z-enriched genes, consistent with our earlier results [16]. The magnitude of these changes is unlikely to substantially contribute to the genome-wide anticorrelation between DNA methylation and H2A.Z, indicating that exclusion of H2A.Z from methylated DNA is the cause (Figure 6). The bisulfite sequencing data also reveal global decreases in CHG and CHH TE methylation in the h2a.z mutant. Changes in TE methylation could be a direct result of H2A.Z loss, or may be caused by a variety of indirect effects. Given the depletion of H2A.Z from methylated transposons and the substantial transcriptional and developmental changes in h2a.z plants, we consider indirect explanations to be more probable. For example, the approximately two-fold downregulation of the DNA methyltransferase CMT3, which catalyzes CHG methylation [77], might be partly responsible for the decreased CHG methylation. (Table S5).

H2A.Z in gene bodies regulates transcription of responsive genes

The significance of H2A.Z enrichment near transcriptional start sites has been a major focus of research [22], [23], [78], [79], [80], but a distinct function for H2A.Z in gene bodies has been recently hypothesized [81]. Consistent with this idea, we previously showed that H2A.Z abundance within gene bodies correlates negatively with transcription in Arabidopsis and the pufferfish Tetraodon nigroviridis, whereas H2A.Z near the TSS is most enriched in moderately transcribed genes in both organisms [16], [24]. Human studies also show that gene body H2A.Z correlates with silencing [12] and that H2A.Z is depleted from the bodies of actively transcribed genes [23]. The presence of this relationship in plants and animals implies that it is an ancient property of eukaryotes. Interestingly, recent studies in yeast have shown that mutation of the IN080 complex causes loss of H2A.Z near the TSS and gain of H2A.Z across the coding region [51], suggesting that competing nucleosome remodelers may shape the genic patterns of H2A.Z.

Here, we show that H2A.Z within gene bodies is correlated with gene responsiveness, consistent with recent yeast data demonstrating that H2A.Z is enriched across coding sequences of genes that are differentially transcribed after environmental stress [26]. Loss of H2A.Z leads to misregulation of Arabidopsis genes with high responsiveness scores, which measure differential expression across both tissue types and environmental conditions. Furthermore, this misregulation occurs despite a lack of change in the DNA methylation profiles of these genes in the h2a.z mutant. Our results are consistent with evidence from many other species, where loss of H2A.Z leads to misregulation of various inducible genes, including environmental response genes in yeast [21], [26], [78] and developmentally regulated and tissue-specific genes in animals [13], [15], [18], [82], [83]. The phenotypes of h2a.z mutants, including altered flowering time, floral homeotic transformations and silique deformation, also strongly imply that developmental regulators are misregulated. We also demonstrate that genes that show little change in transcription in our h2a.z mutant plants tend to have H2A.Z depleted from the gene body, whereas those genes with either strong up- or downregulation tend to have much more gene-body H2A.Z. Taken together, these results indicate that H2A.Z within transcribed sequences is necessary for proper regulation of responsive genes but may antagonize constitutive and high-level expression, and that this relationship is both ancient and well-conserved across many eukaryotic lineages.

Gene body methylation may regulate gene expression by preventing H2A.Z incorporation

The presence of DNA methylation within the bodies of animal and plant genes has been known for some time [84], [85]. Recent genome-wide bisulfite sequencing in various eukaryotic species has revealed that gene body methylation is an ancient and widely conserved feature of eukaryotic chromatin predating the divergence of animals and plants [24], [38], [43], [44], [45], [71], [86], [87]. In animals and plants, gene body methylation exists almost exclusively within the CG context and follows a consistent pattern, with depletion of DNA methylation from the 5′ and 3′ ends of genes. Taken together with the finding that many species of invertebrates have DNA methylation primarily or exclusively within gene bodies [24], [87], [88], these results strongly suggest that genic methylation plays an important and conserved function in at least some eukaryotic lineages [89].

Despite the prevalence of gene body methylation in diverse eukaryotes, its function remains mysterious [90]. A potential clue comes from the correlation between genic methylation and transcription. Gene body methylation is highest in moderately transcribed genes in plants and animals, with the lowest levels of genic methylation at either transcriptional extreme [24], [45], [71]. Additionally, there is an unexplained negative linear correlation between genic methylation and gene responsiveness in Arabidopsis and the honeybee Apis mellifera [70], [71], [91]. High levels of body methylation tend to be found in slowly evolving genes with vital housekeeping functions in honeybee, silkworm (Bombyx mori), sea squirt (Ciona intestinalis), sea anemone (Nematostella vectensis), poplar (Populus tricharpa), and Arabidopsis [71], [88], [92], [93]. These results indicate that DNA methylation of the transcribed region may be important for proper regulation of constitutively expressed genes.

Here, we show that the genome-wide anticorrelation between DNA methylation and H2A.Z is established by the exclusion of H2A.Z from methylated DNA. Because gene body DNA methylation and H2A.Z show opposing correlations with gene responsiveness, and the anticorrelation between DNA methylation and H2A.Z is ancient, we propose that a basal function of genic DNA methylation is the stabilization of constitutive expression patterns within housekeeping genes by antagonizing H2A.Z deposition (Figure 6). As H2A.Z has been linked to the regulation of inducible genes in many organisms, including species such as S. cerevisae and C. elegans that lack DNA methylation [8], [13], [15], [18], [21], [26], [69], [78], [82], [94], and DNA methylation can exclude H2A.Z but not vice versa, we believe that the presence or absence of H2A.Z in the gene body is a better candidate for direct gene regulation than DNA methylation. The functional significance of DNA methylation of constitutive genes may be primarily to prevent incorporation of H2A.Z.

Materials and Methods

Biological materials

The Arabidopsis T-DNA lines hta9-1 (SALK_054814), hta11-1 (SALK_017235), ibm1–6 (SALK_006042), and pie1–5 (SALK_096434) were obtained from the SALK collection (Col-0 ecotype) (http://signal.salk.edu/). The Arabidopsis T-DNA line hta8-1 (FLAG_593B04) was obtained from the INRA (http://www-ijpb.versailles.inra.fr/) collection (WS ecotype). Sequencing of the 5′ promoter region of HTA8 confirmed the T-DNA insertion site for hta8-1 at position 16,220,917 on Chr2 (NC_003071.1), 8 bp downstream of the 5′ end of gene model AT2G38810.2. The Arabidopsis EMS mutant met1–6 is described in [95] (Col-0 ecotype). The h2a.z mutant was generated from crosses of an hta9-1; hta-11 double mutant with the hta8-1 line. In all experiments, the WT control line for each mutant was generated from the nearest WT sibling of the mutant (i.e. the WT control for the h2a.z mutant is an F2 progeny of wild-type genotype at all three H2A.Z loci, derived from an F1 parent that is heterozygous for hta9-1, hta11-1, and hta8-1). In instances where we felt it was necessary to distinguish which WT control we refer to (i.e. Figure S6), we use the non-mutant gene identifier (in all capital letters) preceding the “WT” (i.e., the WT line associated with the pie1 mutant is “PIE1 WT”). For bisulfite sequencing of seedling tissues, seeds were planted on 1× Murashige and Skoog Media with micronutrients and 1.5% Sucrose (Caisson Laboratories) and grown under 16 h light/8 h dark for 14 days in a growth chamber. For bisulfite sequencing of rosette and cauline leaf tissue, seeds were planted on soil and grown in greenhouse conditions with LD 16 h light/8 h dark. For phenotype analysis of the h2a.z mutant, seeds were planted on soil and grown in greenhouse conditions with either 16 h light/8 h dark (LD) or 8 h light/16 h dark (SD). Genotyping of SALK and INRA T-DNA lines was carried out by PCR with primers listed in Table S6. Genotyping of the met1–6 line was carried out by dCAPS-PCR with primers listed in Table S6 and subsequent digestion with BglII.

Transcript analysis of H2A.Z genes

Expression analyses for the h2a.z and hta9-1 mutant lines and for the WT control were performed on total RNA extracted from 4 week post germination rosette leaves grown on soil in LD conditions using the RNeasy Plant Extraction Kit (Qiagen) with the optional on-column DNAse treatment. RT-PCR reactions were carried out on cDNA generated using 1 ug total RNA and the Superscript III Kit (Invitrogen) using gene specific primers listed in Table S6. qPCR was carried out on similarly generated cDNA using EvaGreen Detection chemistry on an ABI 7500 FAST Real-Time PCR System with primers in exons flanking the single intron in HTA9. The gene UBQ5 (AT3G62250) was used as in internal control. Three biological replicates, each with three technical replicates, were averaged.

Bisulfite sequencing

Approximately 100–500 ng genomic DNA was isolated from either seedling, rosette or cauline leaf tissues. Seedling tissue was obtained from 14 days post germination seedlings grown on Murashige and Skoog media in LD (16 h light/8 h dark). Mature rosette leaves and mature cauline leaves were obtained from 4 week post germination mature plants grown on soil in LD (16 h light/8 h dark). In general, multiple biological replicates were generated for each mutant and WT line; a complete list of all generated libraries is available in Table S7. WT datasets for each mutant were generated from plants derived from recent relatives of the relevant mutant. Bisulfite conversion and Illumina library construction and sequencing were performed as described in [96]. We used single ends (SE) Illumina sequencing for bisulfite sequencing on the GAII and HiSeq platforms and sequence alignments were performed using Bowtie [97] and the TAIR8 Genome Annotation (http://www.arabidopsis.org/) as in [24]. The average percent methylation plots were generated as described in [96] and [24].

Locus-specific methylation analysis

For locus-specific bisulfite sequencing (referred to as BS-PCR), data were generated exactly as described previously [16]. Bisulfite-converted genomic DNA from 14 day seedlings was PCR amplified using primers from each of four genotypes, including h2a.z, pie1, H2A.Z WT, and PIE1 WT. Approximately 10–12 clones were sequenced from each genotype (except for At3g22340, in which 6 clones were sequenced for pie1) and percent methylation was determined at the same cytosine sites used for this calculation in our previous publication. Alternatively, percent methylation scores were also calculated by extracting the reads associated with each locus from the relevant whole genome bisulfite sequencing datasets (referred to as BS-Seq). Average percent methylation levels were calculated at the same cytosine sites described above for BS-PCR from available seedling replicates for each genotype (2 for pie1, 3 for PIE1 WT, 2 for h2a.z, 2 for H2A.Z WT).

RNA sequencing

Approximately 30 ug total RNA was isolated from 4 week post germination mature rosette leaves using the RNeasy Plant Extraction Kit (Qiagen) with the optional on-column DNAse treatment. mRNA was purified from total RNA by two treatments of poly-A enrichment using the Oligotex kit (Qiagen #72022), followed by a rRNA removal step using the RiboMinus Plant Kit for RNA sequencing (Invitrogen #A1083702). Illumina library construction and RNA sequencing were performed as described in [24]. We used single ends (SE) Illumina sequencing for RNA sequencing on the GAII platform and sequence alignments were performed using Bowtie and the TAIR8 Genome Annotation and cDNA Annotation (http://www.arabidopsis.org/) as in [24].

H2A.Z ChAP sequencing

H2A.Z-containing nucleosomes were chromatin affinity purified (ChAP) from 4 week post germination Arabidopsis roots of our H2A.Z-BLRP transgenic lines grown in LD conditions as in [16]. Illumina libraries were constructed for IP and input DNA samples and sequenced on the HiSeq 2000 platform, generating 50 bp reads. Sequence alignments were performed using Bowtie and the TAIR8 genome annotation as in [24]. Nucleosomal midpoints were estimated based on an average 150-bp nucleosome length by adding 75 bp to the start position of each read. Differences between IP and input over each single-base window were generated to give an overall genome-wide map of H2A.Z-enrichment.

Differential expression sequence analysis

For differential expression analysis of the RNA sequencing datasets, a strategy was employed to account for expression differences between WS and Col ecotypes. In brief, we used the recently published list of 144,879 SNPs between the WS and Col ecotypes [98] to obtain reads per kilobase of exon model per million reads (RPKM) scores for each gene in h2a.z and WT from either the WS or Col backgrounds.

First, using Bowtie with no tolerance for mismatches, reads from each of the three h2a.z and WT RNA sequencing datasets were mapped to small 75 bp scaffolds containing either the WS or Col SNP around each SNP locus that mapped within an exon of a gene greater than 200 bp in length and with at least 10 mapped reads. We removed all SNPs that were less than one read-length (36 bp) from the end of the exon, which left approximately 5,000 SNPs across the genome. The number of reads mapping to the WS and Col scaffolds were compared at each SNP locus and used to determine whether the region was homozygous for WS, Col or heterozygous for the two ecotypes in each dataset. For SNPs at heterozygous loci, a Read Count Contribution from each WS or Col genome was determined by dividing the number of reads mapping to either WS or Col genome by the total reads mapping to the SNP scaffold for each ecotype. As SNPs within a given heterozygous region generally exhibited similar ratios of WS to Col mapped reads, a rolling 20-window (where the windows are the 5,000 SNPs) smoothing function was applied to these read count contribution values.

Next, the six RNA sequencing datasets were mapped to the TAIR cDNA scaffold, and each cDNA model was assigned a score equal to the number of mapped RPKM. For both the h2a.z and WT datasets, the normalized read counts of the three replicates were partitioned into reads contributed by WS and by Col using the smoothed read count contribution value obtained from the nearest SNP. In this way, approximate WS and Col read count scores were determined for each gene in both h2a.z and WT.

To test for statistical significance of the difference between the h2a.z and WT, we repeated the above partitioning process using read counts normalized to the size of the smallest library, rather than per million of reads. This alternate normalization less drastically underestimates the number of reads per locus, which is important as the statistical significance is dependent on the number of reads. We calculated the probability that a gene's expression deviates from expectation using a Fisher's two-tailed exact test of h2a.z vs. WT scores for each ecotype. Genes were determined to be differentially expressed if for either ecotype they exhibited a two-fold change in expression between h2a.z and WT and had a P-value<0.001, or if for both ecotypes they exhibited a two-fold change in expression and had p-values<0.005. Gene Ontology analysis was performed on the up- and downregulated gene lists using the GO FAT Ontology on the DAVID web server (http://david.abcc.ncifcrf.gov) [99], [100] and categories with P-values<1×10⁻⁵ were considered enriched.

Data deposition

Sequences are deposited in Gene Expression Omnibus (GEO) with accession number GSE39045.

Supporting Information

Attachment 1

Attachment 2

Attachment 3

Attachment 4

Attachment 5

Attachment 6

Attachment 7

Attachment 8

Attachment 9

Attachment 10

Attachment 11

Attachment 12

Attachment 13

Attachment 14

Attachment 15

Attachment 16

Attachment 17

Attachment 18

Attachment 19

Attachment 20

Attachment 21

Attachment 22

Attachment 23

Attachment 24

Attachment 25

Attachment 26

Attachment 27

Attachment 28