Herbicide injury induces DNA methylome alterations in Arabidopsis

Introduction

The development of herbicide-resistant weed populations is a major crisis facing modern agriculture (Bonny, 2016). Herbicide resistance has evolved in populations of at least 470 different weeds, including 35 species that have developed resistance to the herbicide glyphosate (Heap, 2016), the most widely used herbicide in the history of agriculture (Duke & Powles, 2008). Mechanisms of herbicide resistance are easily explained in some cases, as when a mutation occurs in the protein target of the herbicide that reduces herbicide binding (González-Torralva et al., 2012) or even the copy number of herbicide target genes (Gaines et al., 2010). Other cases, termed non-target site resistance, are poorly understood, but have been attributed to quantitative accumulation of minor resistance alleles selected from standing genetic variation at multiple gene loci (Busi, Neve & Powles, 2013; Délye, 2013). Non-target site resistance may involve various mechanisms that affect herbicide metabolism or translocation (González-Torralva et al., 2012). Resistance to glyphosate is interesting in that it often appears to involve non-target site mechanisms, which can emerge and dominate a population after as few as three generations of sub-lethal exposure to the herbicide (Busi & Powles, 2009), and can spread through a population faster than predicted by gene flow or propagule dispersal (Escorial et al., 2011; Espeby, Fogelfors & Milberg, 2011; Okada et al., 2013). Given that herbicides induce a strong abiotic stress, it is likely that weeds respond by activating stress-signaling networks that reprogram gene expression (Busi, Neve & Powles, 2013), and we hypothesize that this involves epigenetic regulation of gene function that could contribute to herbicide resistance as implicated in weed adaptation to other stresses (Verhoeven et al., 2010) and in the response of rice to the pesticide atrazine (Lu et al., 2016).

A prominent mechanism of epigenetic regulation of gene expression is through cytosine methylation. Methylated cytosines (mCs) occur in three different sequence contexts in plants: CG methylation, which is the predominant form of methylation in gene-rich regions, and CHH or CHG (H = A, T, or C) methylation predominantly in transposable elements and repetitive sequences (Cokus et al., 2008). CHH and CHG mCs are much more common in plants than animals (Cokus et al., 2008). Multiple RNA-dependent methylation pathways control de novo methylation in the mustard weed Arabidopsis thaliana (Matzke & Mosher, 2014) with methylation at the different sequence contexts altered and maintained through both overlapping and sequence-context-specific mechanisms. Mutations in the genes of these pathways lead to aberrations in the methylome (Stroud et al., 2013). Additionally, wildtype A. thaliana plants accumulate epimutations (i.e., changes in the methylome) over generations of greenhouse propagation (Becker et al., 2011; Schmitz et al., 2011). Stress exposure can lead to substantial methylome reprogramming as observed in A. thaliana following pathogen attack or salicylic acid treatment (Dowen et al., 2012), phosphate starvation (Secco et al., 2015; Yong-Villalobos et al., 2015), and salinity stress (Jiang et al., 2014; Wibowo et al., 2016). The majority of methylation changes following these stresses were considered transient—not transgenerationally stable—but a subset of stress-responsive methylations is stably fixed in populations over several subsequent generations of exposure to stress (Jiang et al., 2014). Thus, transgenerational epigenetic changes in gene expression could lead to enhanced and adaptive stress tolerance. However, the ability of epigenetic changes to be transgenerationally stable and affect the adaptation of plants in their environment remains contentious (Hagmann et al., 2015; Kawakatsu et al., 2016).

Because the specificity of methylome reprogramming following stress may appear inexact [e.g., compare Secco et al. (2015) to Yong-Villalobos et al. (2015)] and the effect of alterations in methylation levels on gene expression may seem inconsistent (Karan et al., 2012; Kawakatsu et al., 2016; Li et al., 2012), we tested whether there was a dosage effect by glyphosate stress on specific alterations of the plant methylome. Glyphosate has been previously shown to alter DNA methylation in wheat (Nardemir et al., 2015) using relatively low resolution techniques. In addition, the herbicide atrazine was recently shown to induce global changes in methylation in rice (Lu et al., 2016). We used bisulfite sequencing to determine the effect of glyphosate on the Arabidopsis methylome at single base pair resolution to identify specific genetic loci that have altered methylation patterns following glyphosate exposure.

Materials and Methods

Results and Discussion

