Modeling Interactions between Transposable Elements and the Plant Epigenetic Response: A Surprising Reliance on Element Retention

Modeling Interactions between Transposable Elements and the Plant Epigenetic Response: A... Transposable elements (TEs) compose the majority of angiosperm DNA. Plants counteract TE activity by silencing them epigenetically. One form of epigenetic silencing requires 21–22 nt small interfering RNAs that act to degrade TE mRNA and may also trigger DNA methylation. DNA methylation is reinforced by a second mechanism, the RNA-dependent DNA methylation (RdDM) pathway. RdDM relies on 24 nt small interfering RNAs and ultimately establishes TEs in a quiescent state. These host factors interact at a systems level, but there have been no system level analyses of their interactions. Here, we define a deterministic model that represents the propagation of active TEs, aspects of the host response and the accumulation of silenced TEs. We describe general properties of the model and also fit it to biological data in order to explore two questions. The first is why two overlapping pathways are maintained, given that both are likely energetically expensive. Under our model, RdDM silenced TEs effectively even when the initiation of silencing was weak. This rela- tionship implies that only a small amount of RNAi is needed to initiate TE silencing, but reinforcement by RdDM is necessary to efficiently counter TE propagation. Second, we investigated the reliance of the host response on rates of TE deletion. The model predicted that low levels of deletion lead to few active TEs, suggesting that silencing is most efficient when methylated TEs are retained in the genome, thereby providing one explanation for the large size of plant genomes. Key words: methylation, silencing, transposable elements, genome size. Introduction Despite the obvious evolutionary success of TEs, the plant Angiosperm genomes vary >1,000-fold in size, and this var- host checks their proliferation. The two entities engage in a iation correlates strongly with transposable element (TE) con- continuous arms-race, where TEs seek to proliferate and the tent. For plant species with small genomes, like Arabidopsis host attempts to control them (Lisch and Slotkin 2011). In thaliana or Brachypodium distachyon, DNA derived from TEs fact, most—but not all (Li et al. 2010)—TEs are epigenetically constitute 20–30% of the genome (AGI 2000; IBI 2010). silenced under normal conditions (Lisch 2009). The plant host Species with larger genomes have commensurately larger exerts this control by suppressing TE activity both before and proportions of TE-derived DNA. For example, TE-derived after transcription. Posttranscriptional modification relies DNA represents >85% of the barley (Hordeum vulgare)and chiefly on RNAi that recognizes and degrades TE mRNA pro- maize (Zea mays ssp. mays) genomes (Wicker et al. 2004; duced by RNA polymerase II (Pol II). Degradation requires as- Schnable et al. 2009). When one considers that the average sociated factors like RNA-polymerase 6 (RDR6), which size of a diploid angiosperm genome is similar to that of barley converts single-stranded to double-stranded RNA (dsRNA); genome, at 6400 Mb, then it is clear that most extant plant the Dicer-like proteins DCL2 and DCL4 that cleave dsRNAs DNA is derived from TEs (Tenaillon et al. 2010). to produce 21 and 22 nucleotide (nt) small interfering RNAs The Author(s) 2018. Published by Oxford University Press on behalf of the Society for Molecular Biology and Evolution. This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/), which permits non- commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact journals.permissions@oup.com Genome Biol. Evol. 10(3):803–815. doi:10.1093/gbe/evy043 Advance Access publication February 21, 2018 803 Downloaded from https://academic.oup.com/gbe/article-abstract/10/3/803/4895092 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Roessler et al. GBE (siRNAs); and the Argonaute1 (AGO1) protein that guides the RdDM pathway, which preferentially acts on TEs that have siRNAs to mRNAs for cleavage (Fultz et al. 2015). already been targeted for silencing. An important feature of Presumably, 21–22 nt siRNAs can prime multiple cycles of these 24 nt siRNAs is that they can act in trans to guide the mRNA cleavage, but they may have another important func- methylation of TEs that have similar sequence characteristics tion, which is to initiate transcriptional silencing (Nuthikattu to the original TE template (Slotkin et al. 2005; Teixeira et al. et al. 2013; McCue et al. 2015). Hence, 21–22 nt siRNAs can 2009; Ito et al. 2011; Ye et al. 2012; Fultz et al. 2015). Under be seen as dual-purpose, because they are involved in post- this process, 24 nt siRNAs may constitute a kind of “immune transcriptional silencing and also because they initiate DNA memory” that act as a buffer against the possibility of TE methylation (Cuerda-Gil and Slotkin 2016). activity (Fultz et al. 2015). If true, this implies that the strength Transcriptional silencing is achieved through epigenetic of the host epigenetic response is related to the number of modifications like DNA methylation, histone modifications, similar TEs in the genome that have already been silenced. and shifts in nucleosome positioning (Bernatavichute et al. Yet, no studies have explored the potential codependence 2008; Chodavarapu et al. 2010). The first of these, DNA between TE copy numbers and the strength of the host methylation, relies on the RNA-directed DNA methylation response. (RdDM) pathway. RdDM begins when the plant-specific Our final systems-level question concerns a separate pro- RNA polymerase Pol IV transcribes a TE. The resulting single- cess that occurs in cells associated with (but not part of) the stranded RNA is processed into 24 nt siRNAs by RDR2 and germline. In cells such as the pollen vegetative nucleus (Slotkin DCL3, two homologs that are distinct from those employed et al. 2009), some TEs are actively demethylated, expressed, in RNAi. Ultimately, the 24 nt siRNAs guide protein complexes and utilized to produce 21–22 nt siRNAs. These siRNAs are to homologous DNA sequences that are then targeted for then transported to the germline, where they presumably cytosine methylation. Once DNA methylation is established, contribute to stable TE silencing across generations (Slotkin at least two mechanisms act to maintain it. The first is a pos- et al. 2009; Ibarra et al. 2012; Martınez et al. 2016; Martinez itive feedback loop: Pol IV and Pol V, the RNA polymerases and Ko ¨ hler 2017). But what is the systems-level benefit of this involved in RdDM, preferentially act on methylated DNA (Law additional step in the host response, given that there are al- et al. 2013; Johnson et al. 2014), thereby reinforcing silencing ready at least two overlapping pathways dedicated to silenc- (Panda and Slotkin 2013). The second is the maintenance of ing TEs and also that symmetric DNA methylation is typically symmetric CGand CHG(where H¼ A, C, or T) methylation inherited faithfully? during DNA replication and cell division (Law and Jacobsen Here, we address these questions by building a model 2010). Although the switch from RdDM to maintenance is of host: TE interactions based on ordinary differential not well understood (Panda and Slotkin 2013), onceaTEis equations (ODEs). ODE models have been used widely to targeted for DNA methylation the host genome employs study biological phenomena that range from population feedbacks to ensure that the TE reaches and maintains a qui- growth (Malthus 1798), to predator–prey interactions escent state. (Volterra 1926), to the dynamics of viral infection and re- Numerous molecular studies have characterized the RNAi production (Perelson 2002). ODE models have also studied and RdDM pathways (reviewed in Law and Jacobsen 2010; the interactions between TEs and the host response Fultz et al. 2015; Matzke et al. 2015). These have been com- (Abrusan and Krambeck 2006), but without a focus on plemented by evolutionary studies showing that small RNAs plants and with few details of host response mechanisms. are used for TE defense across both prokaryotes and eukar- Our model includes proxies for RNAi, RdDM, and addi- yotes (Blumenstiel 2011)and that most RNAi and RdDM com- tional factors like TE propagation and TE deletion. We ponents are present in early land plant lineages (Huang et al. study properties of the model but also estimate reasonable 2015; Ma et al. 2015; Zhang et al. 2015; Tsuzuki et al. 2016). biological parameters by fitting the model to biological However, several important questions remain about systems- data, specifically from the study of the accumulation of level interactions between TEs and their plant hosts. One ma- the Evade element in an A. thaliana inbred line (Mari- jor question is why the host relies on two mechanisms—that Ordonez et al. 2013). Given these parameter estimates, is, RNAi and RdDM—to silence TEs. Presumably both path- we explore dynamics of the model and address systems- ways are capable of silencing; they are thus overlapping and level questions about host: TE interactions. We focus on potentially redundant. Both require the production of myriad three sets of questions: 1) Are both pre- and posttranscrip- polymerases, methylases and small RNAs and therefore must tional silencing necessary to control TEs? If not, what ad- have some energetic cost (Bousios and Gaut 2016). Why, vantage is gained by having two mechanisms? 2) Given then, are two pathways maintained? One working hypothesis that methylated TEs may be an important source of im- is that they act synergistically, but this hypothesis has yet to be mune memory, does TE deletion affect the dynamics of explored. the host response? And, finally, 3) What is the added ben- A second major question concerns 24 nt siRNAs. As men- efit of a third mechanism for generating 21–22 nt siRNAs tioned earlier, 24 nt siRNAs are predominantly produced by in the germline? 804 Genome Biol. Evol. 10(3):803–815 doi:10.1093/gbe/evy043 Advance Access publication February 21, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/3/803/4895092 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Modeling Interactions between Transposable Elements and the Plant Epigenetic Response GBE X X 2 2 Materials and Methods sqEr ¼ ðE  O Þ þ w  ðE  O Þ : CN CN Exp Exp Equilibria and Stability In this formula, E and O are the expected and observed CN CN Given our ODE model, its equations and its parameters copy number, respectively. The expected copy number was (see Results), we found equilibria by solving for active TEs defined as the sum of aTEs and sTEs obtained from the model. (aTEs), silenced TEs (sTEs) and siRNAs when all equations E and O are the expected and observed values, respec- Exp Exp were equal to zero. The first, trivial equilibrium point was tively, for relative expression. aTE ¼ sTE ¼ siRNA ¼ 0. To derive the stability of this eq eq eq The expected relative expression for generation n was equilbrium, we calculated the Jacobian matrix for the obtained from the model by taking the total expression in ODEs around that equilibrium, which provided: generation 8, which is equal to v multiplied by the aTE copy 2 3 vp  d 00 number at generation 8, and comparing that to the total ex- 6 7 pression at generation n, which is equal to v multiplied by the d 0 6 7 JðÞ 0; 0; 0¼ : 6 7 TE number of aTEs in generation n. Note, however, that our 4 5 measure of relative expression may not correspond perfectly ev 0 1 to that from Mari-Ordonez et al. (2013), because the empir- ical data on relative expression actually compares two genes The resulting characteristic equation is: (Evade and ACT2) within each generation and also because qRT-PCR can be inaccurate, especially when it is used as a detðÞ JðÞ 0; 0; 0 k  I¼ 0 TE 3 3 2 ratio (of ACT2 vs. Evade expression). In the square error (sqER) ¼ k þ vpk þ k vpd  vp þ d equation, we assigned w a weight of 40 to reflect the mag- vpd  d ; nitude of difference in the empirical data, because copy num- (1) ber reached 40 and relative expression plateaued at 1 (fig. 2A). where the solutions of this equation are the eigenvalues. We used a Monte Carlo approach to estimate fitted The equation clearly communicates that stability depends parameters. In this approach, all seven parameters were ini- on a complex relationship among v, p,and d but only on tialized with randomly drawn values from a uniform distribu- these parameters. The critical eigenvalue (i.e., the one tion between 0 and 1, except for v, which was ranged that crosses zero) is in fact (pv – d), and hence the sta- between 0 and 20. We also imposed the constraint that bility of this equilibrium is controlled by this compound pþ e  1.0. Given initial parameters, the sqEr was calculated parameter. See supplementary text, Supplementary as above. A single parameter was then altered, with a step Material online, for additional details based on a rescaled size between 0.1 and 0.1 for all parameters (except v where model. step size was between 1.0 and 1.0). The sqEr was calculated The second equilibrium point is shown in equations (3) and and the iteration moved forward only if sqEr > sqEr ;oth- n nþ1 (4) (see Results) for aTE and sTE ; the corresponding equa- eq eq erwise a new step size would be calculated. All the parame- tion for siRNA is: eq ters (in the following order: v, l, l,p,i,r, e,and d)were iterated through 100 times with 50 steps for each parameter, siRNA ¼ (2) eq r i ed until the final fitted parameters were found with the smallest ð pÞ: sqEr for each run. The initialization and iteration of all param- eters was performed >10,000 times; the lowest sqEr across all We also examined the Jacobian matrix and eigenvalues 10,000 runs was used to define the fitted parameters. We to study stability for this equilibrium point. Because the note, however, that other fitted data sets with low sqEr values stability equation was complex, we analyzed the rescaled produced similar model dynamics (supplementary fig. S6, model for further insights into the stability of this nontrivial Supplementary Material online). equilibrium (see supplementary text, Supplementary Material online). Running the ODE Model Fitted Parameters The ODE model was run using odeint from the scipy.integrate We obtained the data from Mari-Ordonez et al. (2013) by package (https://docs.scipy.org/doc/scipy/reference/integrate. loading their figure 3a onto WebPlotDigitizer (https://autome- html, last accessed February 28, 2018) and python (v2.6.6). ris.io/WebPlotDigitizer/, last accessed February 28, 2018). To Figure 1A was made with draw.oi (https://www.draw.io/,last estimate model parameters that fit the empirical data, we accessed February 28, 2018); all other figures were made with used the sum of least squares method, based on the following R (v. 3.3.2). The heatmaps were made with heatmap2, from formula: the gplots library in R. Genome Biol. Evol. 10(3):803–815 doi:10.1093/gbe/evy043 Advance Access publication February 21, 2018 805 Downloaded from https://academic.oup.com/gbe/article-abstract/10/3/803/4895092 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Roessler et al. GBE FIG.1.—A schematic of the model, with details provided in the text. The dashed arrow represents a step specific to cells that contribute to germline material. Table 1 24 nt siRNA is proportional to the number of silenced TEs Summary of Parameters and Their Fitted Estimates (sTEs). Furthermore, 24 nt siRNAs are considered to be trans-acting and thus may affect numerous TE insertions, in- Parameter Description Fitted cluding active elements. Overall, active TEs (aTEs) may be- Estimate come sTEs through 21–22 nt siRNAs, 24 nt siRNAs, or by a V Amount of Pol II mRNA expressed by active TEs 1.630 combination of both (fig. 1). P Proportion of mRNA that contributes to 0.340 transposition The model includes two additional parameters. The first is Proportion of mRNA that contributes to 0.051 TE deletion from the genome, which occurs at rate d for both 21–22nt siRNA production aTEs and sTEs. The second is the potential for the loss of I The rate at which 21–22nt siRNA initiate 0.062 silencing from TEs over time (e.g., through the loss of meth- methylation ylation), which we assume can lead to reactivation of TEs at R The rate at which 24nt siRNA reinforce 0.025 rate u.When u¼ 0, maintenance of silencing is perfect, but methylation silencing is not maintained when u¼ 1. D The rate of TE deletion per generation 0.161 The model is represented diagrammatically in figure 1 and U The rate of methylation loss per generation 0.000 consists of three differential equations: D The rate of degradation of 21–22nt siRNAs per 0.999 generation dðaTEÞ ¼ðÞ v  p  d  i  siRNA  r  sTE aTE þ u  sTE; dt dðsTEÞ Results ¼ðÞ i  siRNA þ r  sTE aTE ðÞ d þ u sTE; dt A Model of TE Propagation and Silencing Our model assumes that an active TE (aTE) begins as single dðsiRNAÞ ¼ e  v  aTE  d  siRNA: copy and expresses mRNA at rate v (fig. 1 and table 1). dt Among the produced mRNA, a proportion p is transposed into new genomic copies of the TE per host generation. The first equation describes the change in the number of Another proportion, e, of the TE mRNA is processed into aTEs over time; the second describes the change in the 21–22 nt siRNAs. Note that pþ e  1.0 under our model. number of sTEs over time, and the third monitors numbers We assume that the 21–22 nt siRNAs degrade at rate d and of 21–22 nt siRNAs over time. Although these three equa- initiate TE silencing at rate i. Initiation encompasses both post- tions represent our basic model, figure 1 includes a transcriptional silencing (RNAi) and the onset of methylation, dashed arrow representing a fourth process, the epige- following previous models (Nuthikattu et al. 2013; McCue netic remodeling of TEs in the germline. This process will et al. 2015). Finally, 24 nt siRNAs reinforce methylation at be incorporated after we first explore the dynamics of the rate r, representing RdDM. In our model, the amount of basic model. 806 Genome Biol. Evol. 10(3):803–815 doi:10.1093/gbe/evy043 Advance Access publication February 21, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/3/803/4895092 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Modeling Interactions between Transposable Elements and the Plant Epigenetic Response GBE vp We recognize that the model, as presented, is over- sTE ¼ (4) eq parameterized in a mathematical sense. For example, r ie pÞ the parameter v could be eliminated by redefining e and p. We retain the parameter definitions provided with siRNA given by equation (2) (see Materials and eq above throughout the main text, because we believe Methods). These two equations illustrate that aTEs and sTEs their meaning to be biologically intuitive. However, in have similar parameter dependencies. However, equilibrium the supplementary text, Supplementary Material online, values of sTEs depend more explicitly on v and p in the nu- we also present a mathematical treatment that includes merator than does the equilibrium values of aTEs. This is an parameter reduction, rescaling, and more extensive der- interesting observation because v and p are properties of ivations of the model’s analytical properties. We refer to aTEs; it drives home the point that equilibria copy numbers this supplement throughout the main text, where of sTEs relies intricately on the properties of their active coun- appropriate. terparts. The denominator of the two equations clearly indi- cates that increasing r tends to decrease both aTE and sTE . eq eq Model Equilibria Finally, the equations also hint at a complex relationship be- tween equilibrium copy numbers and d, because the latter Once a TE has invaded a host it has three possible fates: it appears twice in the denominator (and once in the nominator may fail to successfully invade and be lost completely; it may for sTEeq). As d increases, these appearances have opposite establish itself and reach an equilibrium number of copies effects on equilibrium values. over time; or it may expand in copy number unabated. An We studied these equilibria using an equivalent rescaled advantage of ODE models is that we can analytically solve model (see supplementary text, Supplementary Material on- the equilibrium points to understand TE invasion behavior line). Our analytical results provided additional insights and parameter dependence. We analyzed equilibria and the about the behavior of the model and particularly the stabil- stability of those equilibria. For these analyses we assumed ity of the nontrivial equilibrium. For example, the nontrivial u ¼ 0and d ¼ 1 for simplicity, but also because it is biolog- equilibrium is stable when (pv  d) is positive and unstable ically reasonable to assume both that maintenance of the when this is negative. Conversely, the trivial equilibrium is silenced state is strong (u ¼ 0), based on the conservation of stable when (pv  d) is negative, which implies that both symmetric methylation, and that siRNAs degrade rapidly equilibria exchange their stability for (pv  d)¼0. (d ¼ 1). We identified two equilibrium points in our system. The first is when there are no TEs and, hence, no 21–22 nt and Fitting the Model to Biological Data 24 nt siRNAs in the host. That is, the equilibrium points for the active copies (aTE It can be difficult to identify biologically reasonable parameter ), silenced copies (sTE ), and eq eq values for ODE models. To address this concern, we fitted the siRNA (siRNA ) are equal to zero. Stability around this eq model to biological data and then perturbed parameter values point provides information as to whether a TE will suc- cessfully invade the genome or be lost. We investigated separately to explore parameter dependencies and to assess effects on TE copy numbers. We fitted the model to experi- stability (see Materials and Methods; see eq. 1)and found mental data from the study of Mari-Ordonez et al. (2013), that it does not rely on any of the parameters associated who characterized the expression and transposition of a with epigenetic processes—that is, i, e,or r. Instead, sta- single-copy of the Evade retroelement that had become bility relies only on the parameters for TE expression, unmethylated in A. thaliana met1-mutant epigenetic recom- propagation, and deletion (v, p,and d; see also supple- mentary text, Supplementary Material online). Although binant inbred lines. By following two lines to generations 14 and 15, they showed that Evade was highly expressed until equation (1) is complex, the Jacobian matrix (see Materials generation 11 and 7, respectively, after which expression and Methods) suggests the intuitive notion that invasion proceeds when expression and propagation (v  p)out- plummeted precipitously, presumably due to host silencing. The number of Evade copies increased rapidly while its expres- competes deletion (d). sion was high, to a maximum of 40 copies after 11 and 7 Once a TE has established its presence in the host, it may increase in number until the second, nontrivial equilibrium generations. To fit our model to their data, we extracted information point (see Materials and Methods). The equilibrium points about Evade copy numbers and relative expression (see for aTEs and sTEs are given by: Materials and Methods). We focused on one inbred line aTE ¼ (3) eq (met1) from their study, because this was the only line for r ie pÞ; which data were sampled for consecutive generations: in total seven generations (from 8 to 14) since the reactivation of the single Evade element. We fitted the model to the Evade data Genome Biol. Evol. 10(3):803–815 doi:10.1093/gbe/evy043 Advance Access publication February 21, 2018 807 Downloaded from https://academic.oup.com/gbe/article-abstract/10/3/803/4895092 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Roessler et al. GBE TE copy number. If TEs were able to grow unabated, there would be an exponential increase at a rate of 0.554 (¼pv) TEs per generation. However, some transcripts are processed into 21–22 nt siRNAs (e ¼ 0.051) that silence TEs at rate i¼ 0.062. These 21–22 nt siRNAs degrade quickly for each host gener- ation (d ¼ 0.999), and therefore any new 21–22 nt siRNAs are not residual, but must be made from active TEs. Once initia- 0 5 10 15 0 5 10 15 tion of methylation has begun as part of i, reinforcement quickly takes hold at rate r¼ 0.025. Eventually, the number Generations Generations of sTEs increases and the number of aTEs decreases, so that aTEs total expression begins to decline. mTEs sTEs TE max As TEs become silenced, they have two fates under our Copy Number model: they can be deleted from the genome or become TE final active again due to loss of silencing (fig. 1). Since loss of si- lencing was very low (u¼ 410 ) inthe fittedparameter set, the main fate of sTEs is to be deleted (d¼ 0.16). As these quiescent TEs are lost, so is the source of reinforcing 24 nt 0 100 200 300 400 500 siRNAs. When reinforcement becomes unreliable, the host Generations loses epigenetic control, the subset of remaining aTEs propa- gate, and the phased cycle begins again. These cycles dissi- FIG.2.—(A) Model fit to the Evade data for total copy number (left) pate in amplitude until equilibria are reached at 20 total TE and relative expression (right). The empirical data from the Evade study are copies, with more sTEs (14) than aTEs (6) (fig. 2B). It is represented by circles; the whiskers indicate SD. The model results based on the fitted parameters (table 1) are represented by the solid line. important to note that the equilibrium is not necessarily static; (B) Long-term behavior of the model, based on the fitted parameters to it can be reached when equal numbers of TEs are created the Evade data. Arrows show TE and TE , which are defined in the max final versus deleted. text. Copy number refers to the summation of aTEs and sTEs. These phased interactions occur with the fitted parame- ters, but decaying oscillations in copy number also occur reg- ularly with other parameter combinations (see supplementary with a Monte Carlo approach that concurrently considered text and fig. S1, Supplementary Material online). Oscillating TE the total TE copy number (i.e., the combined total of aTEs and numbers are not, however, a necessary outcome of the model sTEs) and TE expression. Our set of fitted parameter values are (see examples below and supplementary text, Supplementary reported in table 1. These parameter values produced a good Material online). fit to the copy number data, and a curve of similar shape to the observed relative expression data over time (fig. 2A). (Note that our measure of expression is only a proxy for expression Examining Initiation (i) and Reinforcement (r) measured experimentally; see Materials and Methods.) We We have shown that the model can have complex, oscillating recognize that we have fitted a complex model to relatively dynamics based on parameters inferred from biological data. simple data and that our fitted parameters may represent one These parameters can be modified independently to explore of many potential reasonably fitting parameter sets (but see the importance of various processes. In this section, we assess Materials and Methods). They nonetheless provide a biologi- the effect of perturbing the system by varying either initiation cally plausible foundation for examining model behavior. (i) or reinforcement (r), or both, while holding the remaining parameters to the values estimated from the Evade data. We first set i¼ 0, and the result was both intuitive and trivial. With Model Behavior under Fitted Parameters i¼ 0 silencing never begins. Hence, the number of aTEs Given the fitted parameters, we explored host: TE dynamics trended upward at an exponential rate, with no resulting over 500 generations, monitoring numbers of aTEs, sTEs, and sTEs (fig. 3A). total TE copy number (¼aTEsþ sTEs) (fig. 2B). With these The effectofsetting r to zero was less straightforward, parameter values, the model produces oscillations of all three because i was > 0 and hence silencing was initiated. entities for 200 generations until it reaches an equilibrium. Without reinforcement, copy numbers no longer oscillated, The oscillations of aTEs and sTEs are somewhat out of phase but instead burst and rapidly reached a maximum for both with one another. We interpret these results as reflecting aTEs and sTEs. These copy numbers remained flat, implying a feedbacks in the epigenetic system. When a TE first invades steady state in which silencing was initiated by 21–22 siRNAs a host, the combination of expression (fitted value v¼ 1.63; and there was sufficient transposition to counteract TE dele- table 1) and propagation (p¼ 0.340) create an initial burst in tion. Under these parameter values, the steady state of sTEs 808 Genome Biol. Evol. 10(3):803–815 doi:10.1093/gbe/evy043 Advance Access publication February 21, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/3/803/4895092 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Number of TEs Copy Number 020 40 60 0 1020304050 Relative Expression 0.0 1.0 2.0 3.0 Modeling Interactions between Transposable Elements and the Plant Epigenetic Response GBE FIG.3.—Model behavior with the fitted values for all parameters but initiation (i) and reinforcement (r). (A) Graphs illustrate the effect of setting initiation and reinforcement parameters to zero for active TEs (left) and methylated TEs (right). In both graphs, the gray dashed lines represent the number of TEsbased on the fitted model parameters to the Evade data (see also fig. 2B). (B) Heat maps showing the TE (left) and TE (right) for the total copy number max final (¼ aTEsþ sTEs) based on varied values of initiation (y-axis) and reinforcement (x-axis), with copy number displayed in each cell. The dashed cell in each heat map represents the fitted values (table 1). was higher than that of aTEs (fig. 3A), as with the equilibria when r was low (e.g., r 0.1), the value of i had notable reached with fitted parameters (fig. 2B). effects on both TE and TE For example, when max final. If initiation by 21–22 nt siRNAs is sufficient to reach a r¼ 0.001, TE varied over two orders of magnitude as a max steady state and to control TEs, then what is the advantage function of i. Similarly, TE differed 33-fold when i ranged final of reinforcement by 24 nt siRNAs? Equations (3) and (4) show from 0.001 to 0.99 (fig. 3B). This relationship implies that that the steady-state TE copy numbers depend critically on the reinforcement can counter TE propagation efficiently, even values of r and i, since the denominators become very small when initiation of silencing is weak. This observation held where r and i are small. To explore their interdependencies, true when also adjusting for TE expression (v) and deletion we varied i and r across their parameter ranges and assessed (d)(supplementary fig. S2, Supplementary Material online). total copy numbers (¼aTEsþ sTEs). To help characterize effects, we focused on two descriptive statistics, TE and max The Effects of TE Deletion (d) TE (see fig. 2B) TE is the highest total TE copy number final . max achieved under a set of model parameters, and TE is the Theoretically, high TE deletion rates should be advantageous final total copy number after 5,000 generations, a point by which for the plant host, because they limit opportunities for trans- total aTE and sTE copy numbers have typically reached a position and consequent deleterious mutations. However, steady state. Our analyses show that when r0.5, any high amounts of TE deletion could have consequences for change in i had little effect on TE and TE , so long as immune memory, because quiescent TEs may be a major max final there was at least some initiation (fig. 3B). In contrast, source of trans-acting 24 nt siRNAs (Teixeira et al. 2009; Ito Genome Biol. Evol. 10(3):803–815 doi:10.1093/gbe/evy043 Advance Access publication February 21, 2018 809 Downloaded from https://academic.oup.com/gbe/article-abstract/10/3/803/4895092 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Roessler et al. GBE 0.001 0.01 0.1 0.5 0.99 0 500 1000 1500 Generations 0.001 0.01 0.1 0.5 0.99 FIG.5.—The effect of varying rates of TE deletion (d) on the final 0 500 1000 1500 number of silenced TEs (sTE), active TE (aTE) and total copy number Generations (aTEþ sTE). These calculations used the fitted values for all parameters but d. Note also that d begins at an arbitrarily low value of 0.001; when FIG.4.—Model behavior with the fitted values for all parameters but d¼ 0, the number of sTEs diverges. TE deletion (d), which is varied from 0.001 to 0.99. et al. 2011; Fultz et al. 2015). Hence, high deletion rates may silencing. This supports our observation, based on equilibrium adversely affect the epigenetic response. To illustrate the ef- equations (eqs. 3 and 4), that deletion plays a complex role in fect of deletion on TE copy numbers, we varied the deletion determining aTE and sTE . eq eq parameter d from 0.001 to 0.99 (fig. 4), while holding the remaining parameters to their fitted values (table 1). The model produced four noteworthy results. First, when Additional Parameters TE deletion was very low (d¼0.001), aTEs burst quickly to high copy number (30). After peaking at a total copy num- We also varied values of expression (v), propagation (p), and ber of 80, all TEs were silenced and the population of sTEs loss of silencing (u), while the remaining parameters were declined slowly over time, reflecting the low rate of deletion held at their fitted values. The parameter v was arbitrarily (fig. 4). Throughout this process, there were no aTEs after the ranged between 0 and 5. The chief effect of this range was initial burst. Second, when d increased (0.01 d < 0.5), the on the amplitude and periodicity of TE oscillations. Higher system generated oscillations in the number of aTE and sTEs. expression levels led to more dramatic copy number oscilla- The amplitude, frequency, and equilibrium values (i.e., TE ) tions (supplementary fig. S3, Supplementary Material online). final varied with d. Note that the running average of sTEs exceeded Importantly, at low parameter values (e.g., v0.5) TEs either that of aTEs for these parameter values (fig. 4). Third, when TE did not invade the genome or were maintained at very low deletion was at intermediate levels (d¼ 0.5), aTEs reached a copy numbers (<5 total TEs) over the long term. Varying steady state, but there were very few sTEs. Finally, when the p produced results similar to varying v (supplementary fig. rate of TE deletion was very high (d¼ 0.99), all TEs were re- S4, Supplementary Material online). Increasing p did, how- moved from the genome. ever, tend to lead to higher TE and average copy num- max To further illustrate these dependencies on d,we plotted bers relative to the parameter values we explored for v TE for sTEs, aTEs, and all TEs as a function of d (fig. 5). (supplementary fig. S3, Supplementary Material online). final Overall, these results convey a somewhat counterintuitive This was presumably because there is a trade off with v; idea: if the goal is to have few aTEs, then it is beneficial either as it increases, so does the production of 21–22 nt siRNAs, to have dramatically high rates of TE deletion (e.g., d¼ 0.99) which then potentially affect RNAi. Propagation (p), on the or to have such low (e.g., 0.01–0.1) deletion rates that a res- other hand, contributes only to the proliferation of more ervoir of sTEs is preserved and contributes to reinforcement of TEs. Note that low levels of propagation (p < 0.25) resulted 810 Genome Biol. Evol. 10(3):803–815 doi:10.1093/gbe/evy043 Advance Access publication February 21, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/3/803/4895092 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Silenced TEs (sTEs) Methylated TEs (mTEs) Active TEs (aTEs) 020 40 60 80 010 20 30 40 Modeling Interactions between Transposable Elements and the Plant Epigenetic Response GBE a level proportional to the number of sTEs that were deme- Pollen Reactivation No Pollen Reactivation thylated in the companion cells. That is, dsiRNA ¼ e  v ðÞ aTE þ sTE d  siRNA: dt This equation is represented by the dotted arrow in figure 1. We evaluated the effects of this additional process on the 0 100 200 300 400 500 system with fitted parameter values. The effects were consis- Generations tent: it decreased TE ,TE and the periodicity of copy max final number oscillations (fig. 6). Thus, this additional process yields FIG.6.—TE reactivation in pollen. The black line is based on the model notable decrements in TE copy numbers. with fitted parameters (no pollen reactivation); the dashed line is using the same parameters but including additional feedback for pollen guard cells (pollen reactivation). Both lines indicate total copy numbers Discussion (¼aTEsþ sTEs). The additional mechanism in pollen guard cells is denoted In this study, we have devised an ODE model to examine the by the dashed arrow in figure 1. systems dynamics of TE propagation within the context of the epigenetic response of a plant host (fig. 1). Although there are in no invasion. Hence, TEs cannot invade if expression or clear limitations to our approach, the model has produced at propagation is low. least four fundamental insights. The first is the prediction of Our model also assumes a process of silencing loss (u), for oscillating copy numbers typified by a burst of TE activity, which the most likely example is methylation loss. followed by silencing, deletion, and then reactivity. Despite Methylation loss is known to be low based on empirical these oscillations, the system often reached equilibrium copy data because symmetric methylation is typically maintained numbers (fig. 2). Second, our model emphasizes the impor- faithfully through cell division (Becker et al. 2011). Indeed, tance of reinforcement by RdDM-like processes, because it our fitted parameter estimate was u ¼ 410 , suggesting buffers potential upstream inefficiencies in the initiation of that a very low amount of sTEs become aTEs due to, for silencing (fig. 3). Third, we show that these outcomes are example, leaky maintenance of symmetrical methylation. linked to the rate of TE deletion. Somewhat nonintuitively, Overall, we found that varying the u parameter had little the model predicts that either low or very high levels of dele- effect on model behavior at parameter values < 0.01 (sup- tion lead to more efficient control of the number of aTEs plementary fig. S5, Supplementary Material online). This (figs. 4 and 5). Finally, we show that demethylation within implies that variation in spontaneous demethylation rates germline cells reinforces host defenses by dampening TE is likely to have few effects on the dynamics of host: TE bursts and lowering steady-state copy numbers (fig. 6). interactions unless u varies by several orders of magnitudes Below, we first discuss the caveats of our ODE model before from our fitted estimate. placing our insights into the context of plant genome struc- ture and evolution. TE Reactivation Dampens TE Oscillations Caveats Finally, we incorporated an interesting biological observa- tion—that is, the fact that TEs are activated in some repro- Every model has limitations, and ours is no exception. One ductive tissues, ostensibly to ensure the transmission of a important consideration is that our biological knowledge of complement of siRNAs to egg and sperm (Slotkin et al. the host response is incomplete. For example, the details of 2009; Ibarra et al. 2012; Martınez et al. 2016; Martinez and the initiation of methylation are not yet clear, because there Ko ¨ hler 2017). TEs are known, for example, to be demethy- are at least two competing (but likely nonexclusive) hypothe- lated and reactivated in the pollen vegetative nucleus, which ses as to how the host transitions from RNAi to the RdDM accompanies the sperm cell, but does not contribute DNA response (Mari-Ordonez et al. 2013; Nuthikattu et al. 2013; to the fertilized zygote. The reactivated TEs are sources McCue et al. 2015). Furthermore, some aspects of the host of 21–22 nt sRNAs that are transported to the sperm and pre- response have not been included in our model, such as recent sumably target silencing of TEs in the zygote (Slotkin et al. discoveries that 18–22 nt tRNA fragments (Martinez et al. 2009). The net effect of this process is to increase the numbers 2017; Schorn et al. 2017) and some miRNAs (Creasey et al. of 21–22 nt siRNAs in germline cells; these 21–22 nt sRNA 2014) may interfere with TE replication and propagation. originate not only from aTEs but also from sTEs (see below). However, these additional host mechanisms fit relatively easily We added this mechanism to our model with an equation in our model, because they would likely affect conversion (e) that increases the number of 21–22 nt siRNAs in the system at andinitiation(i)(fig. 1). In this sense, our model already Genome Biol. Evol. 10(3):803–815 doi:10.1093/gbe/evy043 Advance Access publication February 21, 2018 811 Downloaded from https://academic.oup.com/gbe/article-abstract/10/3/803/4895092 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Copy Number 0 1020 3040 Roessler et al. GBE implicitly accounts for some exciting new findings, but other Supplementary Material online). These results likely reflect our new insights may require model modifications. reliance on data from a study in which silencing occurred Another limitation is that we have studied the invasion of rapidly (Mari-Ordonez et al. 2013), but there is other experi- only one TE family. In reality, plant genomes harbor a multi- mental evidence that host defenses react quickly to silence tude of TE types that may interact with each other and also active TEs within a few host generations (Teixeira et al. 2009; vary with respect to the host response. For example, some but Fultz and Slotkin 2017), perhaps even more quickly than the not all TE families in A. thaliana are recognized by endoge- host response to Evade. nous miRNAs (Creasey et al. 2014), and short, nonautono- It is interesting to note that these experimental studies mous DNA elements are methylated less efficiently than contradict numerous genome-wide analyses, which suggest longer, autonomous elements (Hollister and Gaut 2009), per- that TE families experience massive bursts lasting thousands haps in part due to biases in genomic location (Zemach et al. or even millions of years (Piegu et al. 2006; Schnable et al. 2013). Finally, we have used only one data set to fit the 2009; Bousios et al. 2012; Daron et al. 2014). One likely ex- model, which followed the invasion of the Evade TE for a planation for this incongruence may be the difficulty of re- short period of few host generations (Mari-Ordonez et al. solving the occurrence of multiple rounds of episodic bursts 2013). Therelianceon Evade reflects the fact that very few within the expanded timeframes reported by the genome- studies have monitored the copy number and expression of wide studies. Limited resolution may be due to technical issues TEs within a plant genome over time, particularly beginning related to in silico TE identification, accurate age estimation, from recent invasion or reactivation. In short, we recognize and perhaps even heterogeneous rates of TE sequence loss the limitations of the empirical data, but they nonetheless and decay across the genome (Tian et al. 2009). No matter allow a glimpse into model behavior under relevant parame- the cause, the apparent gaps between experimental and ter values. genome-wide studies deserve further thought and consider- ation. Longer term experimental studies that monitor TE copy numbers over time and under different stress conditions Invasion and Oscillations would certainly be welcome contributions to our empirical How long does it take to silence TEs in vivo? Our understand- understanding of host: TE interactions. ing of the duration and intensity of TE amplification bursts remains limited (Bousios and Gaut 2016). In order to be si- Equilibria lenced, a TE must first invade. Based on our model and anal- yses of the stability of the first equilibrium point (where Another question is whether TEs reach long-term equilibria aTEs¼ sTEs¼ siRNA¼ 0), invasion depends on expression within a genome. In our model, the oscillations often reduce (v), propagation (p), and deletion (d) but not on downstream in intensity over time to reach a steady state (figs. 2 and 4 and properties of the host response, such as conversion of TE supplementary figs. S3–S5, Supplementary Material online). transcripts to 21–22 nt siRNAs (e), initiation (i), and reinforce- In this equilibrium, sTEs are found in higher numbers than ment (r). Put simply, pv needs to outpace d for a TE to suc- aTEs whenever (pv)/d > 2(eqs. 1 and 2 and supplementary cessfully invade the host. We also investigated invasion by text, Supplementary Material online). modifying v and p from the fitted parameter values (table 1); Our ODE-based approach regularly predicts two phases of invasion did not occur when expression or propagation were host: TE dynamics: one shaped by oscillating changes in TE low (v < 0.5, supplementary fig. S3, Supplementary Material numbers, and another characterized by an equilibrium. TE online; p < 0.25, supplementary fig. S3, Supplementary evolution has been modeled extensively with population ge- Material online). netic approaches (Charlesworth and Charlesworth 1983; Assuming a TE invades successfully, it has the potential to Charlesworth et al. 1994; Brookfield 2005; Le Rouzic and increase rapidly in copy number. Under our model, we found Deceliere 2005), and the basic models predict that TEs reach that copy numbers often oscillated before reaching an equi- steady-state copy numbers after the first TE invasion through librium (e.g., fig. 2). Mathematically, the prevalence of oscil- either a transposition-selection or transposition-deletion equi- lations is related to the value of the expression ðv  e  iÞ=r (see librium. In other words, they do not predict oscillations prior to supplementary text, Supplementary Material online). an equilibrium. In contrast, some studies have expanded their Oscillations tend to occur when TE expression is low (v), the models to include TE sequence evolution or competition be- proportion of 21–22 nt siRNAs is low (e), initiation (i) is low or tween TEs, and these often predict oscillations in TE copy reinforcement (r)is high (see supplementary text, numbers (Le Rouzic and Capy 2006; Le Rouzic, Boutin, Supplementary Material online). et al. 2007; Le Rouzic, Dupas, et al. 2007). For example, Le Under many parameter values explored in this work, the Rouzic, Boutin et al. (2007) investigated host–parasite inter- maximum duration of a TE burst lasts for only a few dozen actions between autonomous TEs and their nonautonomous generations before they are temporarily silenced and decrease counterparts, and they found oscillations in copy numbers in copy number (figs. 2 and 4 and supplementary figs. S3–S5, between both entities. Notably, the oscillations continued 812 Genome Biol. Evol. 10(3):803–815 doi:10.1093/gbe/evy043 Advance Access publication February 21, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/3/803/4895092 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Modeling Interactions between Transposable Elements and the Plant Epigenetic Response GBE indefinitely; an equilibrium was rarely reached unless there parameter under the Evade model (figs. 4 and 5). If the were very low mutation rates and few adaptive TE insertions. goal is simply to rid the genome of TEs, the most efficient Le Rouzic, Boutin et al (2007) and also Brookfield (2005) have method is to have a very high d (> 0.5) that removes all aTEs argued that equilibria are reached under conditions that are and sTEs. However, deletions are mediated by ectopic recom- probably unrealistic for in vivo TEs. This is because the param- bination and illegitimate recombination (Devos et al. 2002) eters that affect TE dynamics such as selection, transposition, that may introduce a substantial fitness cost due to the po- and deletion are likely to change at faster rates than the time tential for catastrophic mutations (Langley et al. 1988). required to reach an equilibrium. Our model does not include Assuming that high ectopic recombination carries an unac- autonomous and nonautonomous TEs, nor does it allow per- ceptable fitness cost, our model suggests that the next best turbations in subsequent generations. Yet, the focus on active solution to limit the number of aTEs is to have very low rates and sTEs may mimic some characteristics of host–parasite of TE deletion (d 0.01). relationships and may contribute to our observed oscillating Ourargumentis thatthe retention ofsTEs may benefit the dynamics. We must caution, however, that our model is not host by boosting immune memory. In theory, this immune explicitly evolutionary, because it does not consider fitness or memory provides a defense against the invasion of new TEs population variation. that have sequence homology to existing genomic TEs (Fultz and Slotkin 2017) and also against TEs that have escaped silencing and need to be resilenced. Two interesting features The Importance of Overlapping Mechanisms of acquired immune memory are that it is energetically ex- Why do plants maintain two overlapping and energetically pensive but also maintained under frequent cycles of reinfec- costly pathways (RNAi and RdDM) to silence TEs? Here, i tion (Best and Hoyle 2013). Under the parameter values encompasses posttranscriptional silencing and the initiation explored with our model, the system usually reaches a steady of methylation, and r represents RdDM (fig. 1). Our results state in which the copy number of aTEs is >0. To the extent show that only a small amount of i is needed to begin silenc- that these dynamics reflect reality, a nonzero equilibrium of ing of an unrecognized TE, but r is necessary to counter prop- aTEs defines a system in which reinfection is not merely fre- agation efficiently. For example, the host maintains TE copy quent but constant. This observation may explain one feature numbers at low levels even when i is inefficient (e.g., of the selective pressure to maintain RdDM-like mechanisms, i¼ 0.001), so long as r reinforces silencing by a value of even though it seems as if most TEs within plant genomes are r 0.1 (fig. 3B). RNAi is clearly not as efficient at limiting TE effectively silenced. There is also a conjecture that “zombie” copy numbers when there is no RdDM, yet it is essential for TEs are maintained in the genome in order to produce siRNAs silencing TEs (fig. 3A). Hence, to the extent the model is cor- that boost immune memory and can trigger the trans-silenc- rect, it implies that plants must have RNAi to start the process ing of active relatives (Lisch 2009). Indirect in silico evidence of silencing, but RdDM vastly enhances host control over TEs. for the existence of zombie TEs has been recently uncovered The inclusion of another, apparently overlapping mecha- in maize (Bousios et al. 2016). nism—that is, the active demethylation of TEs in cells that Finally, if low rates of TE deletion are somehow beneficial contribute siRNAs to germline cells—further enhances host to the host response, this process could drive genome size silencing (fig. 6). increases over evolutionary time, because each new TE infec- Our data are consistent with the argument that 24 nt tion or TE reactivation adds copies that are silenced, retained, siRNAs are important for buffering TE activity, even though and not quickly deleted. We also note that this is unlikely to be they seem unnecessary because most methylation is main- a run-away process, because there is evidence for selection on tained independently of RdDM in heterochromatic regions genome size (Diez et al. 2013; Bilinski et al. 2017), especially (Zemach et al. 2013). In fact, it was recently shown that these when genome size gets too large (Knight et al. 2005). heterochromatic regions also produce 24 nt siRNAs, albeit to a Nonethless, our model offers a partial explanation for the smallerextent(Li et al. 2015). These findings are consistent high TE contents and sizes of plant genomes. with the idea that 24 nt siRNAs may act as immune memory (Fultz et al. 2015; Fultz and Slotkin 2017), based on evidence that they may play a key role in suppressing reactivated TEs Future Directions (Teixeira et al. 2009; Ito et al. 2011; Fultz et al. 2015; Fultz and This is the first study to explicitly incorporate features of the Slotkin 2017). plant host response into a quantitative model of host: TE dy- namics. We view this model as a foundation for further exten- The Curious Case of TE Deletion sions that will continue to elucidate important features of If 24-nt siRNAs act as a source of immune memory, then the host: TE interactions. One promising avenue will be to extend retention of sTEs may be a benefit to the host, because they our model to include populations, genetic drift, and fitness may be the template for 24 nt siRNA production. This rela- (Szitenberg et al. 2016), perhaps with a potential for rare tionship is implied by our analyses of the deletion (d) beneficial effects (Le Rouzic, Boutin, et al. 2007). Such an Genome Biol. Evol. 10(3):803–815 doi:10.1093/gbe/evy043 Advance Access publication February 21, 2018 813 Downloaded from https://academic.oup.com/gbe/article-abstract/10/3/803/4895092 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Roessler et al. GBE Charlesworth B, Charlesworth D. 1983. The population dynamics of trans- approach is likely to yield more realistic understandings of the posable elements. Genet Res. 42(01):1–27. evolution of host: TE interactions than are available at present. Charlesworth B, Sniegowski P, Stephan W. 1994. The evolutionary dy- It will also be illustrative to model multiple TE families, includ- namics of repetitive DNA in eukaryotes. Nature 371(6494):215–220. ing autonomous and nonautonomous elements, different Chodavarapu RK, Feng S, Bernatavichute YV, et al. 2010. Relationship length and classes of elements, and the possibility of extensive between nucleosome positioning and DNA methylation. Nature 466(7304):388–392. siRNA cross-homologies. Finally, an important future goal will Creasey KM, et al. 2014. miRNAs trigger widespread epigenetically acti- be to mimic reality by introducing stresses and perturbations vated siRNAs from transposons in Arabidopsis. Nature into the model. One potential example of a perturbation is 508(7496):411–415. polyploidy, which is thought to lead to epigenetic repattern- Cuerda-Gil D, Slotkin RK. 2016. Non-canonical RNA-directed DNA meth- ing (Matzke et al. 1999), but for which the causes remain a ylation. Nat Plants 2(11):16163. Daron J, Glover N, Pingault L, et al. 2014. Organization and evolution of mystery. transposable elements along the bread wheat chromosome 3B. Genome Biol. 15(12):546. Supplementary Material Devos KM, Brown JK, Bennetzen JL. 2002. Genome size reduction through illegitimate recombination counteracts genome expansion Supplementary data areavailableat Genome Biology and in Arabidopsis. Genome Res. 12(7):1075–1079. Evolution online. Diez CM, et al. 2013. Genome size variation in wild and cultivated maize along altitudinal gradients. New Phytol. 199:264–276. Fultz D, Choudury SG, Slotkin RK. 2015. Silencing of active transposable elements in plants. Curr Opin Plant Biol. 27:67–76. Acknowledgments Fultz D, Slotkin RK. 2017. Exogenous transposable elements circumvent identity-based silencing, permitting the dissection of expression- We are grateful for feedback from D. Seymour and dependent silencing. Plant Cell 29(2):360–376. D. Wodarz. K.R. is supported by the National Institute of Hollister JD, Gaut BS. 2009. Epigenetic silencing of transposable elements: Biomedical Imaging and Bioengineering, National Research a trade-off between reduced transposition and deleterious effects on neighboring gene expression. Genome Res. 19(8):1419–1428. Service Award EB009418. A.B. is supported by The Royal Huang Y, et al. 2015. Ancient origin and recent innovations of RNA po- Society (Award Numbers UF160222 and RGF\R1\180006). lymerase IV and V. Mol Biol Evol. 32(7):1788–1799. This work was also supported by NSF grant DEB-1655808 Ibarra CA, Feng X, Schoft VK, et al. 2012. Active DNA demethylation in to B.S.G. plant companion cells reinforces transposon methylation in gametes. Science 337(6100):1360–1364. IBI, International Brachypodium Initiative. 2010. Genome sequencing and Literature Cited analysis of the model grass Brachypodium distachyon.Nature Abrusan G, Krambeck HJ. 2006. Competition may determine the diversity 463:763–768. Ito H, et al. 2011. An siRNA pathway prevents transgenerational retro- of transposable elements. Theor Popul Biol. 70:364–375. transposition in plants subjected to stress. Nature 472(7341):115–119. AGI, Arabidopsis Genome Initiative. 2000. Analysis of the genome se- quence of the flowering plant Arabidopsis thaliana.Nature Johnson LM, Du J, Hale CJ, et al. 2014. SRA- and SET-domain-containing proteins link RNA polymerase V occupancy to DNA methylation. 408:796–815. Becker C, et al. 2011. Spontaneous epigenetic variation in the Arabidopsis Nature 507(7490):124–128. thaliana methylome. Nature 480(7376):245–249. Knight CA, Molinari NA, Petrov DA. 2005. The large genome con- straint hypothesis: evolution, ecology and phenotype. Ann Bot Bernatavichute YV, Zhang X, Cokus S, Pellegrini M, Jacobsen SE. 2008. Genome-wide association of histone H3 lysine nine methylation with 95:171–190. Langley CH, Montgomery E, Hudson R, Kaplan N, Charlesworth B. 1988. CHG DNA methylation in Arabidopsis thaliana. PLoS One 3(9):e3156. Best A, Hoyle A. 2013. The evolution of costly acquired immune memory. On the role of unequal exchange in the containment of transposable Ecol Evol. 3(7):2223–2232. element copy number. Genet Res. 52(3):223–235. Law JA, et al. 2013. Polymerase IV occupancy at RNA-directed DNA meth- Bilinski P, et al. 2017. Parallel altitudinal clines reveal adaptive evolution of genome size in zea mays. bioRxiv. Available from: https://www.biorxiv. ylation sites requires SHH1. Nature 498(7454):385–389. Law JA, Jacobsen SE. 2010. Establishing, maintaining and modifying DNA org/content/early/2017/07/13/134528, last accessed February 28, 2018. methylation patterns in plants and animals. Nat Rev Genet. Blumenstiel JP. 2011. Evolutionary dynamics of transposable elements in a 11(3):204–220. Le Rouzic A, Boutin TS, Capy P. 2007. Long-term evolution of transposable small RNA world. Trends Genet. 27(1):23–31. Bousios A, et al. 2012. The turbulent life of Sirevirus retrotransposons and elements. Proc Natl Acad Sci U S A. 104(49):19375–19380. Le Rouzic A, Capy P. 2006. Population genetics models of competi- the evolution of the maize genome: more than ten thousand elements tion between transposable element subfamilies. Genetics tell the story. Plant J. 69(3):475–488. Bousios A, et al. 2016. A role for palindromic structures in the cis-region of 174(2):785–793. Le Rouzic A, Deceliere G. 2005. Models of the population genetics of maize Sirevirus LTRs in transposable element evolution and host epi- genetic response. Genome Res. 26(2):226–237. transposable elements. Genet Res. 85(3):171–181. Le Rouzic A, Dupas S, Capy P. 2007. Genome ecosystem and transposable Bousios A, Gaut BS. 2016. Mechanistic and evolutionary questions about elements species. Gene 390(1–2):214–220. epigenetic conflicts between transposable elements and their plant hosts. Curr Opin Plant Biol. 30:123–133. Li H, Freeling M, Lisch D. 2010. Epigenetic reprogramming during vege- tative phase change in maize. Proc Natl Acad Sci U S A. Brookfield JF. 2005. The ecology of the genome – mobile DNA elements and their hosts. Nat Rev Genet. 6(2):128–136. 107(51):22184–22189. 814 Genome Biol. Evol. 10(3):803–815 doi:10.1093/gbe/evy043 Advance Access publication February 21, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/3/803/4895092 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Modeling Interactions between Transposable Elements and the Plant Epigenetic Response GBE Li S, et al. 2015. Detection of Pol IV/RDR2-dependent transcripts at the expansions in Oryza australiensis, a wild relative of rice. Genome genomic scale in Arabidopsis reveals features and regulation of siRNA Res. 16(10):1262–1269. biogenesis. Genome Res. 25(2):235–245. Schnable PS, Ware D, Fulton RS, et al. 2009. The B73 maize genome: Lisch D. 2009. Epigenetic regulation of transposable elements in plants. complexity, diversity, and dynamics. Science 326(5956):1112–1115. Annu Rev Plant Biol. 60:43–66. Schorn AJ, Gutbrod MJ, LeBlanc C, Martienssen R. 2017. LTR-retrotrans- Lisch D, Slotkin RK. 2011. Strategies for silencing and escape: the ancient poson control by tRNA-derived small RNAs. Cell 170(1):61–71.e11. struggle between transposable elements and their hosts. Int Rev Cell Slotkin RK, et al. 2009. Epigenetic reprogramming and small RNA silencing Mol Biol. 292:119–152. of transposable elements in pollen. Cell 136(3):461–472. Ma L, et al. 2015. Angiosperms are unique among land plant lineages in Slotkin RK, Freeling M, Lisch D. 2005. Heritable transposon silencing initi- the occurrence of key genes in the RNA-directed DNA methylation ated by a naturally occurring transposon inverted duplication. Nat (RdDM) pathway. Genome Biol Evol. 7(9):2648–2662. Genet. 37(6):641–644. Malthus T. 1798. An essay on the principle of population (ed. 1). London: Szitenberg A, et al. 2016. Genetic drift, not life history or RNAi, determine J. Johnson in St Paul’s Church-yard. long-term evolution of transposable elements. Genome Biol Evol. Mari-Ordonez A, et al. 2013. Reconstructing de novo silencing of an active 8(9):2964–2978. plant retrotransposon. Nat Genet. 45(9):1029–1039. Teixeira FK, Heredia F, Sarazin A, et al. 2009. A role for RNAi in the selec- Martinez G, Choudury SG, Slotkin RK. 2017. tRNA-derived small RNAs tive correction of DNA methylation defects. Science target transposable element transcripts. Nucleic Acids Res. 323(5921):1600–1604. 45(9):5142–5152. Tenaillon MI, Hollister JD, Gaut BS. 2010. A triptych of the evolution of Martinez G, Ko ¨ hler C. 2017. Role of small RNAs in epigenetic reprogram- plant transposable elements. Trends Plant Sci. 15(8):471–478. ming during plant sexual reproduction. Curr Opin Plant Biol. Tian Z, et al. 2009. Do genetic recombination and gene density shape the 36:22–28. pattern of DNA elimination in rice LTR-retrotransposons? Genome Mart ınez G, Panda K, Ko ¨ hler C, Slotkin RK. 2016. Silencing in sperm cells is Res. 19(12):2221–2230. directed by RNA movement from the surrounding nurse cell. Nat Tsuzuki M, et al. 2016. Profiling and characterization of small RNAs in the Plants 2:16030. liverwort, Marchantia polymorpha, belonging to the first diverged land Matzke MA, Kanno T, Matzke AJ. 2015. RNA-directed DNA methylation: plants. Plant Cell Physiol. 57(2):359–372. the evolution of a complex epigenetic pathway in flowering plants. Volterra V. 1926. Pages 409-448 in Chapman RN 1931. New York: Annu Rev Plant Biol. 66:243–267. McGraw-Hill. Matzke MA, Scheid OM, Matzke AJ. 1999. Rapid structural and epigenetic Wicker T, et al. 2004. A detailed look at 7 million years of genome evo- changes in polyploid and aneuploid genomes. Bioessays 21(9):761–767. lution in a 439 kb contiguous sequence at the barley Hv-eIF4E locus: McCue AD, et al. 2015. ARGONAUTE 6 bridges transposable element recombination, rearrangements and repeats. Plant J. 41(2):184–194. mRNA-derived siRNAs to the establishment of DNA methylation. Ye R, et al. 2012. Cytoplasmic assembly and selective nuclear import of EMBO J. 34(1):20–35. Arabidopsis Argonaute4/siRNA complexes. Mol Cell 46(6):859–870. Nuthikattu S, et al. 2013. The initiation of epigenetic silencing of active Zemach A, et al. 2013. The Arabidopsis nucleosome remodeler DDM1 transposable elements is triggered by RDR6 and 21-22 nucleotide allows DNA methyltransferases to access H1-containing heterochro- small interfering RNAs. Plant Physiol. 162(1):116–131. matin. Cell 153(1):193–205. Panda K, Slotkin RK. 2013. Proposed mechanism for the initiation of trans- Zhang H, Xia R, Meyers BC, Walbot V. 2015. Evolution, functions, and posable element silencing by the RDR6-directed DNA methylation mysteries of plant ARGONAUTE proteins. Curr Opin Plant Biol. pathway. Plant Signal Behav. 8:8 e25206. 27:84–90. Perelson AS. 2002. Modelling viral and immune system dynamics. Nat Rev Immunol. 2(1):28–36. Piegu B, Guyot R, Picault N, et al. 2006. Doubling genome size without polyploidization: dynamics of retrotransposition-driven genomic Associate editor: Emmanuelle Lerat Genome Biol. Evol. 10(3):803–815 doi:10.1093/gbe/evy043 Advance Access publication February 21, 2018 815 Downloaded from https://academic.oup.com/gbe/article-abstract/10/3/803/4895092 by Ed 'DeepDyve' Gillespie user on 16 March 2018 http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Genome Biology and Evolution Oxford University Press