To determine whether glyphosate injury induces changes in plant methylation, we treated A. thaliana with low concentration glyphosate sprays representative of real-world doses a weed on the margin of a treated field could receive. Because a plant must survive herbicide injury and produce viable seed to be a founder of an herbicide resistant population, we first identified the appropriate sub-lethal doses for glyphosate treatment of A. thaliana. Four-week-old A. thaliana rosettes were exposed to 0%, 5%, 10%, or 15% of a typical field rate of 0.9 kg acid equivalency (ae) ha⁻¹ glyphosate, with the 5% and 10% rates causing visible herbicide injury, but allowing for plant survival and reproduction (Fig. 1A). gDNA was collected from newly formed cauline leaves at silique maturation of four individuals from each treatment (Figs. 1B–1D). MethylC-seq of 12 total libraries (four replicates from each treatment) produced a total of 530,042,668 aligned reads resulting in genome coverages from 48× to 76× for each replicate (Table S1).

The total amount of methylated cytosines (mCs) ranged from 95 to 143 million for all replicates (Table S1) and were the most abundant in the CG sequence context (Fig. 2A), which is consistent with previous analyses of the A. thaliana methylome (Cokus et al., 2008). Only an average of 0.2% mCs was identified on the chloroplast genome, which is not methylated in A. thaliana (Cokus et al., 2008), demonstrating the fidelity of the methylC-seq protocol and data filtering. mCs were distributed across all five nuclear chromosomes. In addition, a small number of mCs were identified in the mitochondrial genome though these mCs are likely false positives because of the insertion of a region of the mitochondrial genome on chromosome 2 of A. thaliana eco. Columbia (Lin et al., 1999). Neither the abundance of mCs nor the relative frequency of mCs in the three different analyzed sequence contexts (CG, CHG, and CHH) significantly changed due to herbicide injury (Fig. 2A). Clustering the Pearson’s correlation coefficient of the mCs did not result in treatment-specific monophyletic clades, demonstrating the relative similarity of the methylomes across all 12 individuals (Fig. S3). Taken together, these results suggest that glyphosate-induced stress did not lead to changes in overall levels of methylation in any of the three sequence contexts.

Because glyphosate did not induce global shifts in methylation levels, differences in specific location and context of DNA methylation can be potentially classified as selective responses to herbicide injury. We used Methylkit (Akalin et al., 2012) to compare mCs in the 5% or 10% glyphosate treatments to the control (0% treatment) to identify differentially methylated cytosines (DmCs). This approach enumerated 17,017 DmCs following 5% glyphosate treatment and 23,341 DmCs following 10% glyphosate treatment after normalizing all four replicates and applying a cutoff to only consider the top 25% highest confidence DmCs (Dataset S1). Interestingly, 10% glyphosate treatment led to substantially more DmCs than 5% glyphosate treatment in the two most abundant sequence contexts CG and CHG (Figs. 2B–2D), suggesting that the degree of herbicide injury is correlated with the context and magnitude of methylation changes. However, it is important to note that the 10% glyphosate-treated plants required an additional six weeks to reach silique maturity. Therefore, a subset of the methylation differences could be due to the increased growing time and not the herbicide treatment.

Differentially methylated regions are more strongly associated with regulatory changes in gene expression than DmCs. We called DMRs induced by both 5% and 10% glyphosate treatment using two independent approaches: a bimodal DmC distribution modeling approach using eDMR (Li et al., 2013), and a curve smoothing approach of DmCs using bsseq (Hansen, Langmead & Irizarry, 2012). All DMRs were defined as CG, CHH, or CHG based on the predominance of the sequence context of DmCs comprising the DMR. The eDMR algorithm identified 1,949 hypomethylated and 1,229 hypermethylated non-redundant DMRs following 5% and 10% glyphosate treatments (Dataset S2). More than 95% of the DMRs called by eDMR were defined by the CG sequence context (Fig. S4). Using the bsseq algorithm, 4,053 hypomethylated DMRs and 5,082 hypermethylated DMRs were identified following glyphosate herbicide injury across all three sequence contexts (Dataset S3). We excluded 70 DMRs identified in both the 5% and 10% treatment groups but conflicting in directional change of methylation (i.e., hypo in one treatment and hyper in the other treatment) and therefore unlikely to represent biologically relevant DMRs. In contrast to the eDMR data, only 78% and 51% of hypomethylated and hypermethylated bsseq DMRs, respectively, were categorized as CG (Fig. 3A). The lack of congruence between the bsseq and eDMR datasets suggests that the choice of algorithm remains a primary driver for the specificity of DMRs identified as previously suggested (Yu & Sun, 2016). We consider the 545 DMRs identified by both eDMR and bsseq as our highest confidence DMRs (Dataset S3). We focused on the bsseq DMRs for subsequent analyses because of the relative abundance of non-CG DMRs in this dataset and the observation that non-CG DMRs are more frequent in plants than they are in animals.

We identified 1,964 DMRs (counting both hypo- and hyper-methylated DMRs) unique to the 5% glyphosate-treated samples, and 5,614 DMRs unique to the 10% samples (Fig. 3A), which are consistent with the trend observed for DmCs in which the severity of herbicide injury correlates with the magnitude of methylome shifts. The 5,614 DMRs present only in the 10% glyphosate-treated samples were tentatively classified as positive dose-dependent DMRs (i.e., the magnitude of changes in methylation levels is positively correlated with severity of herbicide injury) and the 1,964 DMRs unique to the 5% glyphosate-treated samples as inverse dose-dependent DMRs (i.e., more severe herbicide injury correlated with decreases in the magnitude of changes in methylation levels). The 1,557 DMRs present in both the 5% and 10% treated samples were further analyzed for dose-dependency resulting in 224, 300, and 209 DMRs classified with 95% confidence as dose-independent, positive dose-dependent, and inverse dose-dependent respectively, across both hypomethylated and hypermethylated DMRs (Fig. 3B; Dataset S4; Fig. S2; see Methods). We, therefore, conclude that for nearly two-thirds (5,914 out of 9,205) of DMRs that are positive dose-dependent, the magnitude of changes in methylation is dependent on the severity of the stress. We propose that these positive dose-dependent DMRs are the best candidates for identifying biologically relevant genomic loci that are modified in response to glyphosate stress. However, it is important to note that the 10% glyphosate-treated plants were harvested several weeks after the 5% glyphosate-treated plants to compensate for the developmental delay induced in the 10% glyphosate-treated plants (see Methods). It is possible that a number of the identified differential DMRs between the two treatment levels are due to maintenance of these DMRs being uncoupled from development and instead linked to the age of the plants.

We tentatively hypothesize for the smaller number of the inverse dose-dependent DMRs that the severity of herbicide injury following 10% glyphosate treatment prevented the plant from responding to the stress through specific alterations in DNA methylation. Interestingly, the CHH context had the most pronounced inverse dose-dependent trend, with more DmCs in the 5% treatment group than the 10% treatment group (Fig. 2), in stark contrast to the CG and CHG contexts. While some molecular mechanisms that control methylation act on all three sequence contexts [e.g., DRM 1/2 (Cao & Jacobsen, 2002b)], other mechanisms have distinct effects across different sequence contexts with the most pronounced differences between asymmetrical and symmetrical Cs (Cao & Jacobsen, 2002a; Matzke & Mosher, 2014). CHH sites are asymmetrical (no C in the antisense strand) in contrast to CG and CHG sites. Additionally, CHH sites were previously identified as responding uniquely to biotic stress in contrast to CG and CHG sites (Dowen et al., 2012). However, it is impossible to conclude that glyphosate herbicide injury alters different molecular methylation mechanisms specifically until the effects of glyphosate on methylation patterns in known methylation mutant plant lines are empirically assessed.

To elucidate which genomic regions are responding to the glyphosate stress through changes in methylation, we annotated all DMRs. A plurality of DMRs were identified in gene coding sequences and in the CG context (Fig. 4A) as expected (Cokus et al., 2008). The DMRs categorized as CHG or CHH were predominantly associated with transposable elements, similar to previous studies in plants (Cokus et al., 2008; Li et al., 2012). In total, 1,818 transposable elements or transposable element (TE) genes were hypermethylated and 936 were hypomethylated (Fig. 4A; Dataset S3). Transposable element-associated DMRs occurred at varying frequency across superfamily groups in response to both the 5% and 10% glyphosate treatments (Fig. S5; Dataset S5). Even though hypermethylated TEs are twice as common as hypomethylated TEs, the 936 hypomethylated TEs may contribute to genome destabilization due to increased transposon mobility associated with reduced transposon body methylation (Chan, Henderson & Jacobsen, 2005; Mirouze & Vitte, 2014). The number of hypomethylated TEs is less than the number of activated TEs in the DNA methylation A. thaliana mutant ddm1 (Zemach et al., 2013), but still greater than the number of TEs activated in A. thaliana following treatment with DNA demethylating agent azacytidine based on expression data (Griffin, Niederhuth & Schmitz, 2016). We, therefore, conclude that glyphosate induces hypomethylation of TEs on a scale comparable to other known treatments capable of destabilizing the A. thaliana genome through TE hypomethylation. The hypomethylated TEs are interesting to explore in terms of their potential for contributing to the phenomenon of gene amplification observed in several glyphosate-resistant weeds (Gaines et al., 2010).