Modeling Interactions between Transposable Elements and the Plant Epigenetic Response: A Surprising Reliance on Element Retention

Free
13 pages

Loading next page...
 
/lp/ou_press/modeling-interactions-between-transposable-elements-and-the-plant-U30x0Zddck
Publisher
Oxford University Press
Copyright
© The Author(s) 2018. Published by Oxford University Press on behalf of the Society for Molecular Biology and Evolution.
ISSN
1759-6653
eISSN
1759-6653
D.O.I.
10.1093/gbe/evy043
Publisher site
See Article on Publisher Site

Abstract

Transposable elements (TEs) compose the majority of angiosperm DNA. Plants counteract TE activity by silencing them epigenetically. One form of epigenetic silencing requires 21–22 nt small interfering RNAs that act to degrade TE mRNA and may also trigger DNA methylation. DNA methylation is reinforced by a second mechanism, the RNA-dependent DNA methylation (RdDM) pathway. RdDM relies on 24 nt small interfering RNAs and ultimately establishes TEs in a quiescent state. These host factors interact at a systems level, but there have been no system level analyses of their interactions. Here, we define a deterministic model that represents the propagation of active TEs, aspects of the host response and the accumulation of silenced TEs. We describe general properties of the model and also fit it to biological data in order to explore two questions. The first is why two overlapping pathways are maintained, given that both are likely energetically expensive. Under our model, RdDM silenced TEs effectively even when the initiation of silencing was weak. This rela- tionship implies that only a small amount of RNAi is needed to initiate TE silencing, but reinforcement by RdDM is necessary to efficiently counter TE propagation. Second, we investigated the reliance of the host response on rates of TE deletion. The model predicted that low levels of deletion lead to few active TEs, suggesting that silencing is most efficient when methylated TEs are retained in the genome, thereby providing one explanation for the large size of plant genomes. Key words: methylation, silencing, transposable elements, genome size. Introduction Despite the obvious evolutionary success of TEs, the plant Angiosperm genomes vary >1,000-fold in size, and this var- host checks their proliferation. The two entities engage in a iation correlates strongly with transposable element (TE) con- continuous arms-race, where TEs seek to proliferate and the tent. For plant species with small genomes, like Arabidopsis host attempts to control them (Lisch and Slotkin 2011). In thaliana or Brachypodium distachyon, DNA derived from TEs fact, most—but not all (Li et al. 2010)—TEs are epigenetically constitute 20–30% of the genome (AGI 2000; IBI 2010). silenced under normal conditions (Lisch 2009). The plant host Species with larger genomes have commensurately larger exerts this control by suppressing TE activity both before and proportions of TE-derived DNA. For example, TE-derived after transcription. Posttranscriptional modification relies DNA represents >85% of the barley (Hordeum vulgare)and chiefly on RNAi that recognizes and degrades TE mRNA pro- maize (Zea mays ssp. mays) genomes (Wicker et al. 2004; duced by RNA polymerase II (Pol II). Degradation requires as- Schnable et al. 2009). When one considers that the average sociated factors like RNA-polymerase 6 (RDR6), which size of a diploid angiosperm genome is similar to that of barley converts single-stranded to double-stranded RNA (dsRNA); genome, at 6400 Mb, then it is clear that most extant plant the Dicer-like proteins DCL2 and DCL4 that cleave dsRNAs DNA is derived from TEs (Tenaillon et al. 2010). to produce 21 and 22 nucleotide (nt) small interfering RNAs The Author(s) 2018. Published by Oxford University Press on behalf of the Society for Molecular Biology and Evolution. This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/), which permits non- commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact journals.permissions@oup.com Genome Biol. Evol. 10(3):803–815. doi:10.1093/gbe/evy043 Advance Access publication February 21, 2018 803 Downloaded from https://academic.oup.com/gbe/article-abstract/10/3/803/4895092 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Roessler et al. GBE (siRNAs); and the Argonaute1 (AGO1) protein that guides the RdDM pathway, which preferentially acts on TEs that have siRNAs to mRNAs for cleavage (Fultz et al. 2015). already been targeted for silencing. An important feature of Presumably, 21–22 nt siRNAs can prime multiple cycles of these 24 nt siRNAs is that they can act in trans to guide the mRNA cleavage, but they may have another important func- methylation of TEs that have similar sequence characteristics tion, which is to initiate transcriptional silencing (Nuthikattu to the original TE template (Slotkin et al. 2005; Teixeira et al. et al. 2013; McCue et al. 2015). Hence, 21–22 nt siRNAs can 2009; Ito et al. 2011; Ye et al. 2012; Fultz et al. 2015). Under be seen as dual-purpose, because they are involved in post- this process, 24 nt siRNAs may constitute a kind of “immune transcriptional silencing and also because they initiate DNA memory” that act as a buffer against the possibility of TE methylation (Cuerda-Gil and Slotkin 2016). activity (Fultz et al. 2015). If true, this implies that the strength Transcriptional silencing is achieved through epigenetic of the host epigenetic response is related to the number of modifications like DNA methylation, histone modifications, similar TEs in the genome that have already been silenced. and shifts in nucleosome positioning (Bernatavichute et al. Yet, no studies have explored the potential codependence 2008; Chodavarapu et al. 2010). The first of these, DNA between TE copy numbers and the strength of the host methylation, relies on the RNA-directed DNA methylation response. (RdDM) pathway. RdDM begins when the plant-specific Our final systems-level question concerns a separate pro- RNA polymerase Pol IV transcribes a TE. The resulting single- cess that occurs in cells associated with (but not part of) the stranded RNA is processed into 24 nt siRNAs by RDR2 and germline. In cells such as the pollen vegetative nucleus (Slotkin DCL3, two homologs that are distinct from those employed et al. 2009), some TEs are actively demethylated, expressed, in RNAi. Ultimately, the 24 nt siRNAs guide protein complexes and utilized to produce 21–22 nt siRNAs. These siRNAs are to homologous DNA sequences that are then targeted for then transported to the germline, where they presumably cytosine methylation. Once DNA methylation is established, contribute to stable TE silencing across generations (Slotkin at least two mechanisms act to maintain it. The first is a pos- et al. 2009; Ibarra et al. 2012; Martınez et al. 2016; Martinez itive feedback loop: Pol IV and Pol V, the RNA polymerases and Ko ¨ hler 2017). But what is the systems-level benefit of this involved in RdDM, preferentially act on methylated DNA (Law additional step in the host response, given that there are al- et al. 2013; Johnson et al. 2014), thereby reinforcing silencing ready at least two overlapping pathways dedicated to silenc- (Panda and Slotkin 2013). The second is the maintenance of ing TEs and also that symmetric DNA methylation is typically symmetric CGand CHG(where H¼ A, C, or T) methylation inherited faithfully? during DNA replication and cell division (Law and Jacobsen Here, we address these questions by building a model 2010). Although the switch from RdDM to maintenance is of host: TE interactions based on ordinary differential not well understood (Panda and Slotkin 2013), onceaTEis equations (ODEs). ODE models have been used widely to targeted for DNA methylation the host genome employs study biological phenomena that range from population feedbacks to ensure that the TE reaches and maintains a qui- growth (Malthus 1798), to predator–prey interactions escent state. (Volterra 1926), to the dynamics of viral infection and re- Numerous molecular studies have characterized the RNAi production (Perelson 2002). ODE models have also studied and RdDM pathways (reviewed in Law and Jacobsen 2010; the interactions between TEs and the host response Fultz et al. 2015; Matzke et al. 2015). These have been com- (Abrusan and Krambeck 2006), but without a focus on plemented by evolutionary studies showing that small RNAs plants and with few details of host response mechanisms. are used for TE defense across both prokaryotes and eukar- Our model includes proxies for RNAi, RdDM, and addi- yotes (Blumenstiel 2011)and that most RNAi and RdDM com- tional factors like TE propagation and TE deletion. We ponents are present in early land plant lineages (Huang et al. study properties of the model but also estimate reasonable 2015; Ma et al. 2015; Zhang et al. 2015; Tsuzuki et al. 2016). biological parameters by fitting the model to biological However, several important questions remain about systems- data, specifically from the study of the accumulation of level interactions between TEs and their plant hosts. One ma- the Evade element in an A. thaliana inbred line (Mari- jor question is why the host relies on two mechanisms—that Ordonez et al. 2013). Given these parameter estimates, is, RNAi and RdDM—to silence TEs. Presumably both path- we explore dynamics of the model and address systems- ways are capable of silencing; they are thus overlapping and level questions about host: TE interactions. We focus on potentially redundant. Both require the production of myriad three sets of questions: 1) Are both pre- and posttranscrip- polymerases, methylases and small RNAs and therefore must tional silencing necessary to control TEs? If not, what ad- have some energetic cost (Bousios and Gaut 2016). Why, vantage is gained by having two mechanisms? 2) Given then, are two pathways maintained? One working hypothesis that methylated TEs may be an important source of im- is that they act synergistically, but this hypothesis has yet to be mune memory, does TE deletion affect the dynamics of explored. the host response? And, finally, 3) What is the added ben- A second major question concerns 24 nt siRNAs. As men- efit of a third mechanism for generating 21–22 nt siRNAs tioned earlier, 24 nt siRNAs are predominantly produced by in the germline? 804 Genome Biol. Evol. 10(3):803–815 doi:10.1093/gbe/evy043 Advance Access publication February 21, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/3/803/4895092 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Modeling Interactions between Transposable Elements and the Plant Epigenetic Response GBE X X 2 2 Materials and Methods sqEr ¼ ðE  O Þ þ w  ðE  O Þ : CN CN Exp Exp Equilibria and Stability In this formula, E and O are the expected and observed CN CN Given our ODE model, its equations and its parameters copy number, respectively. The expected copy number was (see Results), we found equilibria by solving for active TEs defined as the sum of aTEs and sTEs obtained from the model. (aTEs), silenced TEs (sTEs) and siRNAs when all equations E and O are the expected and observed values, respec- Exp Exp were equal to zero. The first, trivial equilibrium point was tively, for relative expression. aTE ¼ sTE ¼ siRNA ¼ 0. To derive the stability of this eq eq eq The expected relative expression for generation n was equilbrium, we calculated the Jacobian matrix for the obtained from the model by taking the total expression in ODEs around that equilibrium, which provided: generation 8, which is equal to v multiplied by the aTE copy 2 3 vp  d 00 number at generation 8, and comparing that to the total ex- 6 7 pression at generation n, which is equal to v multiplied by the d 0 6 7 JðÞ 0; 0; 0¼ : 6 7 TE number of aTEs in generation n. Note, however, that our 4 5 measure of relative expression may not correspond perfectly ev 0 1 to that from Mari-Ordonez et al. (2013), because the empir- ical data on relative expression actually compares two genes The resulting characteristic equation is: (Evade and ACT2) within each generation and also because qRT-PCR can be inaccurate, especially when it is used as a detðÞ JðÞ 0; 0; 0 k  I¼ 0 TE 3 3 2 ratio (of ACT2 vs. Evade expression). In the square error (sqER) ¼ k þ vpk þ k vpd  vp þ d equation, we assigned w a weight of 40 to reflect the mag- vpd  d ; nitude of difference in the empirical data, because copy num- (1) ber reached 40 and relative expression plateaued at 1 (fig. 2A). where the solutions of this equation are the eigenvalues. We used a Monte Carlo approach to estimate fitted The equation clearly communicates that stability depends parameters. In this approach, all seven parameters were ini- on a complex relationship among v, p,and d but only on tialized with randomly drawn values from a uniform distribu- these parameters. The critical eigenvalue (i.e., the one tion between 0 and 1, except for v, which was ranged that crosses zero) is in fact (pv – d), and hence the sta- between 0 and 20. We also imposed the constraint that bility of this equilibrium is controlled by this compound pþ e  1.0. Given initial parameters, the sqEr was calculated parameter. See supplementary text, Supplementary as above. A single parameter was then altered, with a step Material online, for additional details based on a rescaled size between 0.1 and 0.1 for all parameters (except v where model. step size was between 1.0 and 1.0). The sqEr was calculated The second equilibrium point is shown in equations (3) and and the iteration moved forward only if sqEr > sqEr ;oth- n nþ1 (4) (see Results) for aTE and sTE ; the corresponding equa- eq eq erwise a new step size would be calculated. All the parame- tion for siRNA is: eq ters (in the following order: v, l, l,p,i,r, e,and d)were iterated through 100 times with 50 steps for each parameter, siRNA ¼ (2) eq r i ed until the final fitted parameters were found with the smallest ð pÞ: sqEr for each run. The initialization and iteration of all param- eters was performed >10,000 times; the lowest sqEr across all We also examined the Jacobian matrix and eigenvalues 10,000 runs was used to define the fitted parameters. We to study stability for this equilibrium point. Because the note, however, that other fitted data sets with low sqEr values stability equation was complex, we analyzed the rescaled produced similar model dynamics (supplementary fig. S6, model for further insights into the stability of this nontrivial Supplementary Material online). equilibrium (see supplementary text, Supplementary Material online). Running the ODE Model Fitted Parameters The ODE model was run using odeint from the scipy.integrate We obtained the data from Mari-Ordonez et al. (2013) by package (https://docs.scipy.org/doc/scipy/reference/integrate. loading their figure 3a onto WebPlotDigitizer (https://autome- html, last accessed February 28, 2018) and python (v2.6.6). ris.io/WebPlotDigitizer/, last accessed February 28, 2018). To Figure 1A was made with draw.oi (https://www.draw.io/,last estimate model parameters that fit the empirical data, we accessed February 28, 2018); all other figures were made with used the sum of least squares method, based on the following R (v. 3.3.2). The heatmaps were made with heatmap2, from formula: the gplots library in R. Genome Biol. Evol. 10(3):803–815 doi:10.1093/gbe/evy043 Advance Access publication February 21, 2018 805 Downloaded from https://academic.oup.com/gbe/article-abstract/10/3/803/4895092 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Roessler et al. GBE FIG.1.—A schematic of the model, with details provided in the text. The dashed arrow represents a step specific to cells that contribute to germline material. Table 1 24 nt siRNA is proportional to the number of silenced TEs Summary of Parameters and Their Fitted Estimates (sTEs). Furthermore, 24 nt siRNAs are considered to be trans-acting and thus may affect numerous TE insertions, in- Parameter Description Fitted cluding active elements. Overall, active TEs (aTEs) may be- Estimate come sTEs through 21–22 nt siRNAs, 24 nt siRNAs, or by a V Amount of Pol II mRNA expressed by active TEs 1.630 combination of both (fig. 1). P Proportion of mRNA that contributes to 0.340 transposition The model includes two additional parameters. The first is Proportion of mRNA that contributes to 0.051 TE deletion from the genome, which occurs at rate d for both 21–22nt siRNA production aTEs and sTEs. The second is the potential for the loss of I The rate at which 21–22nt siRNA initiate 0.062 silencing from TEs over time (e.g., through the loss of meth- methylation ylation), which we assume can lead to reactivation of TEs at R The rate at which 24nt siRNA reinforce 0.025 rate u.When u¼ 0, maintenance of silencing is perfect, but methylation silencing is not maintained when u¼ 1. D The rate of TE deletion per generation 0.161 The model is represented diagrammatically in figure 1 and U The rate of methylation loss per generation 0.000 consists of three differential equations: D The rate of degradation of 21–22nt siRNAs per 0.999 generation dðaTEÞ ¼ðÞ v  p  d  i  siRNA  r  sTE aTE þ u  sTE; dt dðsTEÞ Results ¼ðÞ i  siRNA þ r  sTE aTE ðÞ d þ u sTE; dt A Model of TE Propagation and Silencing Our model assumes that an active TE (aTE) begins as single dðsiRNAÞ ¼ e  v  aTE  d  siRNA: copy and expresses mRNA at rate v (fig. 1 and table 1). dt Among the produced mRNA, a proportion p is transposed into new genomic copies of the TE per host generation. The first equation describes the change in the number of Another proportion, e, of the TE mRNA is processed into aTEs over time; the second describes the change in the 21–22 nt siRNAs. Note that pþ e  1.0 under our model. number of sTEs over time, and the third monitors numbers We assume that the 21–22 nt siRNAs degrade at rate d and of 21–22 nt siRNAs over time. Although these three equa- initiate TE silencing at rate i. Initiation encompasses both post- tions represent our basic model, figure 1 includes a transcriptional silencing (RNAi) and the onset of methylation, dashed arrow representing a fourth process, the epige- following previous models (Nuthikattu et al. 2013; McCue netic remodeling of TEs in the germline. This process will et al. 2015). Finally, 24 nt siRNAs reinforce methylation at be incorporated after we first explore the dynamics of the rate r, representing RdDM. In our model, the amount of basic model. 806 Genome Biol. Evol. 10(3):803–815 doi:10.1093/gbe/evy043 Advance Access publication February 21, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/3/803/4895092 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Modeling Interactions between Transposable Elements and the Plant Epigenetic Response GBE vp We recognize that the model, as presented, is over- sTE ¼ (4) eq parameterized in a mathematical sense. For example, r ie pÞ the parameter v could be eliminated by redefining e and p. We retain the parameter definitions provided with siRNA given by equation (2) (see Materials and eq above throughout the main text, because we believe Methods). These two equations illustrate that aTEs and sTEs their meaning to be biologically intuitive. However, in have similar parameter dependencies. However, equilibrium the supplementary text, Supplementary Material online, values of sTEs depend more explicitly on v and p in the nu- we also present a mathematical treatment that includes merator than does the equilibrium values of aTEs. This is an parameter reduction, rescaling, and more extensive der- interesting observation because v and p are properties of ivations of the model’s analytical properties. We refer to aTEs; it drives home the point that equilibria copy numbers this supplement throughout the main text, where of sTEs relies intricately on the properties of their active coun- appropriate. terparts. The denominator of the two equations clearly indi- cates that increasing r tends to decrease both aTE and sTE . eq eq Model Equilibria Finally, the equations also hint at a complex relationship be- tween equilibrium copy numbers and d, because the latter Once a TE has invaded a host it has three possible fates: it appears twice in the denominator (and once in the nominator may fail to successfully invade and be lost completely; it may for sTEeq). As d increases, these appearances have opposite establish itself and reach an equilibrium number of copies effects on equilibrium values. over time; or it may expand in copy number unabated. An We studied these equilibria using an equivalent rescaled advantage of ODE models is that we can analytically solve model (see supplementary text, Supplementary Material on- the equilibrium points to understand TE invasion behavior line). Our analytical results provided additional insights and parameter dependence. We analyzed equilibria and the about the behavior of the model and particularly the stabil- stability of those equilibria. For these analyses we assumed ity of the nontrivial equilibrium. For example, the nontrivial u ¼ 0and d ¼ 1 for simplicity, but also because it is biolog- equilibrium is stable when (pv  d) is positive and unstable ically reasonable to assume both that maintenance of the when this is negative. Conversely, the trivial equilibrium is silenced state is strong (u ¼ 0), based on the conservation of stable when (pv  d) is negative, which implies that both symmetric methylation, and that siRNAs degrade rapidly equilibria exchange their stability for (pv  d)¼0. (d ¼ 1). We identified two equilibrium points in our system. The first is when there are no TEs and, hence, no 21–22 nt and Fitting the Model to Biological Data 24 nt siRNAs in the host. That is, the equilibrium points for the active copies (aTE It can be difficult to identify biologically reasonable parameter ), silenced copies (sTE ), and eq eq values for ODE models. To address this concern, we fitted the siRNA (siRNA ) are equal to zero. Stability around this eq model to biological data and then perturbed parameter values point provides information as to whether a TE will suc- cessfully invade the genome or be lost. We investigated separately to explore parameter dependencies and to assess effects on TE copy numbers. We fitted the model to experi- stability (see Materials and Methods; see eq. 1)and found mental data from the study of Mari-Ordonez et al. (2013), that it does not rely on any of the parameters associated who characterized the expression and transposition of a with epigenetic processes—that is, i, e,or r. Instead, sta- single-copy of the Evade retroelement that had become bility relies only on the parameters for TE expression, unmethylated in A. thaliana met1-mutant epigenetic recom- propagation, and deletion (v, p,and d; see also supple- mentary text, Supplementary Material online). Although binant inbred lines. By following two lines to generations 14 and 15, they showed that Evade was highly expressed until equation (1) is complex, the Jacobian matrix (see Materials generation 11 and 7, respectively, after which expression and Methods) suggests the intuitive notion that invasion proceeds when expression and propagation (v  p)out- plummeted precipitously, presumably due to host silencing. The number of Evade copies increased rapidly while its expres- competes deletion (d). sion was high, to a maximum of 40 copies after 11 and 7 Once a TE has established its presence in the host, it may increase in number until the second, nontrivial equilibrium generations. To fit our model to their data, we extracted information point (see Materials and Methods). The equilibrium points about Evade copy numbers and relative expression (see for aTEs and sTEs are given by: Materials and Methods). We focused on one inbred line aTE ¼ (3) eq (met1) from their study, because this was the only line for r ie pÞ; which data were sampled for consecutive generations: in total seven generations (from 8 to 14) since the reactivation of the single Evade element. We fitted the model to the Evade data Genome Biol. Evol. 10(3):803–815 doi:10.1093/gbe/evy043 Advance Access publication February 21, 2018 807 Downloaded from https://academic.oup.com/gbe/article-abstract/10/3/803/4895092 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Roessler et al. GBE TE copy number. If TEs were able to grow unabated, there would be an exponential increase at a rate of 0.554 (¼pv) TEs per generation. However, some transcripts are processed into 21–22 nt siRNAs (e ¼ 0.051) that silence TEs at rate i¼ 0.062. These 21–22 nt siRNAs degrade quickly for each host gener- ation (d ¼ 0.999), and therefore any new 21–22 nt siRNAs are not residual, but must be made from active TEs. Once initia- 0 5 10 15 0 5 10 15 tion of methylation has begun as part of i, reinforcement quickly takes hold at rate r¼ 0.025. Eventually, the number Generations Generations of sTEs increases and the number of aTEs decreases, so that aTEs total expression begins to decline. mTEs sTEs TE max As TEs become silenced, they have two fates under our Copy Number model: they can be deleted from the genome or become TE final active again due to loss of silencing (fig. 1). Since loss of si- lencing was very low (u¼ 410 ) inthe fittedparameter set, the main fate of sTEs is to be deleted (d¼ 0.16). As these quiescent TEs are lost, so is the source of reinforcing 24 nt 0 100 200 300 400 500 siRNAs. When reinforcement becomes unreliable, the host Generations loses epigenetic control, the subset of remaining aTEs propa- gate, and the phased cycle begins again. These cycles dissi- FIG.2.—(A) Model fit to the Evade data for total copy number (left) pate in amplitude until equilibria are reached at 20 total TE and relative expression (right). The empirical data from the Evade study are copies, with more sTEs (14) than aTEs (6) (fig. 2B). It is represented by circles; the whiskers indicate SD. The model results based on the fitted parameters (table 1) are represented by the solid line. important to note that the equilibrium is not necessarily static; (B) Long-term behavior of the model, based on the fitted parameters to it can be reached when equal numbers of TEs are created the Evade data. Arrows show TE and TE , which are defined in the max final versus deleted. text. Copy number refers to the summation of aTEs and sTEs. These phased interactions occur with the fitted parame- ters, but decaying oscillations in copy number also occur reg- ularly with other parameter combinations (see supplementary with a Monte Carlo approach that concurrently considered text and fig. S1, Supplementary Material online). Oscillating TE the total TE copy number (i.e., the combined total of aTEs and numbers are not, however, a necessary outcome of the model sTEs) and TE expression. Our set of fitted parameter values are (see examples below and supplementary text, Supplementary reported in table 1. These parameter values produced a good Material online). fit to the copy number data, and a curve of similar shape to the observed relative expression data over time (fig. 2A). (Note that our measure of expression is only a proxy for expression Examining Initiation (i) and Reinforcement (r) measured experimentally; see Materials and Methods.) We We have shown that the model can have complex, oscillating recognize that we have fitted a complex model to relatively dynamics based on parameters inferred from biological data. simple data and that our fitted parameters may represent one These parameters can be modified independently to explore of many potential reasonably fitting parameter sets (but see the importance of various processes. In this section, we assess Materials and Methods). They nonetheless provide a biologi- the effect of perturbing the system by varying either initiation cally plausible foundation for examining model behavior. (i) or reinforcement (r), or both, while holding the remaining parameters to the values estimated from the Evade data. We first set i¼ 0, and the result was both intuitive and trivial. With Model Behavior under Fitted Parameters i¼ 0 silencing never begins. Hence, the number of aTEs Given the fitted parameters, we explored host: TE dynamics trended upward at an exponential rate, with no resulting over 500 generations, monitoring numbers of aTEs, sTEs, and sTEs (fig. 3A). total TE copy number (¼aTEsþ sTEs) (fig. 2B). With these The effectofsetting r to zero was less straightforward, parameter values, the model produces oscillations of all three because i was > 0 and hence silencing was initiated. entities for 200 generations until it reaches an equilibrium. Without reinforcement, copy numbers no longer oscillated, The oscillations of aTEs and sTEs are somewhat out of phase but instead burst and rapidly reached a maximum for both with one another. We interpret these results as reflecting aTEs and sTEs. These copy numbers remained flat, implying a feedbacks in the epigenetic system. When a TE first invades steady state in which silencing was initiated by 21–22 siRNAs a host, the combination of expression (fitted value v¼ 1.63; and there was sufficient transposition to counteract TE dele- table 1) and propagation (p¼ 0.340) create an initial burst in tion. Under these parameter values, the steady state of sTEs 808 Genome Biol. Evol. 10(3):803–815 doi:10.1093/gbe/evy043 Advance Access publication February 21, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/3/803/4895092 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Number of TEs Copy Number 020 40 60 0 1020304050 Relative Expression 0.0 1.0 2.0 3.0 Modeling Interactions between Transposable Elements and the Plant Epigenetic Response GBE FIG.3.—Model behavior with the fitted values for all parameters but initiation (i) and reinforcement (r). (A) Graphs illustrate the effect of setting initiation and reinforcement parameters to zero for active TEs (left) and methylated TEs (right). In both graphs, the gray dashed lines represent the number of TEsbased on the fitted model parameters to the Evade data (see also fig. 2B). (B) Heat maps showing the TE (left) and TE (right) for the total copy number max final (¼ aTEsþ sTEs) based on varied values of initiation (y-axis) and reinforcement (x-axis), with copy number displayed in each cell. The dashed cell in each heat map represents the fitted values (table 1). was higher than that of aTEs (fig. 3A), as with the equilibria when r was low (e.g., r 0.1), the value of i had notable reached with fitted parameters (fig. 2B). effects on both TE and TE For example, when max final. If initiation by 21–22 nt siRNAs is sufficient to reach a r¼ 0.001, TE varied over two orders of magnitude as a max steady state and to control TEs, then what is the advantage function of i. Similarly, TE differed 33-fold when i ranged final of reinforcement by 24 nt siRNAs? Equations (3) and (4) show from 0.001 to 0.99 (fig. 3B). This relationship implies that that the steady-state TE copy numbers depend critically on the reinforcement can counter TE propagation efficiently, even values of r and i, since the denominators become very small when initiation of silencing is weak. This observation held where r and i are small. To explore their interdependencies, true when also adjusting for TE expression (v) and deletion we varied i and r across their parameter ranges and assessed (d)(supplementary fig. S2, Supplementary Material online). total copy numbers (¼aTEsþ sTEs). To help characterize effects, we focused on two descriptive statistics, TE and max The Effects of TE Deletion (d) TE (see fig. 2B) TE is the highest total TE copy number final . max achieved under a set of model parameters, and TE is the Theoretically, high TE deletion rates should be advantageous final total copy number after 5,000 generations, a point by which for the plant host, because they limit opportunities for trans- total aTE and sTE copy numbers have typically reached a position and consequent deleterious mutations. However, steady state. Our analyses show that when r0.5, any high amounts of TE deletion could have consequences for change in i had little effect on TE and TE , so long as immune memory, because quiescent TEs may be a major max final there was at least some initiation (fig. 3B). In contrast, source of trans-acting 24 nt siRNAs (Teixeira et al. 2009; Ito Genome Biol. Evol. 10(3):803–815 doi:10.1093/gbe/evy043 Advance Access publication February 21, 2018 809 Downloaded from https://academic.oup.com/gbe/article-abstract/10/3/803/4895092 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Roessler et al. GBE 0.001 0.01 0.1 0.5 0.99 0 500 1000 1500 Generations 0.001 0.01 0.1 0.5 0.99 FIG.5.—The effect of varying rates of TE deletion (d) on the final 0 500 1000 1500 number of silenced TEs (sTE), active TE (aTE) and total copy number Generations (aTEþ sTE). These calculations used the fitted values for all parameters but d. Note also that d begins at an arbitrarily low value of 0.001; when FIG.4.—Model behavior with the fitted values for all parameters but d¼ 0, the number of sTEs diverges. TE deletion (d), which is varied from 0.001 to 0.99. et al. 2011; Fultz et al. 2015). Hence, high deletion rates may silencing. This supports our observation, based on equilibrium adversely affect the epigenetic response. To illustrate the ef- equations (eqs. 3 and 4), that deletion plays a complex role in fect of deletion on TE copy numbers, we varied the deletion determining aTE and sTE . eq eq parameter d from 0.001 to 0.99 (fig. 4), while holding the remaining parameters to their fitted values (table 1). The model produced four noteworthy results. First, when Additional Parameters TE deletion was very low (d¼0.001), aTEs burst quickly to high copy number (30). After peaking at a total copy num- We also varied values of expression (v), propagation (p), and ber of 80, all TEs were silenced and the population of sTEs loss of silencing (u), while the remaining parameters were declined slowly over time, reflecting the low rate of deletion held at their fitted values. The parameter v was arbitrarily (fig. 4). Throughout this process, there were no aTEs after the ranged between 0 and 5. The chief effect of this range was initial burst. Second, when d increased (0.01 d < 0.5), the on the amplitude and periodicity of TE oscillations. Higher system generated oscillations in the number of aTE and sTEs. expression levels led to more dramatic copy number oscilla- The amplitude, frequency, and equilibrium values (i.e., TE ) tions (supplementary fig. S3, Supplementary Material online). final varied with d. Note that the running average of sTEs exceeded Importantly, at low parameter values (e.g., v0.5) TEs either that of aTEs for these parameter values (fig. 4). Third, when TE did not invade the genome or were maintained at very low deletion was at intermediate levels (d¼ 0.5), aTEs reached a copy numbers (<5 total TEs) over the long term. Varying steady state, but there were very few sTEs. Finally, when the p produced results similar to varying v (supplementary fig. rate of TE deletion was very high (d¼ 0.99), all TEs were re- S4, Supplementary Material online). Increasing p did, how- moved from the genome. ever, tend to lead to higher TE and average copy num- max To further illustrate these dependencies on d,we plotted bers relative to the parameter values we explored for v TE for sTEs, aTEs, and all TEs as a function of d (fig. 5). (supplementary fig. S3, Supplementary Material online). final Overall, these results convey a somewhat counterintuitive This was presumably because there is a trade off with v; idea: if the goal is to have few aTEs, then it is beneficial either as it increases, so does the production of 21–22 nt siRNAs, to have dramatically high rates of TE deletion (e.g., d¼ 0.99) which then potentially affect RNAi. Propagation (p), on the or to have such low (e.g., 0.01–0.1) deletion rates that a res- other hand, contributes only to the proliferation of more ervoir of sTEs is preserved and contributes to reinforcement of TEs. Note that low levels of propagation (p < 0.25) resulted 810 Genome Biol. Evol. 10(3):803–815 doi:10.1093/gbe/evy043 Advance Access publication February 21, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/3/803/4895092 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Silenced TEs (sTEs) Methylated TEs (mTEs) Active TEs (aTEs) 020 40 60 80 010 20 30 40 Modeling Interactions between Transposable Elements and the Plant Epigenetic Response GBE a level proportional to the number of sTEs that were deme- Pollen Reactivation No Pollen Reactivation thylated in the companion cells. That is, dsiRNA ¼ e  v ðÞ aTE þ sTE d  siRNA: dt This equation is represented by the dotted arrow in figure 1. We evaluated the effects of this additional process on the 0 100 200 300 400 500 system with fitted parameter values. The effects were consis- Generations tent: it decreased TE ,TE and the periodicity of copy max final number oscillations (fig. 6). Thus, this additional process yields FIG.6.—TE reactivation in pollen. The black line is based on the model notable decrements in TE copy numbers. with fitted parameters (no pollen reactivation); the dashed line is using the same parameters but including additional feedback for pollen guard cells (pollen reactivation). Both lines indicate total copy numbers Discussion (¼aTEsþ sTEs). The additional mechanism in pollen guard cells is denoted In this study, we have devised an ODE model to examine the by the dashed arrow in figure 1. systems dynamics of TE propagation within the context of the epigenetic response of a plant host (fig. 1). Although there are in no invasion. Hence, TEs cannot invade if expression or clear limitations to our approach, the model has produced at propagation is low. least four fundamental insights. The first is the prediction of Our model also assumes a process of silencing loss (u), for oscillating copy numbers typified by a burst of TE activity, which the most likely example is methylation loss. followed by silencing, deletion, and then reactivity. Despite Methylation loss is known to be low based on empirical these oscillations, the system often reached equilibrium copy data because symmetric methylation is typically maintained numbers (fig. 2). Second, our model emphasizes the impor- faithfully through cell division (Becker et al. 2011). Indeed, tance of reinforcement by RdDM-like processes, because it our fitted parameter estimate was u ¼ 410 , suggesting buffers potential upstream inefficiencies in the initiation of that a very low amount of sTEs become aTEs due to, for silencing (fig. 3). Third, we show that these outcomes are example, leaky maintenance of symmetrical methylation. linked to the rate of TE deletion. Somewhat nonintuitively, Overall, we found that varying the u parameter had little the model predicts that either low or very high levels of dele- effect on model behavior at parameter values < 0.01 (sup- tion lead to more efficient control of the number of aTEs plementary fig. S5, Supplementary Material online). This (figs. 4 and 5). Finally, we show that demethylation within implies that variation in spontaneous demethylation rates germline cells reinforces host defenses by dampening TE is likely to have few effects on the dynamics of host: TE bursts and lowering steady-state copy numbers (fig. 6). interactions unless u varies by several orders of magnitudes Below, we first discuss the caveats of our ODE model before from our fitted estimate. placing our insights into the context of plant genome struc- ture and evolution. TE Reactivation Dampens TE Oscillations Caveats Finally, we incorporated an interesting biological observa- tion—that is, the fact that TEs are activated in some repro- Every model has limitations, and ours is no exception. One ductive tissues, ostensibly to ensure the transmission of a important consideration is that our biological knowledge of complement of siRNAs to egg and sperm (Slotkin et al. the host response is incomplete. For example, the details of 2009; Ibarra et al. 2012; Martınez et al. 2016; Martinez and the initiation of methylation are not yet clear, because there Ko ¨ hler 2017). TEs are known, for example, to be demethy- are at least two competing (but likely nonexclusive) hypothe- lated and reactivated in the pollen vegetative nucleus, which ses as to how the host transitions from RNAi to the RdDM accompanies the sperm cell, but does not contribute DNA response (Mari-Ordonez et al. 2013; Nuthikattu et al. 2013; to the fertilized zygote. The reactivated TEs are sources McCue et al. 2015). Furthermore, some aspects of the host of 21–22 nt sRNAs that are transported to the sperm and pre- response have not been included in our model, such as recent sumably target silencing of TEs in the zygote (Slotkin et al. discoveries that 18–22 nt tRNA fragments (Martinez et al. 2009). The net effect of this process is to increase the numbers 2017; Schorn et al. 2017) and some miRNAs (Creasey et al. of 21–22 nt siRNAs in germline cells; these 21–22 nt sRNA 2014) may interfere with TE replication and propagation. originate not only from aTEs but also from sTEs (see below). However, these additional host mechanisms fit relatively easily We added this mechanism to our model with an equation in our model, because they would likely affect conversion (e) that increases the number of 21–22 nt siRNAs in the system at andinitiation(i)(fig. 1). In this sense, our model already Genome Biol. Evol. 10(3):803–815 doi:10.1093/gbe/evy043 Advance Access publication February 21, 2018 811 Downloaded from https://academic.oup.com/gbe/article-abstract/10/3/803/4895092 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Copy Number 0 1020 3040 Roessler et al. GBE implicitly accounts for some exciting new findings, but other Supplementary Material online). These results likely reflect our new insights may require model modifications. reliance on data from a study in which silencing occurred Another limitation is that we have studied the invasion of rapidly (Mari-Ordonez et al. 2013), but there is other experi- only one TE family. In reality, plant genomes harbor a multi- mental evidence that host defenses react quickly to silence tude of TE types that may interact with each other and also active TEs within a few host generations (Teixeira et al. 2009; vary with respect to the host response. For example, some but Fultz and Slotkin 2017), perhaps even more quickly than the not all TE families in A. thaliana are recognized by endoge- host response to Evade. nous miRNAs (Creasey et al. 2014), and short, nonautono- It is interesting to note that these experimental studies mous DNA elements are methylated less efficiently than contradict numerous genome-wide analyses, which suggest longer, autonomous elements (Hollister and Gaut 2009), per- that TE families experience massive bursts lasting thousands haps in part due to biases in genomic location (Zemach et al. or even millions of years (Piegu et al. 2006; Schnable et al. 2013). Finally, we have used only one data set to fit the 2009; Bousios et al. 2012; Daron et al. 2014). One likely ex- model, which followed the invasion of the Evade TE for a planation for this incongruence may be the difficulty of re- short period of few host generations (Mari-Ordonez et al. solving the occurrence of multiple rounds of episodic bursts 2013). Therelianceon Evade reflects the fact that very few within the expanded timeframes reported by the genome- studies have monitored the copy number and expression of wide studies. Limited resolution may be due to technical issues TEs within a plant genome over time, particularly beginning related to in silico TE identification, accurate age estimation, from recent invasion or reactivation. In short, we recognize and perhaps even heterogeneous rates of TE sequence loss the limitations of the empirical data, but they nonetheless and decay across the genome (Tian et al. 2009). No matter allow a glimpse into model behavior under relevant parame- the cause, the apparent gaps between experimental and ter values. genome-wide studies deserve further thought and consider- ation. Longer term experimental studies that monitor TE copy numbers over time and under different stress conditions Invasion and Oscillations would certainly be welcome contributions to our empirical How long does it take to silence TEs in vivo? Our understand- understanding of host: TE interactions. ing of the duration and intensity of TE amplification bursts remains limited (Bousios and Gaut 2016). In order to be si- Equilibria lenced, a TE must first invade. Based on our model and anal- yses of the stability of the first equilibrium point (where Another question is whether TEs reach long-term equilibria aTEs¼ sTEs¼ siRNA¼ 0), invasion depends on expression within a genome. In our model, the oscillations often reduce (v), propagation (p), and deletion (d) but not on downstream in intensity over time to reach a steady state (figs. 2 and 4 and properties of the host response, such as conversion of TE supplementary figs. S3–S5, Supplementary Material online). transcripts to 21–22 nt siRNAs (e), initiation (i), and reinforce- In this equilibrium, sTEs are found in higher numbers than ment (r). Put simply, pv needs to outpace d for a TE to suc- aTEs whenever (pv)/d > 2(eqs. 1 and 2 and supplementary cessfully invade the host. We also investigated invasion by text, Supplementary Material online). modifying v and p from the fitted parameter values (table 1); Our ODE-based approach regularly predicts two phases of invasion did not occur when expression or propagation were host: TE dynamics: one shaped by oscillating changes in TE low (v < 0.5, supplementary fig. S3, Supplementary Material numbers, and another characterized by an equilibrium. TE online; p < 0.25, supplementary fig. S3, Supplementary evolution has been modeled extensively with population ge- Material online). netic approaches (Charlesworth and Charlesworth 1983; Assuming a TE invades successfully, it has the potential to Charlesworth et al. 1994; Brookfield 2005; Le Rouzic and increase rapidly in copy number. Under our model, we found Deceliere 2005), and the basic models predict that TEs reach that copy numbers often oscillated before reaching an equi- steady-state copy numbers after the first TE invasion through librium (e.g., fig. 2). Mathematically, the prevalence of oscil- either a transposition-selection or transposition-deletion equi- lations is related to the value of the expression ðv  e  iÞ=r (see librium. In other words, they do not predict oscillations prior to supplementary text, Supplementary Material online). an equilibrium. In contrast, some studies have expanded their Oscillations tend to occur when TE expression is low (v), the models to include TE sequence evolution or competition be- proportion of 21–22 nt siRNAs is low (e), initiation (i) is low or tween TEs, and these often predict oscillations in TE copy reinforcement (r)is high (see supplementary text, numbers (Le Rouzic and Capy 2006; Le Rouzic, Boutin, Supplementary Material online). et al. 2007; Le Rouzic, Dupas, et al. 2007). For example, Le Under many parameter values explored in this work, the Rouzic, Boutin et al. (2007) investigated host–parasite inter- maximum duration of a TE burst lasts for only a few dozen actions between autonomous TEs and their nonautonomous generations before they are temporarily silenced and decrease counterparts, and they found oscillations in copy numbers in copy number (figs. 2 and 4 and supplementary figs. S3–S5, between both entities. Notably, the oscillations continued 812 Genome Biol. Evol. 10(3):803–815 doi:10.1093/gbe/evy043 Advance Access publication February 21, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/3/803/4895092 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Modeling Interactions between Transposable Elements and the Plant Epigenetic Response GBE indefinitely; an equilibrium was rarely reached unless there parameter under the Evade model (figs. 4 and 5). If the were very low mutation rates and few adaptive TE insertions. goal is simply to rid the genome of TEs, the most efficient Le Rouzic, Boutin et al (2007) and also Brookfield (2005) have method is to have a very high d (> 0.5) that removes all aTEs argued that equilibria are reached under conditions that are and sTEs. However, deletions are mediated by ectopic recom- probably unrealistic for in vivo TEs. This is because the param- bination and illegitimate recombination (Devos et al. 2002) eters that affect TE dynamics such as selection, transposition, that may introduce a substantial fitness cost due to the po- and deletion are likely to change at faster rates than the time tential for catastrophic mutations (Langley et al. 1988). required to reach an equilibrium. Our model does not include Assuming that high ectopic recombination carries an unac- autonomous and nonautonomous TEs, nor does it allow per- ceptable fitness cost, our model suggests that the next best turbations in subsequent generations. Yet, the focus on active solution to limit the number of aTEs is to have very low rates and sTEs may mimic some characteristics of host–parasite of TE deletion (d 0.01). relationships and may contribute to our observed oscillating Ourargumentis thatthe retention ofsTEs may benefit the dynamics. We must caution, however, that our model is not host by boosting immune memory. In theory, this immune explicitly evolutionary, because it does not consider fitness or memory provides a defense against the invasion of new TEs population variation. that have sequence homology to existing genomic TEs (Fultz and Slotkin 2017) and also against TEs that have escaped silencing and need to be resilenced. Two interesting features The Importance of Overlapping Mechanisms of acquired immune memory are that it is energetically ex- Why do plants maintain two overlapping and energetically pensive but also maintained under frequent cycles of reinfec- costly pathways (RNAi and RdDM) to silence TEs? Here, i tion (Best and Hoyle 2013). Under the parameter values encompasses posttranscriptional silencing and the initiation explored with our model, the system usually reaches a steady of methylation, and r represents RdDM (fig. 1). Our results state in which the copy number of aTEs is >0. To the extent show that only a small amount of i is needed to begin silenc- that these dynamics reflect reality, a nonzero equilibrium of ing of an unrecognized TE, but r is necessary to counter prop- aTEs defines a system in which reinfection is not merely fre- agation efficiently. For example, the host maintains TE copy quent but constant. This observation may explain one feature numbers at low levels even when i is inefficient (e.g., of the selective pressure to maintain RdDM-like mechanisms, i¼ 0.001), so long as r reinforces silencing by a value of even though it seems as if most TEs within plant genomes are r 0.1 (fig. 3B). RNAi is clearly not as efficient at limiting TE effectively silenced. There is also a conjecture that “zombie” copy numbers when there is no RdDM, yet it is essential for TEs are maintained in the genome in order to produce siRNAs silencing TEs (fig. 3A). Hence, to the extent the model is cor- that boost immune memory and can trigger the trans-silenc- rect, it implies that plants must have RNAi to start the process ing of active relatives (Lisch 2009). Indirect in silico evidence of silencing, but RdDM vastly enhances host control over TEs. for the existence of zombie TEs has been recently uncovered The inclusion of another, apparently overlapping mecha- in maize (Bousios et al. 2016). nism—that is, the active demethylation of TEs in cells that Finally, if low rates of TE deletion are somehow beneficial contribute siRNAs to germline cells—further enhances host to the host response, this process could drive genome size silencing (fig. 6). increases over evolutionary time, because each new TE infec- Our data are consistent with the argument that 24 nt tion or TE reactivation adds copies that are silenced, retained, siRNAs are important for buffering TE activity, even though and not quickly deleted. We also note that this is unlikely to be they seem unnecessary because most methylation is main- a run-away process, because there is evidence for selection on tained independently of RdDM in heterochromatic regions genome size (Diez et al. 2013; Bilinski et al. 2017), especially (Zemach et al. 2013). In fact, it was recently shown that these when genome size gets too large (Knight et al. 2005). heterochromatic regions also produce 24 nt siRNAs, albeit to a Nonethless, our model offers a partial explanation for the smallerextent(Li et al. 2015). These findings are consistent high TE contents and sizes of plant genomes. with the idea that 24 nt siRNAs may act as immune memory (Fultz et al. 2015; Fultz and Slotkin 2017), based on evidence that they may play a key role in suppressing reactivated TEs Future Directions (Teixeira et al. 2009; Ito et al. 2011; Fultz et al. 2015; Fultz and This is the first study to explicitly incorporate features of the Slotkin 2017). plant host response into a quantitative model of host: TE dy- namics. We view this model as a foundation for further exten- The Curious Case of TE Deletion sions that will continue to elucidate important features of If 24-nt siRNAs act as a source of immune memory, then the host: TE interactions. One promising avenue will be to extend retention of sTEs may be a benefit to the host, because they our model to include populations, genetic drift, and fitness may be the template for 24 nt siRNA production. This rela- (Szitenberg et al. 2016), perhaps with a potential for rare tionship is implied by our analyses of the deletion (d) beneficial effects (Le Rouzic, Boutin, et al. 2007). Such an Genome Biol. Evol. 10(3):803–815 doi:10.1093/gbe/evy043 Advance Access publication February 21, 2018 813 Downloaded from https://academic.oup.com/gbe/article-abstract/10/3/803/4895092 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Roessler et al. GBE Charlesworth B, Charlesworth D. 1983. The population dynamics of trans- approach is likely to yield more realistic understandings of the posable elements. Genet Res. 42(01):1–27. evolution of host: TE interactions than are available at present. Charlesworth B, Sniegowski P, Stephan W. 1994. The evolutionary dy- It will also be illustrative to model multiple TE families, includ- namics of repetitive DNA in eukaryotes. Nature 371(6494):215–220. ing autonomous and nonautonomous elements, different Chodavarapu RK, Feng S, Bernatavichute YV, et al. 2010. Relationship length and classes of elements, and the possibility of extensive between nucleosome positioning and DNA methylation. Nature 466(7304):388–392. siRNA cross-homologies. Finally, an important future goal will Creasey KM, et al. 2014. miRNAs trigger widespread epigenetically acti- be to mimic reality by introducing stresses and perturbations vated siRNAs from transposons in Arabidopsis. Nature into the model. One potential example of a perturbation is 508(7496):411–415. polyploidy, which is thought to lead to epigenetic repattern- Cuerda-Gil D, Slotkin RK. 2016. Non-canonical RNA-directed DNA meth- ing (Matzke et al. 1999), but for which the causes remain a ylation. Nat Plants 2(11):16163. Daron J, Glover N, Pingault L, et al. 2014. Organization and evolution of mystery. transposable elements along the bread wheat chromosome 3B. Genome Biol. 15(12):546. Supplementary Material Devos KM, Brown JK, Bennetzen JL. 2002. Genome size reduction through illegitimate recombination counteracts genome expansion Supplementary data areavailableat Genome Biology and in Arabidopsis. Genome Res. 12(7):1075–1079. Evolution online. Diez CM, et al. 2013. Genome size variation in wild and cultivated maize along altitudinal gradients. New Phytol. 199:264–276. Fultz D, Choudury SG, Slotkin RK. 2015. Silencing of active transposable elements in plants. Curr Opin Plant Biol. 27:67–76. Acknowledgments Fultz D, Slotkin RK. 2017. Exogenous transposable elements circumvent identity-based silencing, permitting the dissection of expression- We are grateful for feedback from D. Seymour and dependent silencing. Plant Cell 29(2):360–376. D. Wodarz. K.R. is supported by the National Institute of Hollister JD, Gaut BS. 2009. Epigenetic silencing of transposable elements: Biomedical Imaging and Bioengineering, National Research a trade-off between reduced transposition and deleterious effects on neighboring gene expression. Genome Res. 19(8):1419–1428. Service Award EB009418. A.B. is supported by The Royal Huang Y, et al. 2015. Ancient origin and recent innovations of RNA po- Society (Award Numbers UF160222 and RGF\R1\180006). lymerase IV and V. Mol Biol Evol. 32(7):1788–1799. This work was also supported by NSF grant DEB-1655808 Ibarra CA, Feng X, Schoft VK, et al. 2012. Active DNA demethylation in to B.S.G. plant companion cells reinforces transposon methylation in gametes. Science 337(6100):1360–1364. IBI, International Brachypodium Initiative. 2010. Genome sequencing and Literature Cited analysis of the model grass Brachypodium distachyon.Nature Abrusan G, Krambeck HJ. 2006. Competition may determine the diversity 463:763–768. Ito H, et al. 2011. An siRNA pathway prevents transgenerational retro- of transposable elements. Theor Popul Biol. 70:364–375. transposition in plants subjected to stress. Nature 472(7341):115–119. AGI, Arabidopsis Genome Initiative. 2000. Analysis of the genome se- quence of the flowering plant Arabidopsis thaliana.Nature Johnson LM, Du J, Hale CJ, et al. 2014. SRA- and SET-domain-containing proteins link RNA polymerase V occupancy to DNA methylation. 408:796–815. Becker C, et al. 2011. Spontaneous epigenetic variation in the Arabidopsis Nature 507(7490):124–128. thaliana methylome. Nature 480(7376):245–249. Knight CA, Molinari NA, Petrov DA. 2005. The large genome con- straint hypothesis: evolution, ecology and phenotype. Ann Bot Bernatavichute YV, Zhang X, Cokus S, Pellegrini M, Jacobsen SE. 2008. Genome-wide association of histone H3 lysine nine methylation with 95:171–190. Langley CH, Montgomery E, Hudson R, Kaplan N, Charlesworth B. 1988. CHG DNA methylation in Arabidopsis thaliana. PLoS One 3(9):e3156. Best A, Hoyle A. 2013. The evolution of costly acquired immune memory. On the role of unequal exchange in the containment of transposable Ecol Evol. 3(7):2223–2232. element copy number. Genet Res. 52(3):223–235. Law JA, et al. 2013. Polymerase IV occupancy at RNA-directed DNA meth- Bilinski P, et al. 2017. Parallel altitudinal clines reveal adaptive evolution of genome size in zea mays. bioRxiv. Available from: https://www.biorxiv. ylation sites requires SHH1. Nature 498(7454):385–389. Law JA, Jacobsen SE. 2010. Establishing, maintaining and modifying DNA org/content/early/2017/07/13/134528, last accessed February 28, 2018. methylation patterns in plants and animals. Nat Rev Genet. Blumenstiel JP. 2011. Evolutionary dynamics of transposable elements in a 11(3):204–220. Le Rouzic A, Boutin TS, Capy P. 2007. Long-term evolution of transposable small RNA world. Trends Genet. 27(1):23–31. Bousios A, et al. 2012. The turbulent life of Sirevirus retrotransposons and elements. Proc Natl Acad Sci U S A. 104(49):19375–19380. Le Rouzic A, Capy P. 2006. Population genetics models of competi- the evolution of the maize genome: more than ten thousand elements tion between transposable element subfamilies. Genetics tell the story. Plant J. 69(3):475–488. Bousios A, et al. 2016. A role for palindromic structures in the cis-region of 174(2):785–793. Le Rouzic A, Deceliere G. 2005. Models of the population genetics of maize Sirevirus LTRs in transposable element evolution and host epi- genetic response. Genome Res. 26(2):226–237. transposable elements. Genet Res. 85(3):171–181. Le Rouzic A, Dupas S, Capy P. 2007. Genome ecosystem and transposable Bousios A, Gaut BS. 2016. Mechanistic and evolutionary questions about elements species. Gene 390(1–2):214–220. epigenetic conflicts between transposable elements and their plant hosts. Curr Opin Plant Biol. 30:123–133. Li H, Freeling M, Lisch D. 2010. Epigenetic reprogramming during vege- tative phase change in maize. Proc Natl Acad Sci U S A. Brookfield JF. 2005. The ecology of the genome – mobile DNA elements and their hosts. Nat Rev Genet. 6(2):128–136. 107(51):22184–22189. 814 Genome Biol. Evol. 10(3):803–815 doi:10.1093/gbe/evy043 Advance Access publication February 21, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/3/803/4895092 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Modeling Interactions between Transposable Elements and the Plant Epigenetic Response GBE Li S, et al. 2015. Detection of Pol IV/RDR2-dependent transcripts at the expansions in Oryza australiensis, a wild relative of rice. Genome genomic scale in Arabidopsis reveals features and regulation of siRNA Res. 16(10):1262–1269. biogenesis. Genome Res. 25(2):235–245. Schnable PS, Ware D, Fulton RS, et al. 2009. The B73 maize genome: Lisch D. 2009. Epigenetic regulation of transposable elements in plants. complexity, diversity, and dynamics. Science 326(5956):1112–1115. Annu Rev Plant Biol. 60:43–66. Schorn AJ, Gutbrod MJ, LeBlanc C, Martienssen R. 2017. LTR-retrotrans- Lisch D, Slotkin RK. 2011. Strategies for silencing and escape: the ancient poson control by tRNA-derived small RNAs. Cell 170(1):61–71.e11. struggle between transposable elements and their hosts. Int Rev Cell Slotkin RK, et al. 2009. Epigenetic reprogramming and small RNA silencing Mol Biol. 292:119–152. of transposable elements in pollen. Cell 136(3):461–472. Ma L, et al. 2015. Angiosperms are unique among land plant lineages in Slotkin RK, Freeling M, Lisch D. 2005. Heritable transposon silencing initi- the occurrence of key genes in the RNA-directed DNA methylation ated by a naturally occurring transposon inverted duplication. Nat (RdDM) pathway. Genome Biol Evol. 7(9):2648–2662. Genet. 37(6):641–644. Malthus T. 1798. An essay on the principle of population (ed. 1). London: Szitenberg A, et al. 2016. Genetic drift, not life history or RNAi, determine J. Johnson in St Paul’s Church-yard. long-term evolution of transposable elements. Genome Biol Evol. Mari-Ordonez A, et al. 2013. Reconstructing de novo silencing of an active 8(9):2964–2978. plant retrotransposon. Nat Genet. 45(9):1029–1039. Teixeira FK, Heredia F, Sarazin A, et al. 2009. A role for RNAi in the selec- Martinez G, Choudury SG, Slotkin RK. 2017. tRNA-derived small RNAs tive correction of DNA methylation defects. Science target transposable element transcripts. Nucleic Acids Res. 323(5921):1600–1604. 45(9):5142–5152. Tenaillon MI, Hollister JD, Gaut BS. 2010. A triptych of the evolution of Martinez G, Ko ¨ hler C. 2017. Role of small RNAs in epigenetic reprogram- plant transposable elements. Trends Plant Sci. 15(8):471–478. ming during plant sexual reproduction. Curr Opin Plant Biol. Tian Z, et al. 2009. Do genetic recombination and gene density shape the 36:22–28. pattern of DNA elimination in rice LTR-retrotransposons? Genome Mart ınez G, Panda K, Ko ¨ hler C, Slotkin RK. 2016. Silencing in sperm cells is Res. 19(12):2221–2230. directed by RNA movement from the surrounding nurse cell. Nat Tsuzuki M, et al. 2016. Profiling and characterization of small RNAs in the Plants 2:16030. liverwort, Marchantia polymorpha, belonging to the first diverged land Matzke MA, Kanno T, Matzke AJ. 2015. RNA-directed DNA methylation: plants. Plant Cell Physiol. 57(2):359–372. the evolution of a complex epigenetic pathway in flowering plants. Volterra V. 1926. Pages 409-448 in Chapman RN 1931. New York: Annu Rev Plant Biol. 66:243–267. McGraw-Hill. Matzke MA, Scheid OM, Matzke AJ. 1999. Rapid structural and epigenetic Wicker T, et al. 2004. A detailed look at 7 million years of genome evo- changes in polyploid and aneuploid genomes. Bioessays 21(9):761–767. lution in a 439 kb contiguous sequence at the barley Hv-eIF4E locus: McCue AD, et al. 2015. ARGONAUTE 6 bridges transposable element recombination, rearrangements and repeats. Plant J. 41(2):184–194. mRNA-derived siRNAs to the establishment of DNA methylation. Ye R, et al. 2012. Cytoplasmic assembly and selective nuclear import of EMBO J. 34(1):20–35. Arabidopsis Argonaute4/siRNA complexes. Mol Cell 46(6):859–870. Nuthikattu S, et al. 2013. The initiation of epigenetic silencing of active Zemach A, et al. 2013. The Arabidopsis nucleosome remodeler DDM1 transposable elements is triggered by RDR6 and 21-22 nucleotide allows DNA methyltransferases to access H1-containing heterochro- small interfering RNAs. Plant Physiol. 162(1):116–131. matin. Cell 153(1):193–205. Panda K, Slotkin RK. 2013. Proposed mechanism for the initiation of trans- Zhang H, Xia R, Meyers BC, Walbot V. 2015. Evolution, functions, and posable element silencing by the RDR6-directed DNA methylation mysteries of plant ARGONAUTE proteins. Curr Opin Plant Biol. pathway. Plant Signal Behav. 8:8 e25206. 27:84–90. Perelson AS. 2002. Modelling viral and immune system dynamics. Nat Rev Immunol. 2(1):28–36. Piegu B, Guyot R, Picault N, et al. 2006. Doubling genome size without polyploidization: dynamics of retrotransposition-driven genomic Associate editor: Emmanuelle Lerat Genome Biol. Evol. 10(3):803–815 doi:10.1093/gbe/evy043 Advance Access publication February 21, 2018 815 Downloaded from https://academic.oup.com/gbe/article-abstract/10/3/803/4895092 by Ed 'DeepDyve' Gillespie user on 16 March 2018

Journal

Genome Biology and EvolutionOxford University Press

Published: Mar 1, 2018

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 12 million articles from more than
10,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Unlimited reading

Read as many articles as you need. Full articles with original layout, charts and figures. Read online, from anywhere.

Stay up to date

Keep up with your field with Personalized Recommendations and Follow Journals to get automatic updates.

Organize your research

It’s easy to organize your research with our built-in tools.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

Monthly Plan

  • Read unlimited articles
  • Personalized recommendations
  • No expiration
  • Print 20 pages per month
  • 20% off on PDF purchases
  • Organize your research
  • Get updates on your journals and topic searches

$49/month

Start Free Trial

14-day Free Trial

Best Deal — 39% off

Annual Plan

  • All the features of the Professional Plan, but for 39% off!
  • Billed annually
  • No expiration
  • For the normal price of 10 articles elsewhere, you get one full year of unlimited access to articles.

$588

$360/year

billed annually
Start Free Trial

14-day Free Trial