An advantage to studying plant DNA methylation changes in response to glyphosate is that the biochemical mechanism of action of this herbicide has been extensively characterized. Glyphosate stops the flow of carbon through the shikimate pathway by inhibiting activity of the enzyme 5-enolpyruvylshikimate 3-phosphate (EPSP) synthase (Amrhein et al., 1980). Two of the seven genes in the shikimate pathway of Arabidopsis (Tzin & Galili, 2010) were differentially methylated in the CG context in response to glyphosate stress: 3-deoxy-d-arabino-heptulosonate-7-phosphate (DAHP) synthase (At4g33510), which encodes the first enzyme (and major regulatory step for the pathway), and both forms of shikimate kinase (SK1, At2g21940 and SK2, At4g39540) were all hypermethylated in the 10% glyphosate-treated plants. Notably, all three of these genes were previously classified as unmethylated in A. thaliana (Niederhuth et al., 2016). EPSP synthase (At2g45300), which is also considered unmethylated in A. thaliana populations (Niederhuth et al., 2016), was not differentially methylated due to glyphosate injury (Dataset S3). Several biological processes and molecular functions were enriched in the GO terms linked to the DMR-associated genes (genes that contain DMRs in their coding sequence) (Figs. 4B and 4C; Dataset S6). The enrichment of genes associated with phosphate (including phosphorylation, phosphotransferase, protein kinase, and pyrophosphatase) was striking. This could be related to glyphosate induction of microRNAs that regulate phosphate transport pathways observed in maize (Żywicki et al., 2015). The glyphosate molecule is recognized by phosphate transporters (Denis & Delrot, 1993) and some cases of glyphosate resistance have been attributed to alterations in herbicide transport and sequestration (Sammons & Gaines, 2014).

As expected, the majority of the enriched GO terms were associated with genic DMRs in the CG context (Fig. S6; Dataset S6), and these included many of the terms related to phosphate. The GO terms enriched in the CHG and CHH contexts were distinct. In the CHG context, glyphosate herbicide injury was associated with alterations in methylation of cell wall-associated and hydrolase genes. Genes associated with signal reception and transduction were specifically enriched in this analysis for genic DMRs in the CHH context.

Formation and maintenance of DMRs is driven by changes in gene expression (Secco et al., 2015), and DMRs are sometimes associated with downstream changes in gene expression (Li et al., 2012) though often have no effect on gene expression (Kawakatsu et al., 2016). We, therefore, cannot conclude whether glyphosate is directly altering the methylation patterns of the genes associated with phosphate metabolism and the shikimate pathway or if glyphosate is altering transcription of associated genes that in turn modify the methylation states of the underlying genomic loci. Previous work showed that glyphosate caused only minor changes in immediate transcriptional activity in A. thaliana eco. Columbia (Das et al., 2010), but induced major effects on the transcriptome of A. thaliana eco. Lansberg erecta (Faus et al., 2015). Future work including paired transcriptome analyses are required to correlate methylation patterns with gene expression.

Next, we asked whether the DMRs induced by glyphosate treatment are indicative of a general stress response or potentially at least partially glyphosate specific. We compared the 3,680 non-redundant DMR-associated genes in either the 5% or 10% treatment group to DMR-associated genes that were previously identified following salicylic acid (SA) treatment, which mimics biotic stress (Dowen et al., 2012), or phosphate starvation stress (Yong-Villalobos et al., 2015). Even considering the use of different approaches and cutoffs for calling DMRs, we hypothesized that the three stresses would induce overlapping methylome alterations as a result of the plants’ general response to environmental stresses and would therefore share some DMR-associated genes. Only 13 genes were identified as being differentially methylated following all three stress treatments (Fig. 4D; Dataset S7). These 13 DMR-associated genes are top candidates for genes that are epigenetically regulated as part of a general stress response. An additional 109 and 130 glyphosate-induced DMRs overlapped with DMRs associated with biotic stress and phosphate starvation stress, respectively. In total, 3,428 (93%) of the glyphosate-induced DMR-associated genes were unique to glyphosate exposure. While the number of overlapping DMR-associated genes among the three analyzed stresses seems small, they represent significantly more than would occur by random selection from all A. thaliana genes (see Methods). Moreover, χ² tests of independence of DMR-associated genes in three stresses show that stress responses are not independent of one another (χ² = 179.39, df = 4, p = 1.01e⁻³⁷), and there is indeed a higher than expected number of genes involved in response to all three stresses, indicating the existence of a common methylome reprogramming pathway in Arabidopsis regardless of the stressors. Nevertheless, the majority of identified DMRs appear specific to glyphosate injury.

Conclusion

This work identifies a large set of genes and other genomic regions epigenetically regulated in response to glyphosate herbicide injury. Determining the extent to which stresses induce patterned versus random epimutations is critical for understanding the role of DNA methylation in plant adaptation to stresses. We favor the hypothesis that the methylome changes are, at least partially, stress-specific over the hypothesis that methylome alterations are universally random because: (1) A majority of the DMRs exhibit dose-sensitive response patterns; (2) all analyzed glyphosate-induced DMRs were identified in four independent biological replicates; and (3) DMRs were enriched in gene pathways known to be affected by glyphosate exposure. Identification of transgenerationally stable DMRs and confirmation of specific DMRs directly correlated with glyphosate resistance will further clarify the role of methylome reprogramming in the evolution of herbicide resistance.

Supplemental Information