
==== Front
Genetics
Genetics
genetics
Genetics
0016-6731
1943-2631
Oxford University Press US

38733620
10.1093/genetics/iyae080
iyae080
Investigation
Population and Evolutionary Genetics
AcademicSubjects/SCI01180
AcademicSubjects/SCI01140
Deleterious mutation/epimutation–selection balance with and without inbreeding: a population (epi)genetics model
https://orcid.org/0009-0006-8826-8643
Chernomas Gregory Department of Integrative Biology, University of Guelph, Guelph, Ontario N1G 2W1, Canada

Griswold Cortland K Department of Integrative Biology, University of Guelph, Guelph, Ontario N1G 2W1, Canada

Agrawal A Editor
Corresponding author: Department of Integrative Biology, University of Guelph, Guelph, Ontario N1G 2W1, Canada. Email: gchernom@uoguelph.ca
Conflicts of interest The authors declare no conflict of interest.

7 2024
11 5 2024
11 5 2024
227 3 iyae08026 3 2024
02 5 2024
27 5 2024
© The Author(s) 2024. Published by Oxford University Press on behalf of The Genetics Society of America.
2024
https://creativecommons.org/licenses/by/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Epigenetics in the form of DNA methylation and other processes is an established property of genotypes and a focus of empirical research. Yet, there remain fundamental gaps in the evolutionary theory of epigenetics. To support a comprehensive understanding of epigenetics, this paper investigates theoretically the combined effects of deleterious mutation and epimutation with and without inbreeding. Both spontaneous epimutation and paramutation are considered to cover a broader range of epigenetic phenomena. We find that inbreeding generally reduces the amount of segregating deleterious genetic and epigenetic variation at equilibrium, although interestingly inbreeding can also increase the amount of deleterious genetic or epigenetic variation. Furthermore, we also demonstrate that epimutation indirectly can cause increased or decreased deleterious genetic variation at equilibrium relative to classic expectations, which is particularly evident when paramutation is occurring. With the addition of deleterious epimutation, there may be significantly increased purging of deleterious variation in more inbred populations and a significantly increased amount of segregating deleterious variation in more outbred populations, with notable exceptions.

population genetics
epigenetics
epimutation
paramutation
inbreeding
Natural Sciences and Engineering Research Council of Canada 10.13039/501100000038
==== Body
pmcIntroduction

There is a growing literature whereby epigenetic phenomena are incorporated into theoretical evolutionary models. For epigenetic inheritance to be of evolutionary significance, it has to be heritable, variable in populations, and to have nonnegligible effects on fitness (Bonduriansky and Day 2018). Several previous works have presented results pertaining to the adaptive significance of epialleles and conditions that would favor epigenetic variation as a source of adaptive variation (Klironomos et al. 2013; Uller et al. 2015; Kronholm and Collins 2016; Greenspoon and Spencer 2018; Smithson et al. 2019). This sort of work is of great interest; however, a key additional consideration is the deleterious effects of epimutation, and with explicit predictions on the resulting amounts of segregating deleterious epigenetic variation in populations. Some studies have accounted for the possibility that epialleles can be adaptive and/or deleterious (Stenøien and Pedersen 2005; Geoghegan and Spencer 2012, 2013a, 2013b, 2013c; Furrow 2014; Webster 2024); see discussion for details. Yet, the consequences of inbreeding on deleterious epigenetic variation and its interaction with other processes, such as mutation and paramutation, have not been systematically studied, theoretically.

The effects of inbreeding and breeding systems in general may be of special relevance to K. Plantae, as they are a clade more likely to house heritable epigenetic variation (Charlesworth et al. 2017), and where breeding systems are quite variable (Charlesworth 2006). Several epialleles have been identified in plants (Jacobsen and Meyerowitz 1997; Cubas et al. 1999; Soppe et al. 2000; Manning et al. 2006; Martin et al. 2009; Miura et al. 2009; Silveira et al. 2013; Quadrana et al. 2014). These studies indicate a range of observations of the fitness effects of epialleles, including effects that are likely to be deleterious and others which may be a source of deleterious or adaptive fitness variation. Nevertheless, more direct fitness measurements of naturally occurring, heritable, nonneutral epialleles are needed to characterize the distribution of fitness effects of epigenetic variation in natural populations.

Two specific forms of epimutation are spontaneous epimutations, which are analogous to random mutations, and paramutation, which is analogous to gene conversion (Geoghegan and Spencer 2013c). Specifically, spontaneous epimutation refers to a random change in the state of an epigenetic modification, such as the methylation status of a particular nucleotide (van der Graaf et al. 2015). Paramutation is the conversion of the gene expression state of one allele to an identical or similar expression state as its homologous allele through an in trans interaction (Geoghegan and Spencer 2013c). Of note, paramutation can occur at higher rates than rates of mutation and spontaneous epimutation (Springer and McGinnis 2015; Hollick 2017).

Paramutation (or similar phenomena) has been discovered in a diverse set of plant species including maize (Brink 1956; Coe 1966; Pilu et al. 2009; Sidorenko et al. 2024), Arabidopsis (Luff et al. 1999), petunia (van Houwelingen et al. 1999), and tomato (Gouil and Baulcombe 2018). These examples relate to endogenous genes (Sidorenko et al. 2024). There have also been discoveries of paramutation occurring at transgenic loci (Mittelsten Scheid et al. 2003; Khaitová et al. 2011; Xue et al. 2012; Gao and Zhao 2013; Bente et al. 2021). This highlights the potential importance of paramutation in relation to genetically modified crops. In addition, paramutation-like processes have been observed to occur in mice (Rassoulzadegan et al. 2006) and Drosophila (Capovilla et al. 2017; Dorador et al. 2022), suggesting broader evolutionary implications for paramutation across distantly related clades. Other studies also suggest paramutation may be a more common phenomenon than currently realized and a cause of offspring having distinct gene expression patterns than either parent at a given locus (He et al. 2010; Groszmann et al. 2011; Barber et al. 2012; Shen et al. 2012; Li et al. 2013; Springer and McGinnis 2015). This paper focuses on the potential deleterious effects of paramutation through conversion of a wild-type gene expression pattern to an aberrant expression state (i.e. a deleterious epiallele).

With regard to breeding systems, one expects spontaneous epimutations to behave in a similar manner to classic work on deleterious mutation–selection balance (Crow and Kimura 2010), whereby inbreeding should increase the effectiveness of selection and reduce the amount of segregating deleterious epigenetic variation imposed by spontaneous epimutations. Inbreeding is also expected to reduce the amount of deleterious epigenetic variation imposed by the process of paramutation, but due to an additional reason (Springer and McGinnis 2015). Paramutation depends on an (epi)heterozygous pairing of a paramutagenic (epi)allele and a paramutable allele, and therefore, inbreeding (which decreases heterozygosity) is expected to reduce the effective paramutation rate and decrease the amount of deleterious epialleles generated through this process (Springer and McGinnis 2015).

A key aspect of epimutation in terms of distinctiveness relative to mutation is an increased rate of occurrence (van der Graaf et al. 2015) and highly variable reversion rates. With regard to reversion rates, epialleles can be highly unstable, reverting within a few generations (Bell and Hellmann 2019), or can stably segregate for many generations (Vastenhouw et al. 2006; Johannes et al. 2009; Hollick 2017). Both possibilities are considered in this paper. We also consider the degree of susceptibility to epimutation for different genetic variants (Springer and McGinnis 2015; Johannes and Schmitz 2019), such that epimutation may have an asymmetric fitness effect on a given genetic variant vs another.

To our knowledge, this paper is the first to provide a systematic theoretical treatment of deleterious epimutation–selection balance (both spontaneous mutation and epimutation) with the inclusion of inbreeding. As part of this systematic treatment, we incorporate paramutation to broaden the results for different epigenetic processes that may occur in nature. The structure of the paper is such that first a model of deleterious epimutation–selection balance is derived and then analytical results of equilibria states are presented, with complementary numerical analysis.

Overview of genetic and epigenetic states

The model covers a general context where there is a most fit wild-type allele, with, for example, an optimized level of gene expression. In line with classical work (Crow and Kimura 2010), mutation occurs, acting as a deleterious force. In addition, epimutation in the form of both spontaneous epimutation and paramutation acts against selection for the wild-type allele as it modifies the expression state of the wild-type allele, converting it to an epiallelic variant. The epiallele is assumed to have a lower fitness than the wild-type allele, in line with the general expectation that epimutation is typically neutral or deleterious (Charlesworth et al. 2017).

All cases analyzed have two allelic states and one epiallelic state of the wild-type allele. This situation is expected if one allele has a sequence difference that makes it susceptible to epimutation, while the other allele is not epimutable (Springer and McGinnis 2015; Quadrana and Colot 2016; Johannes and Schmitz 2019). Specifically, for example, this may occur if allelic variants differ in the number and/or placement of cytosine residues in their sequences allowing for potentially different methylation patterns. Another possibility is if genetic variants differ as a result of a transposon or repetitive element mutation (Ong-Abdullah et al. 2015; Springer and McGinnis 2015; Quadrana and Colot 2016; Johannes and Schmitz 2019), such that epimutation only occurs on the genetic variant that possesses the transposon/repetitive element. An alternative context for the model is whereby spontaneous epimutation occurs at both alleles, but has a neutral fitness effect on one allele and that same allele is not paramutagenic (Springer and McGinnis 2015). In line with this, it is assumed the deleterious allele is not directly affected by the epimutation process, but the wild-type allele is directly impacted across all cases analyzed.

Model

Consider a population that is infinite in size, diploid, and sexually reproducing. At a locus, there are three (epi)allelic states, including a wild-type allele (A), a deleterious allele (a), and an epiallelic variant of the wild-type allele (B). Table 1 presents all variables and parameters used in the model with symbols indicated.

Table 1. Names and symbols for variables and parameters.

Variable/Parameter	Symbol	
Frequency of allele (A)	pA	
Frequency of allele (a)	pa	
Frequency of epiallele (B)	pB	
Relative fitness of genotype i	wi	
Inbreeding coefficient	f	
Paramutation rate	m	
Forward spontaneous epimutation rate	t 1	
Reverse spontaneous epimutation rate	t 2	
Mutation rate of A ⇒ a	u	
Mutation rate of a ⇒ A	z	
Mutation rate of B ⇒ a	u′	
All variables and parameters are restricted to be between 0 and 1 in magnitude.

Forward spontaneous epimutation of the wild-type allele to the epiallele occurs at rate t1 and the reverse occurs at rate t2.

A→t1B

A←t2B

Paramutation involving only (epi)genotype A/B occurs at rate m (Geoghegan and Spencer 2013c):

A/B→mB/B

The assumed sequence of events is mating, paramutation, natural selection [normalized by mean fitness (w¯)], spontaneous epimutation, and mutation (see Fig. 1). This is consistent with classic work and assumptions whereby selection occurs after mating, and mutation is assumed to occur after selection (Crow and Kimura 2010). Following this order of events, the epiallele homozygote genotype frequency (pB, B) is

pB,B(1)=pB2(1−f)+pBf(followingmating)

pB,B(2)=pB,B(1)+pA,B(1)m(followingparamutation)

pB,B(3)=w4pB,B(2)w¯(postselection)

Mating, the possibility of paramutation, and then selection were modeled for all genotypes, such that

pB(3)=pBB(3)+12pAB(3)+12pBa(3)andpA(3)=pAA(3)+12pAB(3)+12pAa(3).

Fig. 1. Life cycle diagram outlining the sequence of events for the model. Beginning at the top, the initial event is mating/pairing of gametes to generate zygotes. Zygotes then undergo paramutation at rate “m” if the genotype is A/B, but no other (epi)genotypes are affected. Selection then occurs based on assumed genotype fitness, to generate the postselection gamete frequencies. This is followed by forward and reverse spontaneous epimutation of gametes. Mutation of gametes is the final event, which generates the next generation of gamete frequencies prior to the life cycle repeating.

Spontaneous epimutation occurs next, which for the epiallele results in a frequency of

pB(4)=(1−t2)pB(3)+t1pA(3)

Mutation occurs as the final event in the life cycle. We assume mutation has the same effect on the wild-type allele and epiallele. In the case of the epiallele, mutation results in conversion to the “a” state, such that

pB(5)=pB(4)(1−u′)=pB′

The derivation of the full system of recursion equations consisting of the wild-type allele, epiallele, and deleterious allele frequencies, respectively, is provided in Supplementary File S1.

Relative fitnesses for all genotypes are given in Table 2. Due to the constraint pA + pB + pa = 1, the number of dimensions is reduced to 2, and only the deleterious allele and epiallele were directly tracked.

Table 2. Relative fitness symbols for all (epi)genotypes.

Allele/epiallele	A	B	a	
A	w1	w2	w3	
B	w2	w4	w5	
a	w3	w5	w6	

Analytical techniques used to obtain equilibria and perform local stability analysis are presented in the Appendix, and detailed calculations are provided in the Supplementary Material (Supplementary File S2). For analytical tractability and based on reasonable biological assumptions, rates of forward spontaneous epimutation and mutation were assumed to be small in magnitude [i.e. u, u′, z, t1 ∼ O(ζ)] for all cases. This assumption allowed for Taylor's approximations in determining equilibria and stability terms (see Appendices A1 and A2 for details).

Despite the approximations, the eigenvalues were generally still complicated making analytical local stability interpretations difficult. Accordingly, numerical assessments of the approximated eigenvalues and complimentary numerical simulations of the recursion equations were performed. We find only one equilibrium is simultaneously biologically valid and locally stable for a given area of parameter space for cases A–D.

Results

For the results of the main text, incomplete dominance on fitness was assumed (Table 3) for all cases (A–E). Supplementary File S3 examines a complete dominance on fitness scenario where the wild-type allele is completely dominant when present in an (epi)genotype and the epiallele is completely dominant when paired with the deleterious allele. All cases in the main text and Supplementary File S3 allow for two scenarios. One parameter scenario assumes all mutation and spontaneous epimutation rates are small in magnitude [∼ O(ζ)] (see Appendix A1 for details), whereas the second scenario assumes the reverse spontaneous epimutation is not small in magnitude (t2 >> ζ). This accounts for both relatively stable and unstable epialleles, both of which are known to occur (Vastenhouw et al. 2006; Hollick 2017; Bell and Hellmann 2019; Johannes and Schmitz 2019).

Table 3. Assumed fitness of (epi)genotypes for cases A–D.

Allele/epiallele	A	B	a	
A	1	1−h2s2	1−hs	
B	1−h2s2	1−s2	1−h2s2−hs+s3	
a	1−hs	1−h2s2−hs+s3	1 – s	

Cases A and B with no paramutation (m = 0)

Case A: random mating (f = 0)

This case assumes random mating and no paramutation to give baseline predictions of the allele and epiallele frequencies at equilibrium. There is only one biologically valid and stable equilibrium, which to first order is

p^a≈uhs+O[ζ2]

p^B≈t1h2s2+O[ζ2]

(1) p^A≈1−uhs−t1h2s2+O[ζ2]

This equilibrium assumes selection is a stronger force than mutation and spontaneous epimutation causing the wild-type allele to be at a high frequency with the deleterious allele and epiallele at a low frequency at equilibrium. The deleterious allele equilibrium frequency is at the classic expectation of mutation–selection balance with incomplete dominance (Crow and Kimura 2010). The equilibrium frequency of the epiallele has the same form as the deleterious genetic allele. As selection against the epiallele increases or its dominance over the wild-type allele increases, it reaches a lower equilibrium frequency.

Overall, the process of spontaneous epimutation in addition to mutation results in increased segregating deleterious variation at equilibrium, relative to the classic two-allele case (Crow and Kimura 2010). Since, generally t1 > u (van der Graaf et al. 2015), the deleterious epiallele may be at a higher frequency than the deleterious allele, causing an even larger amount of segregating deleterious variation.

Next, we allow the reverse spontaneous epimutation rate to not be small in magnitude (i.e. t2 >> ζ), assuming the forward rate and mutation rates are still relatively small in magnitude (i.e. u, z, u′, t1 ∼ O(ζ)). This corresponds to when the epiallele is unstable. We find one stable equilibrium occurs (eq. 2) that has a similar form to the previous context (eq. 1).

p^a≈uhs+O[ζ2]

p^B≈t1h2s2(1−t2)+t2+O[ζ2]

(2) p^A≈1−uhs−t1h2s2(1−t2)+t2+O[ζ2]

Unlike eq. 1, t2 now affects the wild-type and epiallele frequencies. A higher reverse rate lowers the equilibrium frequency of the epiallele and increases the wild-type frequency compared to eq. 1. As t2 approaches one (complete epiallele resetting), the epiallele equilibrium frequency is completely determined by the forward spontaneous epimutation rate t1 (see Table 4 for a summary of the case A equilibria). Of note, Stenøien and Pedersen (2005) gave explicit analytical results of deleterious epimutation–selection balance with incomplete dominance for diploids, with similar results to eq. 2 (see discussion for details).

Table 4. Summary of case A equilibria [no paramutation (m = 0) and no inbreeding (f = 0)].

Equilibria conditions	p^a	p^B	p^A	Eq. #	
w AA > wAa > waa, wAA > wAB > wBB,	very low	low	high	1	
u << hs,	
t 1 << h2s2	
Epiallele reversion rate is not small (i.e. t2 >> ζ)	
w AA > wAa > waa, wAA > wAB > wBB,	very low	very low–low	high–very high	2	
u << hs,	
t 1 << h2s2, t2	

Separately, we also analyze a distinct fitness scenario where there is complete dominance of the wild-type allele and complete dominance of the epiallele when paired with deleterious allele. A major result of interest is a significantly lower frequency of the deleterious allele at equilibrium compared with the classic two allele expectations with complete dominance. This result likely stems from a higher frequency of the epiallele, due to masking by the wild-type allele, which then allows for more pairings of the epiallele with the deleterious allele, lowering its marginal fitness (see Supplementary File S3, case A for details).

Case B: with inbreeding (f > 0)

Here, we examine the joint effects of incomplete dominance and inbreeding on the equilibrium frequencies of the alleles and epiallele. In line with classic expectations, inbreeding acts to lower p^a and p^B through increased efficiency of selection with a corresponding increase in p^A. The equilibrium approximated to first order is

p^a≈ufs+hs(1−f)+O[ζ2]

p^B≈t1fs2+h2s2(1−f)+O[ζ2]

(3) p^A≈1−ufs+hs(1−f)−t1fs2+h2s2(1−f)+O[ζ2]

Overall, with inbreeding (f > 0), the amount of segregating deleterious genetic and epigenetic variation at equilibrium is reduced in comparison with when random mating occurs (f = 0).

Next, we consider the parameter scenario when the reverse spontaneous epimutation rate is not small (t2 >> ζ), resulting in the equilibrium

p^a≈ufs+hs(1−f)+O[ζ2]

p^B≈t1(fs2+h2s2(1−f))(1−t2)+t2+O[ζ2]

(4) p^A≈1−t1(fs2+h2s2(1−f))(1−t2)+t2−ufs+hs(1−f)+O[ζ2]

Similar to eq. 2, as the reverse spontaneous epimutation rate increases, the effects of incomplete dominance, selection, and, additionally, inbreeding have less of an impact on the equilibrium frequency of the epiallele (see Table 5 for a summary of the case B equilibria).

Table 5. Summary of case B equilibria [no paramutation (m = 0) and with inbreeding (f > 0)].

Equilibria conditions	p^a	p^B	p^A	Eq. #	
w AA > wAa > waa, wAA > wAB > wBB,	very low	low	high	3	
u << hs, fs,	
t 1 << h2s2, fs2	
Epiallele reversion rate is not small (i.e. t2 >> ζ)	
w AA > wAa > waa, wAA > wAB > wBB,	very low	very low–low	high–very high	4	
u << hs, fs,	
t 1 << t2, h2s2, fs2	

Separately, we analyze the context of complete dominance in fitness. Overall, the results for the equilibria are similar as with incomplete dominance, with marginal increases in the epiallele and deleterious allele due to increased masking by the wild-type allele (see Supplementary File S3, case B for details).

Cases C and D with paramutation (m > 0)

Case C: random mating (f = 0)

Paramutation acts as an additional force that can generate epialleles. Importantly, the rate of paramutation can be high compared to genetic mutation and forward spontaneous epimutation (Hollick 2017). Overall, we identify five distinct equilibria when paramutation is incorporated (see Table 6 for a summary). The equilibria occur within mutually exclusive parameter constraints based on biological validity and/or local stability conditions. The distinct parameter conditions for the equilibria are largely a result of the relative fitness differences of the deleterious (epi)genotypes and the relative strength of paramutation compared to selection. Two equilibria are analytically tractable and are presented first; following this, numerical results for the other three equilibria are discussed.

Table 6. Summary of case C equilibria [with paramutation, (m > 0) and no inbreeding (f = 0)].

Equilibria conditions	p^a	p^B	p^A	Eq. #	
w AA > wAa > waa, wAA > wAB > wBB,
u << wAA − wAa,
t1 << wAA − wAB + m(wAB − 2wBB),
m(2wBB − wAB) < (wAA − wAB)
(equilibrium 1)	very low	low	high	5	
w AA > wAa, wAB > wBB > waa, wBa and 2wBB > wAB,
u, u′, z < t1, t2,
(wAB − wBB)/wAB > m > (wAA − wAB)/(2wBB − wAB)
(equilibrium 2)	low	low–high	low–high	AN
(e.g. see Supplementary File S2, fig. C1)	
w AA > wAa, wAB > wBa > waa,  wBB,
u, u′, z < t1, t2,
m  ≳  wAA − wBB, wAA − waa
(equilibrium 3)	moderate–high	moderate–high	low	AN
(e.g. see Fig. 2)	
w AA > wAa > waa,
wAA > wAB > wBB > wBa, wAB(1 − m),
and 2wBB > wAB,
u′ << wBB − wBa,
t2 << wBB − wAB(1 − m),
m > (wAB − wBB)/wAB
(equilibrium 4)	very low	high	low	6	
w AA > wAa > waa > wBa,  wBB,
u, u′, z < t1, t2,
m  ≳  wAA − wBB, wAA − waa
(equilibrium 5)	moderate–high	low–moderate	low	AN (e.g. see Fig. 3)	
AN, analyzed numerically, approximation assumptions: u, u′, z, t1, t2 ∼ O[ζ].

Interestingly, two analytically tractable and biologically relevant equilibria occur roughly depending on whether paramutation is higher or lower than the strength of selection against the epiallele. Specifically, when m(2wBB − wAB) < (wAA − wAB), the following equilibrium arises where the wild-type allele is at a high frequency (approximated to first order):

p^a≈u1−wAa+O[ζ2]

p^B≈t11−wAB+m(wAB−2wBB)+O[ζ2]

(5) p^A≈1−u1−wAa−t11−wAB+m(wAB−2wBB)+O[ζ2]

where wAa=1−hs,wAB=1−h2s2,wBB=1−s2.

The deleterious allele is at classic two-allele expectations of mutation–selection balance with incomplete dominance (Crow and Kimura 2010). Paramutation affects p^B and p^A through the term m(wAB − 2wBB). In particular, if wAB < 2wBB, paramutation increases p^B above the case when m = 0. In contrast, if wAB > 2wBB, then paramutation decreases  p^B.

A second equilibrium arises when paramutation is of sufficient magnitude to cause fixation of the epiallele which is balanced by reverse spontaneous epimutation to maintain the wild-type allele as well as mutation of the epiallele to maintain the deleterious allele. The equilibrium approximated to first order is

p^a≈u′wBBwBB−wBa+O[ζ2]

p^B≈1−t2wBBwBB−wAB(1−m)−u′wBBwBB−wBa+O[ζ2]

(6) p^A≈t2wBBwBB−wAB(1−m)+O[ζ2]

where wBa=1−hs−h2s2+s3,wAB=1−h2s2,wBB=1−s2.

For the equilibrium to be biologically valid, two general constraints must hold: wBB > wBa and wBB. > wAB(1 − m). The first constraint indicates the B/B (epi)genotype is required to have a higher relative fitness than the B/a (epi)genotype. Regarding the second constraint, the model assumes wBB<wAB, and then, paramutation effectively reduces the contribution of the relative fitness of the A/B (epi)genotype such that it can be less than the relative fitness of the epiallele homozygote (B/B).

p^A and p^a are distinct from the classic two-allele case of mutation–selection balance because the epiallele is at the highest frequency. p^A and p^a largely depend on the difference in relative fitness of the epiallele homozygote and the relevant heterozygous pairing with the epiallele.

When the rate of paramutation is sufficiently high, such that the equilibrium corresponding to eq. 5 does not arise biologically [i.e. m(2wBB − wAB) > (wAA − wAB)], but where paramutation is still weaker than the strength of selection against the epiallele (specifically m < (wAB − wBB)/wAB), an equilibrium arises where either the wild-type or epiallele can be at high frequency, or with moderate frequencies of both. Due to analytical complexity, this equilibrium was analyzed numerically. Based on the parameter values examined, the fitness constraint for this equilibrium is wBB > waa, wBa. Paramutation directly increases the epiallele frequency at this equilibrium (see Supplementary File S2, pg. 23, figure C1), but we did not identify any notable parameter effects or interactions for this equilibrium. Based on numerical assessments, a separate equilibrium appears to occur within the fitness constraint: wBa > waa, wBB which implies a compensatory fitness effect when the deleterious allele and epiallele are joined in the same genotype. Of note, the fitness constraints for this equilibrium are mutually exclusive to the previous two equilibria. In terms of biological context, this could arise, for example, if a regulatory mutation and an epimutation both result in a deleterious change in gene expression relative to the wild-type pattern but where one increases gene expression and the other decreases gene expression. This could result in a partial recovery in fitness of the wild-type expression state. Also, in addition, paramutation is near or greater than the strength of selection for this equilibrium.

Figure 2 demonstrates a parameter scenario for this equilibrium with the observed impact on pa and pB as the relative values of wBB and waa vary. It appears that paramutation causes the epiallele to nearly replace the wild-type allele initially (see peaks in Fig. 2) and later maintains the wild-type allele at low frequency. Interestingly, as paramutation maintains the wild-type allele at low frequency, a balance is reached at equilibrium between the deleterious allele and epiallele based on the difference between waa and wBB, which determines whether one will be higher or lower, or if they reach near equal frequencies.

Fig. 2. Plots of numerical simulations of the recursion equations for the deleterious allele frequency (solid line) and epiallele frequency (dotted line) across time. An equilibrium occurs approximately under parameter conditions: m > s, s2 and wB, a > wa, a, wB, B. The fitness relationships presented are: wB, B > wa, a (top panel a, s = 0.02 and s2 = 0.01), wB, B = wa, a (middle panel b, s = s2 = 0.02), and wB, B < wa, a (bottom panel c, s = 0.01 and s2 = 0.02). Other parameters were held constant: s3 = 0.012, t1 = 10−5, t2 = 10−4, u = 10−9, z = 10−9, u′ = 10−9, h = 0.5, h2 = 0.5, m = 0.05.

Another equilibrium arises when waa > wBB > wBa and when paramutation is near or stronger than the strength of selection against the deleterious allele or when waa > wBa > wBB and when paramutation is stronger than selection (for the parameter range examined). Focusing on the second fitness condition, here the a/a genotype is of higher fitness than the B/a and B/B (epi)genotypes. Interestingly, pa and pB cycle as the equilibrium is approached (Fig. 3). Reverse spontaneous epimutation has a clear effect of reducing the deleterious allele equilibrium frequency and increasing the epiallele equilibrium frequency while also reducing the number of generations cycling occurs (for the given parameter values). This is a counterintuitive result, but may be due to the reverse epimutation rate reducing the effective paramutation rate which appears to increase the deleterious allele and decrease the epiallele frequency for this equilibrium (not shown). A possible reason for this effect of paramutation is that it simultaneously reduces the wild-type allele and increases the epiallele frequency in terms of its direct affect but also inputs epialleles into a relatively deleterious (epi)genotype (for this equilibriums' fitness constraints). Selection then causes the epiallele to reach an overall lower equilibrium frequency and a higher frequency of the deleterious allele. Since reverse spontaneous epimutation reduces the effective paramutation rate, this partly mitigates this effect.

Fig. 3. Top panel (a) shows deleterious allele frequency and bottom panel (b) shows epiallele frequency across time as the reverse spontaneous epimutation rate increases going from the gray line (t2 = 10−3), to the black line (t2 = 10−2.7), to the orange line (t2 = 10−2.4) with equilibria indicated by dashed lines. Other parameters were held constant: s = 0.01, s2 = 0.02, s3 = 0.001, t1 = 10−5, u = 10−9, z = 10−9, u′ = 10−9, h = 0.5, h2 = 0.5, m = 0.05. Note, the orange line asymptotes at the lowest frequency in panel a) and highest frequency in panel b).

Next, relaxing the assumption that the reverse spontaneous epimutation rate is small in magnitude (i.e. t2 >> ζ), an equilibrium occurs analogous to that of eq. 5:

p^a≈u1−wAa+O[ζ2]

p^B≈t11−(1−t2)(wAB−m(wAB−2wBB))+O[ζ2]

(7) p^A≈1−u1−wAa−t11−(1−t2)(wAB−m(wAB−2wBB))+O[ζ2]

where wAa=1−hs,wAB=1−h2s2,wBB=1−s2.

The reverse spontaneous epimutation rate directly reduces the impact of paramutation and tends to decrease p^B with a corresponding increase in p^A. As a compliment to the analytical results in eq. 7, we plot the equilibrium frequencies of the epiallele as the paramutation rate varies and for different relationships of wBB and waa (see Fig. 4). When 0 < m < 0.25, an equilibrium arises where the wild type is at a high frequency (eq. 7), and whether wBB or waa is greater does not affect which equilibrium occurs biologically. Once paramutation is more than marginally greater than the epiallele reversion rate (i.e. m > 0.25, when t2 = 0.20, see Fig. 4), a second distinct equilibrium arises where the epiallele ranges from a moderate to high frequency as paramutation increases, at the expense of the wild-type allele. At this equilibrium, wBB and waa do not effect which equilibrium arises until m becomes sufficiently greater than t2. This critical point depends on the magnitudes of the other parameters, but once it occurs, a distinct equilibrium can arise. Once this point is reached, if waa < wBB, the second equilibrium continues to hold, and the epiallele frequency increases with the paramutation rate, with lower frequencies of the wild-type allele and deleterious allele at equilibrium. Alternatively, if waa > wBB, the epiallele remains at a moderate equilibrium frequency and decreases as the paramutation rate increases (as observed in Fig. 4). In contrast, the deleterious allele is also at a moderate frequency at this equilibrium and increases with the paramutation rate (not shown). Overall, these results demonstrate three equilibria arise for case C (when t2  >>  ζ) and the parameter conditions that allow for each.

Fig. 4. Plot of the epiallele equilibrium frequency against the paramutation rate (m) with different relative fitness relationships of the deleterious allele and epiallele. This demonstrates the interaction of parameters that allow for three distinct equilibria to arise biologically. The dashed line indicates the value of the epiallele reversion rate (t2) showing that paramutation must be greater than this rate for two of the equilibria to arise. Other parameters values were: wAA = 1, wAa = 0.995 (gray line), wAa = 0.99 (black line), wAB = 0.99 (gray line), wAB = 0.995 (black line), wBa = 0.985, t1 = 10−5, t2 = 0.2, u = 10−9, z = 10−9, u′ = 10−9.

As a separate analysis, we obtained numerical equilibria results for case C in the context of complete dominance in fitness. A key result is that when the wild-type allele is at a high frequency, the deleterious allele is significantly below classic two-allele expectations with complete dominance. This likely occurs as result of the epiallele reaching a higher frequency, through masking by the wild-type allele, causing more pairings of the epiallele and deleterious allele, reducing the marginal fitness of the deleterious allele (see Supplementary File S3, case C for details and presentation of the other equilibria).

Case D: with inbreeding (f > 0)

This case offers insights into interactions between inbreeding and paramutation on deleterious epimutation–selection balance. Analogous to the previous case with random mating (f = 0), we identify five distinct equilibria (see Table 7 for a summary). The equilibria occur in mutually exclusive parameter conditions, such that more than one biologically valid and locally stable equilibrium is not observed for a given parameter scenario. In addition to the key impacts of relative fitness conditions and the paramutation rate on equilibrium outcomes, the levels of inbreeding also play a major role in determining which equilibrium arises. Two equilibria are analytically tractable and are presented first. Then, numerical results are presented and discussed for the other three equilibria.

Table 7. Summary of case D equilibria [with paramutation (m > 0) and inbreeding (f > 0)].

Equilibria conditions	p^a	p^B	p^A	Eq. #	
w AA > wAa > waa, wAA > wAB > wBB,
u << wAA − fwaa − (1 − f) wAa,
t1 << wAA − fwBB −
(1 − f)(wAB − m(wAB − 2wBB)),
m(1 − f) (2wBB − wAB) <
(wAA − fwBB − (1 − f)wAB)
(equilibrium 1)	very low	low	high	8	
w AA > wAa, wAB > wBB > waa, wBa,
and 2wBB > wAB,
u, u′, z < t1, t2,
(f − wBB − (1 − f) wAB)/wAB >
m (1 − f) >
(wAA − fwBB − (1 − f)wAB)/(2wBB − wAB)
(equilibrium 2)	low	moderate–high	moderate–high	AN
(e.g. see Supplementary File S2, Fig. D1)	
w AA > wAa, wAB > wBa > waa,  wBB,
u, u′, z < t1, t2,
m (1 − f) ≳  wAA − wBB, wAA − waa
(equilibrium 3)	low–high	low–high	low–moderate	AN
(e.g. see Fig. 5)	
w AA > wAa, wAB > wBB > waa, wBa,
and 2wBB > wAB,
u′ << wBB − wBa − f(waa − wBa),
t2 << wBB − f − (1 − f)(1 − m) wAB,
m(1 − f) > (f + (1 − f) wAB − wBB)/wAB
(equilibrium 4)	very low–low	high	low	9	
w AA > wAa, wAB > waa > wBa > wBB,
u, u′, z < t1, t2,
m (1 − f) ≳  wAA − wBB, wAA − waa
(equilibrium 5)	moderate–high	low–moderate	low–moderate	AN (e.g. see Fig. 6)	
AN, analyzed numerically, approximation assumptions: u, u′, z, t1, t2 ∼ O[ζ].

The first equilibrium (eq. 8) arises when the combination of selection and inbreeding is approximately of greater magnitude than the effective paramutation rate.

p^a≈u1−fwaa−(1−f)wAa+O[ζ2]

p^B≈t11−fwB,B−(1−f)(wA,B−m(wA,B−2wB,B))+O[ζ2]

(8) p^A≈1−u1−fwaa−(1−f)wAa−t11−fwB,B−(1−f)(wA,B−m(wA,B−2wB,B))+O[ζ2]

where waa=1−s,wAa=1−hs,wAB=1−h2s2,wBB=1−s2.

The key distinction relative to when f = 0 is the twofold effect of inbreeding, increasing the efficiency of selection through increased homozygosity and reducing the effective paramutation rate. Overall, at this equilibrium, inbreeding decreases the amount of segregating deleterious genetic and epigenetic variation through more effective selection and lowering the effective paramutation rate.

A second equilibrium occurs when the paramutation rate is high enough to outweigh the combined effects of selection and inbreeding, causing the epiallele to reach a high frequency. The equilibrium approximated to the first order is

p^a≈u′wBBwB,B−wB,a+f(wB,a−wa,a)+O[ζ2]

p^B≈1−wBB(t2wBB−f−(1−f)(1−m)wAB+u′wB,B−wB,a+f(wB,a−wa,a))+O[ζ2]

(9) p^A≈t2wBBwBB−f−(1−f)(1−m)wAB+O[ζ2]

where

waa=1−s,wAa=1−hs,wAB=1−h2s2,wBB=1−s2,wBa=1−hs−h2s2+s3.

Two constraints must be met for the equilibrium to arise in a biologically valid form. One constraint that must hold is wBB>wBa+f(waa−wBa). This constraint will always hold if wBB > waa, wBa. Interestingly, when waa > wBa, inbreeding increases  p^a and vice versa when wBa > waa. Therefore, inbreeding can either increase or decrease the amount of deleterious genetic variation at this equilibrium whereby p^a may be higher or lower relative to when f = 0. The second constraint wBB>f+(1−f)(1−m)wAB can never hold if wBB <f. Intuitively, this implies high levels of inbreeding can prevent the epiallele from being at the highest frequency. When this equilibrium (eq. 9) is biologically valid, inbreeding always causes a net reduction in the amount of segregating deleterious epigenetic variation relative to when f = 0.

Based on numerical assessments, another equilibrium occurs with the parameter scenario: (f − wBB − (1 − f) wAB)/wAB > m (1 − f) (2wBB − wAB) > (wAA − fwBB – (1 − f)wAB) and wBB > waa, wBa. The wild-type allele and epiallele were at moderate frequencies and the deleterious allele at low frequency. The key result is that inbreeding reduces the effective paramutation rate which decreases the epiallele frequency at equilibrium and increases the wild-type allele (see Supplementary File S2, pg. 39, figure D1).

A separate equilibrium arises under the mutually exclusive fitness constraints wBa > waa, wBB. As the case with random mating (case C), a prerequisite is that there is a compensatory fitness interaction between the epiallele and the deleterious allele in the same (epi)genotype. With the fitness condition: wBa > waa > wBB, for the parameter values examined, inbreeding increases the deleterious allele equilibrium frequency and decreases that of the epiallele (Fig. 5).

Fig. 5. Plots of numerical simulations of the recursion equations for the deleterious allele frequency (top panel a) and epiallele frequency (bottom panel b) across time. The inbreeding coefficient increases going from the solid line (f = 0), to the dotted line (f = 0.1), to the dashed line (f = 0.2). Other parameters were held constant: s = 0.01, s2 = 0.02, s3 = 0.012, t1 = 10−5, t2 = 10−4, u = 10−9, z = 10−9, u′ = 10−9, h = 0.5, h2 = 0.5, m = 0.05.

Alternatively, when wBa > wBB > waa, inbreeding increases the epiallele frequency at equilibrium and decreases that of the deleterious allele (data not shown). Analogous to the random mating case, paramutation appears to maintain the wild-type allele at a low frequency, and this allows for an equilibrium to unfold based on the relative fitness of the epiallele and deleterious allele homozygotes. As inbreeding increases, the B/a (epi)genotype is reduced with corresponding increases in homozygosity. In turn, the difference in the (epi)homozygote relative fitnesses plays a larger role in determining the equilibrium frequencies. Of note, if inbreeding levels are sufficiently high, this equilibrium no longer occurs, most likely due to a reduction in the effective paramutation rate allowing for selection to cause the wild-type allele to reach a high frequency (i.e. eq. 8).

Analogous to the previous case (without inbreeding, f = 0), cyclical behavior of pa and pB occurs under the fitness condition: waa > wBa > wBB and when paramutation was stronger than selection (m > s, s2). We see the impacts of inbreeding on this dynamic in Fig. 6. Inbreeding lowers the deleterious allele equilibrium frequency and increases that of the epiallele and also increases the length of time that cycling occurs. Importantly, once inbreeding levels are adequately high, the efficiency of selection outweighs the effective paramutation rate and the wild-type allele reaches a high frequency with more marginal amounts of the deleterious variants (as observed in eq. 8).

Fig. 6. Plots of numerical simulations of the recursion equations for the deleterious allele frequency (top panel a) and the epiallele frequency (bottom panel b) across time, demonstrating cycling dynamics. The inbreeding coefficient increases going from the gray line (f = 0), to the black line (f = 0.25), to the orange line (f = 0.5), with equilibria indicated by the dashed lines. Other parameters were held constant: s = 0.01, s2 = 0.02, s3 = 0.001, t1 = 10−5, t2 = 10−3, u = 10−9, z = 10−9, u′ = 10−9, h = 0.5, h2 = 0.5, m = 0.05. Note, the orange line asymptotes at the lowest frequency in panel a) and highest frequency in panel b).

Relaxing the assumption of a small reverse spontaneous epimutation rate (i.e. t2 >> ζ), an equilibrium arises (eq. 10), where the wild-type allele is at a high frequency.

p^a≈u1−fwaa−(1−f)wAa+O[ζ2]

p^B≈t11−(1−t2)(fwBB+(1−f)(wAB−m(wAB−2wBB)))+O[ζ2]

(10) p^A≈1−u1−fwaa−(1−f)wAa−t11−(1−t2)(fwBB+(1−f)(wAB−m(wAB−2wBB)))+O[ζ2]

The key distinction from eq. 8 is the presence of t2 in the denominators of p^B and p^A. Interestingly, as the epiallele reversion rate increases, both inbreeding and paramutation have less of an impact on the equilibrium frequencies of the wild-type allele and epiallele.

As a compliment to the analytical results, we plot the equilibrium frequencies of the epiallele as the paramutation rate varies and for different relationships of wBB and waa (see Fig. 7). Overall, the results are very similar to case C (when t2 >> ζ, Fig. 4). Three equilibria arise, depending on the strength of m relative to t2 and the relative magnitudes of wBB and waa, and the general patterns mirror that of the random mating context (as shown in Fig. 7). The key distinction is that with inbreeding incorporated, there is increased parameter space where the equilibrium arises with the wild-type allele at a high frequency (eq. 10) and a corresponding decrease in the parameter space for the equilibria where the epiallele is at a moderate or high frequency (as demonstrated in Fig. 7). This is a direct result of inbreeding decreasing the effective paramutation rate requiring m to be larger in magnitude for these equilibria to arise.

Fig. 7. Plot of the epiallele equilibrium frequency against the paramutation rate (m) with inbreeding and with different relative fitness relationships of the deleterious allele and epiallele. This figure shows how inbreeding reduces the effective paramutation rate and decreases the parameter space where the epiallele is at an appreciable frequency. This is observed by comparing the f = 0 cases to the f = 0.25 cases. The dashed line indicates the value of the epiallele reversion rate (t2) showing that the effective paramutation rate must be greater than this rate for two of the equilibria to arise. Other parameters were held constant: wAA = 1, wAa = 0.995 (gray line), wAa = 0.99 (black line), wAB = 0.99 (gray line), wAB = 0.995 (black line), wBa = 0.985, t1 = 10−5, t2 = 0.2, u = 10−9, z = 10−9, u′ = 10−9.

Separately, with complete dominance in fitness, the analytical equilibria results generally mirror the incomplete dominance context in terms of equilibria constraints and properties. One key exception is that only four equilibria arise biologically instead of five. Specifically, the equilibrium dependent on the fitness constraint wBa > wBB, waa cannot arise with complete dominance in fitness (see Supplementary File S3, case D, for the analytical and numerical results).

Case E: the deleterious allele is paramutagenic

Following work by Geoghegan and Spencer (2013c), we now allow for a second paramutation rate caused by the presence of the deleterious allele. This property of paramutation has been observed in mice (Rassoulzadegan et al. 2006). Specifically, genotype A/a is now converted to (epi)genotype B/a at rate m2.

A/a→m2B/a

We label the previous paramutation rate m1, and otherwise, the previous assumptions and parameters are kept the same as case D. An equilibrium arises (eq. 11), and overall, this equilibrium has similar properties as eq. 8 (case D) with some key distinctions (see Table 8 for a summary of the case E equilibria).

p^a≈u1−fwa,a−(1−f)(wA,a−m2(wA,a−wB,a))+O[ζ2]

p^B≈t11−fwB,B−(1−f)(wA,B−m1(wA,B−2wB,B))+(1−f)um2wB,a(1−fwa,a−(1−f)(wA,a−m2(wA,a−wB,a)))(1−fwB,B−(1−f)(wA,B−m1(wA,B−2wB,B))))+O[ζ2]

(11) p^A≈1−u1−fwa,a−(1−f)(wA,a−m2(wA,a−wB,a))−t11−fwB,B−(1−f)(wA,B−m1(wA,B−2wB,B))−(1−f)um2wB,a(1−fwa,a−(1−f)(wA,a−m2(wA,a−wB,a)))(1−fwB,B−(1−f)(wA,B−m1(wA,B−2wB,B))))+O[ζ2]

Table 8. Summary of case E equilibria [with paramutation (m1, m2 > 0) and inbreeding (f > 0)].

Equilibria conditions	p^a	p^B	p^A	Eq. #	
w AA > wAa > waa, wAA > wAB > wBB,
u << 1 − fwaa − (1 − f)(wAa − m2(wAa − wBa)),
t1 << 1 − fwBB − (1 − f)(wAB − m1(wAB − 2wBB)),
m1(1 − f) (2wBB − wAB) <
(wAA − fwBB − (1 − f)wAB)
(equilibrium 1)	very low	low	high	11	
w AA > wAa, wAB > wBB > waa, wBa,
u′ << wBB − wBa − f(waa − wBa),
t2 << wBB − f − (1 − f)(1 − m1) wAB,
m1(1 − f) > (f + (1 − f) wAB − wBB)/wAB
(equilibrium 2)	very low–low	high	low	Identical to case D, eq. 9	
w AA > wAa, wAB > waa > wBB,  wBa
z << waa − f − (1 − f)(1 − m2)wAa,
z << waa − (1 − f)wBa − fwBB,
m2 (1 − f) > (f + (1 − f) wAa − waa)/wAa
(equilibrium 3)	very high	very low	very low	12	
AN, analyzed numerically, approximation assumptions: u, u′, z, t1, t2 ∼ O[ζ].

The deleterious allele equilibrium frequency decreases relative to when m2 = 0. This is caused by a reduction in the marginal fitness of the deleterious allele due to the additional paramutation rate converting (epi)genotype A/a to B/a (since by assumption wAa > wBa). Next, focusing on the epiallele equilibrium frequency, the first term is the same result as the equilibrium when m2 = 0 (case D, eq. 8) and therefore has the same interpretation. The second term is new and expected to be positive when the equilibrium is biologically valid and therefore the epiallele equilibrium increases with the additional paramutation rate (m2). The increase is expected to be very small, due to being dependent on the mutation rate (u) in the numerator. However, this term does imply that the epiallele can be maintained at a nonzero frequency even if forward spontaneous is completely absent t1 = 0. This is in contrast to cases A–D for the equilibria where the wild-type is at a high frequency; the epiallele would not be present at equilibrium if t1 = 0.

A second equilibrium occurs where the epiallele is at high frequency (not shown) and is identical to that of case D (eq. 9). The additional source of paramutation does not impact the epiallele (or the alleles) at this equilibrium. This is likely because both the wild-type and the deleterious allele are at a low equilibrium frequency and the A/a genotype (where m2 is relevant) has a negligible impact to the degree of approximation.

Interestingly, a third analytically tractable equilibrium occurs where the deleterious allele is at a high frequency (eq. 12):

p^a≈1−zwa,a(wa,a−(1−f)(1−m2)wB,a−fwB,B)(wa,a−f−(1−f)(1−m2)wA,a)(wa,a−(1−f)wB,a−fwB,B)+O[ζ2]

p^B≈(1−f)zm2wa,awB,a(wa,a−f−(1−f)(1−m2)wA,a)(wa,a−(1−f)wB,a−fwB,B)+O[ζ2]

(12) p^A≈zwa,a(wa,a−f−(1−f)(1−m2)wA,a)+O[ζ2]

It is clear that m2 allows for this equilibrium to occur in a biologically valid form by reducing the fitness contribution of the A/a heterozygote, which is observed through the term (1 − m2)wAa. Without m2, this equilibrium is never biologically valid, since by assumption of the model, wAa > waa always holds. In addition, there are constraints on the fitness of the deleterious allele homozygote. One minimum constraint that must hold is waa > f. Also, the equilibrium is never biologically valid if waa < wBa, wBB, which is intuitive for an equilibrium where the deleterious allele is at a very high frequency.

Interestingly, inbreeding will always increase the epiallele at this equilibrium when wBB > wBa as the marginal fitness of the epiallele increases from less pairings with the deleterious allele and more B/B homozygotes. In contrast, the deleterious allele will always decrease due to inbreeding under the same parameter constraints, whereas inbreeding always increases the wild-type allele at this equilibrium, consistent with classic expectations.

In contrast to previous cases, the constraints on this equilibrium are not mutually exclusive with that of eq. 11 where the wild type is at a high frequency. Therefore, when there is an overlap in parameter space where they are both valid, the equilibrium that arises will depend on the initial frequencies.

As a reminder, case D had three equilibria which were assessed numerically, with interesting properties such as inbreeding causing an increase in the deleterious allele (Fig. 5) and cycling of the (epi)allele frequencies (Fig. 6 and see Table 7 for a summary). As assessed numerically, three equilibria with very similar properties arise for case E (data not shown). Also, when the assumption of a small epiallele reversion rate is relaxed for this case (i.e. t2 >> ζ), we note similar properties at equilibria to previous cases such that the impacts of the other parameters that affect the epiallele are reduced. The net effect is generally a reduction in the epiallele frequency at equilibrium with a corresponding increase in the wild-type allele (data not shown).

Discussion

Overview of results

This paper has presented fundamental theory on deleterious epimutation–selection balance. The equilibrium properties with epimutation can be similar to a three-allele genetic model (discussed below), but unique aspects of epigenetic processes can give rise to conditions when it does not. In particular, when paramutation is occurring and when it is near or stronger than natural selection, unexpected equilibria arise and there are nonintuitive effects of certain parameters. Specifically, we find the deleterious allele or the epiallele can reach high frequencies under the condition of a high paramutation rate and inbreeding can increase deleterious genetic or epigenetic variation. Furthermore, the work gives rise to straightforward directions for empirical research, be it characterizing or accounting for deleterious epimutation or examining uncertainties that limit theoretical and therefore conceptual/general understanding of epigenetics.

Comparison of results to previous population genetics models

A previous paper considered models of deleterious mutation–selection balance with multiple deleterious genetic alleles and incomplete dominance (Clark 1998). Cases A and B (which lack paramutation) could be viewed as a three-allele model with two deleterious alleles, as opposed to one deleterious allele and one deleterious epiallele. However, the distinct biology of epialleles does not allow for some of the assumptions used to obtain analytical results of Clark (1998). Specifically, the assumptions of a single mutation rate of the wild-type allele to distinct deleterious alleles and that the reversion rate is negligible are not generally valid with epialleles. These assumptions were relaxed in Stenøien and Pedersen (2005) which account for the distinct biology of epialleles compared to genetic alleles (see next section for details).

There are similar population genetics contexts to cases C, D, and E (which include paramutation). From a theoretical standpoint, paramutation is analogous to gene conversion (Geoghegan and Spencer 2013c). A previous study looked at the effects of mating systems on gene conversion (Glémin 2010). Inbreeding decreases heterozygous pairings of alleles and epialleles, respectively, and therefore is expected to be of significant importance for both gene conversion and paramutation. These conceptual expectations were confirmed analytically for gene conversion in Glémin (2010) and for paramutation in the current paper. Glémin (2010) highlights how what was previously assumed to be genomic signatures of selection could in fact be gene conversion, such that it represents an additional evolutionary force with its importance underestimated. Paramutation could also be an unaccounted force in shaping epigenetic and genetic variation (Geoghegan and Spencer 2013a).

Although gene conversion and paramutation have similar properties and can have similar outcomes, there are important differences. One finding by Glémin (2010) was an increase in deleterious allele frequencies at equilibrium as a result of gene conversion. We find a similar possibility with the current model. However, in Glémin (2010), gene conversion can cause a direct increase in the deleterious allele, whereas in the case of paramutation the cause is indirect and depends on specific fitness relationships, as a point of distinction. In addition, with paramutation, only one of the alleles may be paramutable and the other is not (Springer and McGinnis 2015). This can put limitations on the ability for paramutation to increase the amount of a deleterious epialleles in comparison with gene conversion increasing the amount of a deleterious allele. Furthermore and interestingly, gene conversion has particularly significant consequences in recombination hotspots, indicating genomic context can be a significant factor (Glémin 2010). In a similar framework, paramutation can depend on the genomic sequence context to occur, but this pattern does not always hold (Springer and McGinnis 2015).

As a final point, our understanding is the deleterious allele with gene conversion is more stable than a deleterious epiallele maintained by paramutation. This distinction arises because of a high epiallele reversion rate. A high epiallele reversion rate can partially or completely mitigate the effects of paramutation (for example see eqs. 7 and 10), something that is unlikely to occur with gene conversion. This suggests the evolution of epigenetic resetting mechanisms can protect against the deleterious potential of paramutation, whereas there may not be a parallel possibility to protecting against the potential accumulation of deleterious alleles due to gene conversion.

Comparison of results to previous theoretical epigenetic models

A previous study derived analytical equilibrium results for deleterious epimutation–selection balance with incomplete dominance in haploids and diploids (Stenøien and Pedersen 2005). Their results for diploids were similar to our case A results in the context where the epiallele reversion rate was not small in magnitude (i.e. t2 >> ζ, see eq. 2). The frequency of the deleterious allele was identical to ours, which was the same in form as classic results of two-allele deleterious mutation–selection balance with incomplete dominance. Also, very similar to our case A (eq. 2) results, the epiallele equilibrium frequency was determined by a balance of selection and the epiallele reversion rate acting against it and the forward epimutation rate maintaining it at a nonzero frequency. In a similar modeling framework, Webster and Phillips (2024) developed a model of mutation/spontaneous epimutation–selection balance, but allowing epialleles to also be adaptive. Interestingly, they found that when an epiallele is adaptive it can compensate for a deleterious allele, which causes an increase in the frequency of the deleterious allele from classic mutation–selection balance. They also highlight that the deleterious allele can reach a high equilibrium frequency with the addition of epigenetic variation, but this occurs in a limited/rare area of a parameter space and requires that the epiallele is adaptive, as well as a very low reverse epimutation rate. In our study, the epiallele is always comparatively deleterious to the wild-type allele, yet the deleterious genetic allele can also reach moderate to very high frequencies through a different mechanism, namely as a result of paramutation. We also find the opposite can occur whereby the deleterious allele is below classic mutation–selection balance (e.g. case E, eq. 11). Generally, we have broadened these modeling frameworks by incorporating inbreeding (cases B, D, E) and paramutation (cases C–E).

There are also theoretical studies that directly relate to cases C–E (with paramutation). Spencer and Geoghegan have made seminal contributions to the field of population epigenetics including a series of papers analyzing environmentally influenced epiallele stability and epigenetic contributions to phenotypic variation alongside genetic variation (Geoghegan and Spencer 2012, 2013a, 2013b). The study most similar to the current paper was Geoghegan and Spencer (2013c). This was the first thorough theoretical treatment of paramutation. Similar to our work, their model has two alleles and an epiallele; however, Geoghegan and Spencer (2013c) do not restrict the fitness assumptions, such that one allele and the epiallele are not always assumed deleterious. They identified parameter conditions of equilibria, including fixation of either allele, a polymorphic equilibrium with a high frequency of the epiallele, as well as a polymorphic equilibrium with appreciable genetic and epigenetic variation. Due to the incorporation of mutation, fixation of an allele was not possible for our model; nevertheless, we observe similar results for the equilibria where one of the variants is at a very high frequency. With regard to polymorphic equilibria, as assessed numerically we also identify parameter conditions where all three variants are at appreciable frequencies. Another finding from Geoghegan and Spencer (2013c) was that increases in the epiallele reversion rate increase the possibility of a unique equilibrium arising. Similarly, we find as the epiallele reversion rate increases to a sufficiently high rate we find only one equilibrium is possible because paramutation becomes negligible.

In contrast to the previous work, our current study does not allow for the possibility that the epiallele is adaptive in addition to being deleterious, and so, we have restricted our analysis to deleterious epimutation–selection balance. However, in our model genetic variation is regeneratable by mutation, as well as epigenetic variation by spontaneous forward epimutation. Furthermore, we studied inbreeding to allow for basic insights into the impacts of breeding systems on deleterious epimutation–selection balance and its effects on the phenomenon of paramutation.

Geoghegan and Spencer (2013c) found the epiallele reversion rate can reduce the number of locally stable equilibria by mitigating the effects of paramutation, which was required for some equilibria to arise. We find inbreeding acts as another factor that can restrict the number of locally stable equilibria, by reducing heterozygosity (which is a prerequisite for paramutation), such that equilibria which depend on a high paramutation rate are unlikely to arise. Our analysis allows for both genetic variation and epigenetic variation to be present at all equilibria due to the addition of mutation and forward spontaneous epimutation rates. Specifically, forward spontaneous epimutation can be an important consideration if epialleles are able to form independent of the (epi)genotype (i.e. regardless of paramutation). However, if spontaneous forward spontaneous epimutation is not occurring (i.e. t1 = 0), we show (in case E, eqs. 11 and 12) the epiallele can be maintained at equilibrium independent of the forward spontaneous epimutation rate. It is maintained by paramutation caused by the deleterious allele and mutation, which highlights an important contribution of incorporating mutation into the model. Also, of note, we allowed for the epiallele to have a distinct mutation rate from the wild-type allele and this can have nontrivial effects on the equilibrium frequency of the deleterious allele, as epialleles can differ by orders of magnitude in their mutation rates (Yi and Goodisman 2021). Following this, the epiallele could have important effects on the amount of segregating deleterious genetic variation even if the epiallele has nearly neutral deleterious effects.

The fitness interaction of the epiallele and deleterious allele is a key biological factor for equilibrium outcomes

Epimutations can be deleterious through changes to gene expression that are random with respect to adaptation. In contrast, mutations can be deleterious either through changes to gene expression or protein sequence that are random with respect to adaptation. Therefore, one possibility is that when both an epimutation and a mutation result in a change in gene expression, there could be a compensatory fitness interaction of the deleterious variants based on how each has modified gene expression. In contrast, this could be less likely when there is a deleterious mutational change in protein sequence combined with a deleterious epimutational change in gene expression. Therefore, the type of mutation (i.e. regulatory vs protein) that occurs could be a relevant factor in nature when epimutation is also occurring.

In line with this discussion, in the parameter scenario where paramutation is occurring and is of sufficient magnitude, we find that the relative fitnesses of the a/a, B/B, and B/a (epi)genotypes are a major determining factor in which equilibrium arises. Three distinct and mutually exclusive (i.e. not simultaneously locally stable) equilibria can occur based on the relative magnitudes of waa, wBa, and wBB. If there is a compensatory fitness interaction of the deleterious allele and the epiallele, this allows for the equilibrium where either variant can be at high frequency or both can be at intermediate equilibrium frequencies based on the relative magnitudes of waa and wBB. Alternatively, if when paired together they are more deleterious (waa, wBB  > wBa), then the deleterious allele can reach high frequencies when waa > wBB and this is also when cycling occurs as the equilibrium is approached. Alternatively, if waa < wBB, this allows for the equilibrium where the epiallele is at a high frequency, essentially taking the place of the wild-type allele in the population (see eqs. 6 and 9). Overall, whether a given mutation or epimutation at a genetic locus is more deleterious and depending on their fitness interaction, this can result in quite different equilibrium outcomes and distinct paths to reaching equilibria.

Epimutation and mating systems

Based on the demonstrated results on the importance of inbreeding effects on deleterious epigenetic variation, there may be patterns in nature reflecting the coevolution of epialleles and mating systems. Comparisons could be performed on sister species (such as Arabidopsis species), for epigenetic modifications such as methylation, and molecular signatures of paramutation. One would expect when comparing highly inbreeding vs randomly mating sister species there would be significantly different epiallele frequencies. Populations with high epimutation rates such as through the process of paramutation may result in near fixation of epialleles in a random mating species, whereas highly inbred sister species would be much more protected from high epimutation rates, reducing the frequency of deleterious epialleles at equilibria. However, we note parameter scenarios where there is a more unexpected effect of the degree of inbreeding. For example, when the epiallele is very deleterious (wBB < 0.5wAB, eg. see eq. 5), paramutation reduces its frequency at equilibrium and this purging effect of paramutation is stronger with random mating. Also, unexpectedly, we find inbreeding can increase the equilibrium frequencies of the deleterious allele or the epiallele under some conditions. This is in direct contrast to the expectations of the increased efficiency of selection against deleterious variation as inbreeding increases within the classic modeling framework. The more unexpected impacts of inbreeding occur when it is at low to moderate levels and when paramutation is stronger than selection. When inbreeding reaches high levels (such that paramutation has a relatively minimal effect), both deleterious genetic variation and epigenetic variation are simultaneously decreased by inbreeding, in line with classic expectations. Overall, these results suggest more nuanced relationships of breeding systems and segregating deleterious variation than anticipated.

General implications of deleterious epimutation–selection balance

Epimutation is expected to have neutral or deleterious effects when impacting phenotypically relevant genetic loci, analogous to classical understanding of mutation (Bonduriansky and Day 2018). As shown, deleterious epimutations may cause significant additional segregating deleterious variation. However, in order for these effects to impact populations, epimutations and the corresponding epialleles must be heritable, variable among individuals, and have nonnegligible effects on fitness (Bonduriansky and Day 2018).

It is possible that as an overall pattern, when epimutation is viewed as a largely deleterious process (when nonneutral, like mutation), but its causal mechanisms are nonetheless invaluable for certain biological roles in many clades, such as transposon methylation and phenotypic plasticity, natural selection minimizes its deleterious consequences through different paths such as generating resetting mechanisms (Charlesworth et al. 2017), increasing reversion rates over forward (van der Graaf et al. 2015), and possibly even reorganizing genomes (Choi and Yuh 2020). Alternatively, for vulnerable populations that possess evolutionarily relevant epigenetic variation, there may be a significant additional risk to their survival, and the co-evolution of epimutation with mating systems will be of prime importance.

Supplementary Material

iyae080_Supplementary_Data

Acknowledgments

The authors thank Dr. Agrawal and the anonymous reviewers for very helpful comments and feedback. The authors also sincerely thank the NSERC for funding this work.

Data availability

The authors affirm that all data necessary for confirming the conclusions of the article are present within the article, figures, tables, and Supplementary material. Supplementary File S1 goes through the detailed analytical derivation of the recursion equations used for the model. Supplementary File S2 contains a detailed calculations and analysis corresponding to the main text. Supplementary File S3 demonstrates the results of the complete dominance context with detailed calculations and analysis with supporting figures.

Supplemental Material

Supplemental material available at GENETICS online.

Funding

This work was funded by the Natural Sciences and Engineering Research Council of Canada Discovery Program.

Appendix

 

A1: determining equilibria

Represented as general functions, the systems of recursion equations for the deleterious allele and epiallele are pa′=f1(pa,pB;n1,1,n1,2,…,n1,x) and pB′=f2(pa,pB;n2,1,n2,2,…,n2,y), such that ni is a general symbol for a parameter. “x” and “y” represent the number of parameters present in the recursion equations for the deleterious allele and epiallele, respectively. The multivariate equilibrium condition is

p^a=f1(p^a,p^B;n1,1,n1,2,…,n1,x)

p^B=f2(p^a,p^B;n2,1,n2,2,…,n2,y)

For each fitness and inbreeding case, a particular set of parameters were assumed to be of order (ζ), such that any parameter assumed to be of order (ζ) for a given case is transformed: ni −> ni~ζ (Otto and Day 2007). Furthermore, the equilibria (p^a, p^B), are assumed to be perturbations from equilibria. Define (pa,0, pB,0) to be zeroth-order equilibria and (pa,1, pB,1) to be first-order perturbations. Accordingly, the following substitutions are made into equations for the equilibrium frequencies: p^a = pa,0 + pa,1ζ +… and pB^ = pB,0 + pB,1ζ +…resulting in f1(ζ)=pa,0 + pa,1ζ + …, f2(ζ)=pB,0 + pB,1ζ +…. f1(ζ) and f2(ζ) are expressions for the recursion equations after the above substitutions for the variables and parameters of order (ζ) have been made. The first-order Taylor's approximations of f1(ζ) and f2(ζ) about ζ and evaluated at ζ = 0 result in the following system of equations:

(f1(ζ)−pa,0)|ζ=0+(df1(ζ)dζ−pa,1)|ζ=0(ζ)+O[ζ2]=0

(f2(ζ)−pB,0)|ζ=0+(df2(ζ)dζ−pB,1)|ζ=0(ζ)+O[ζ2]=0

The above system is satisfied if the terms (fi(ζ)−pj,k)|ζ=0 and (dfi(ζ)dζ−pj,l)|ζ=0 equal zero such that the zeroth-order terms satisfy:

0=(f1(ζ)−pa,0)|ζ=0

0=(f2(ζ)−pB,0)|ζ=0.

Solving for pa,0 and pB,0 simultaneously gives all zeroth-order equilibria for the deleterious allele and epiallele frequencies. The system for the first-order terms is

0=(df1(ζ)dζ−pa,1)|ζ=0

0=(df2(ζ)dζ−pB,1)|ζ=0

Letting (pa,0,i,pB,0,i) represent an arbitrary equilibrium for the zeroth-order terms, substituting this into the above equations, and solving for pa,1andpB,1 simultaneously gives the first-order terms that satisfy the multivariate equilibrium condition. Therefore, a first-order approximation is

p^a,i≈pa,0,i+pa,1,iζ+O[ζ2],

p^B,i≈pB,0,i+pB,1,iζ+O[ζ2],

p^A,i≈1−p^a,i−p^B,i+O[ζ2].

Original parameters are then substituted back that were assumed to be of order ζ, ni~ζ = ni, to obtain the first-order equilibrium approximation.

A2: local stability analysis

Recursion equations for the deleterious allele and epiallele frequencies are

pa′=f1(pa,pB)

pB′=f2(pa,pB).

The Jacobian matrix for the recursion equations is

(∂(f1(pa,pB))∂pa∂(f1(pa,pB))∂pB∂(f2(pa,pB))∂pa∂(f2(pa,pB))∂pB)

The eigenvalues of this Jacobian matrix evaluated at an equilibrium indicate its local stability (Otto and Day 2007). Let an eigenvalue be approximated as λ≈λ0+λ1ζ + λ2ζ2 …Then, the characteristic equation g(ζ) is

g(ζ)=det((∂(f1(pa,pB))∂pa∂(f1(pa,pB))∂pB∂(f2(pa,pB))∂pa∂(f2(pa,pB))∂pB)|(pa=pa,0+pa,1ζ,pB=pB,0+pB,1ζ,ni=ni~ζ,…)−(λ0+λ1ζ+λ2ζ2…00λ0+λ1ζ+λ2ζ2…))

Performing a Taylor series with respect to ζ around the point (ζ=0) on the characteristic equation approximated to the first order is

g(ζ)|ζ=0+(dg(ζ)dζ)|ζ=0(ζ)+O[ζ2]=0

The characteristic equation has the form a(λ0) + b(λ0, λ1) ζ = 0. Solving the zeroth-order term [a(λ0)] for eigenvalues gives both zeroth-order eigenvalues (λ0,1, λ0,2). Separately substituting both these solutions into the first-order term and then solving for the first-order term give both first-order terms (λ1,1, λ1,2). This is followed by substituting in the original parameters that were assumed to be of order (ζ), ni~ζ = ni, giving the approximated eigenvalues

λ(1)≈λ0,1+λ1,1+O[ζ2]

λ(2)≈λ0,2+λ1,2+O[ζ2]

The approximated eigenvalues can then be used to determine local stability such that |λ(1)| < 1 and |λ(2)| < 1 (Otto and Day 2007).

A3: example of an equilibrium calculation and local stability analysis:

To demonstrate the analysis procedure, go through the above methods to determine the equilibria and perform a local stability analysis on case A as an example. First, following assumptions of this case, letting the paramutation rate equal zero (m = 0), and inputting the fitness assumptions of Table 1 into the original recursion equations give the recursion equations for the deleterious allele and epiallele, respectively (see Supplementary File S2, case A for full recursion equations):

pa′=f1(pa,pB)

pB′=f2(pa,pB)

Then, setting the equilibrium condition (see Supplementary File S2, case A for full equations),

f1(p^a,p^B)−p^a=0

f1(p^a,p^B)−p^B=0

It is assumed that spontaneous epimutation and mutation are of order zeta, such that t1=t1~ζ,t2=t2~ζ,u=u~ζ, =z~ζ, c=c~ζ and equilibria are

p^a=pa,0+pa,1ζ+…

p^B=pB,0+pB,1ζ+…

Next, a Taylor series is performed with respect to ζ around the point (ζ=0) on both equations separately, and including up to the first-order terms. The equilibrium condition for the zeroth-order terms is

(−pa,0+−((−1+s)pa,02)+(1−hs−h2s2+s3)pa,0pB,0+(−1+hs)pa,0(−1+pa,0+pB,0)1−2hspa,0−spa,02+2hspa,02−2h2s2pB,0+2s3pa,0pB,0−s2pB,02+2h2s2pB,02)=0

(−pB,0−((−1+h2s2)pB,0(−1+pa,0+pB,0))+pB,0((−1+hs+h2s2−s3)pa,0+(−1+s2)pB,0)−1+2hspa,0+spa,02−2hspa,02+2h2s2pB,0−2s3pa,0pB,0+s2pB,02−2h2s2pB,02)=0

Solving for pa,0, and pB,0 gives three sets of zeroth-order terms. However, after further analysis, only one leads to a biologically valid equilibrium, (pa,0=0,pB,0=0). Substituting these values for the zeroth-order frequencies into the first-order terms (see Supplementary File S2, case A) for both equations gives the equilibrium condition for the first-order terms,

u−pa,1+(1−hs)pa,1=0

t1−pB,1−(−1+h2s2)pB,1=0

Solving for pa,1 and pB,1 and incorporating with the zeroth-order terms give the first-order approximation for the allele and epiallele frequencies:

p^a≈u~hsζ+O[ζ2]

p^B≈t1~h2s2ζ+O[ζ2]

p^A≈1−u~hsζ−t1~h2s2ζ+O[ζ2]

Then, back substituting the original parameters of order (ζ), (t1~ζ=t1,u~ζ=u), gives the first-order approximation.

Next, a local stability analysis is performed for this equilibrium. The Jacobian matrix is calculated and expressions for the approximated equilibrium terms and the mutation and spontaneous epimutation rates are substituted and assumed to be of order (ζ), such that pB=pB,0+pB,1ζ,pa=pa,0+pa,1ζ,t1=t1~ζ,t2=t2~ζ,u=u~ζ, z=z~ζ, c=c~ζ. After the substitutions, the characteristic equation was calculated, and then, its first-order Taylor's approximation with respect to ζ around the point (ζ = 0) was evaluated. Next, substituting in explicit equilibrium terms for this equilibrium (pa,0=0,pB,0=0,pa,1=uhs,pB,1=t1h2s2) results in the first-order approximation for the characteristic equation. Solving

(1−hs−λ0)(1−h2s2−λ0)=0

for the zeroth-order terms gives

λ0,1=1−hs,λ0,2=1−h2s2.

Substituting λ0,1=1−hs into the first-order term above and solving for λ1,1 gives

λ1,1=(5−2h−3hs)u+(−1+hs)z+(1−2hs+s3h2s2)t1

Then, making the substitution λ0,2=1−h2s2 into the first-order term and solving for λ1,2, results in

λ1,2=−c+u+us3hs+5t1−2t1h2−t2+h2s2(c−2u−3t1+t2)+O[ζ2]

Bringing the zero- and first-order terms together and simplifying give the first-order approximation for the eigenvalues

λ(1)≈1−hs+(5−2h−3hs)u+(−1+hs)z+(1−2hs+s3h2s2)t1+O[ζ2]

λ(2)≈1−h2s2−c+u+us3hs+5t1−2t1h2−t2+h2s2(c−2u−3t1+t2)+O[ζ2]

Note: All analysis was performed using Wolfram Mathematica, version 12.3.1.0.
==== Refs
Literature cited

Barber  WT, Zhang  W, Win  H, Varala  KK, Dorweiler  JE, Hudson  ME, Moose  SP. 2012. Repeat associated small RNAs vary among parents and following hybridization in maize. Proc Natl Acad Sci U S A. 109 (26 ):10444–10449. doi:10.1073/pnas.1202073109.22689990
Bell  AM, Hellmann  JK. 2019. An integrative framework for understanding the mechanisms and multigenerational consequences of transgenerational plasticity. Ann Rev Ecol Evol Syst. 50 (1 ):97–118. doi:10.1146/annurev-ecolsys-110218-024613.36046014
Bente  H, Foerster  AM, Lettner  N, Scheid  OM. 2021. Polyploidy-associated paramutation in Arabidopsis is determined by small RNAs, temperature, and allele structure. PLoS Genet. 17 (3 ):e1009444. doi:10.1371/journal.pgen.1009444.33690630
Bonduriansky  R, Day  T. 2018. Extended Heredity: A New Understanding of Inheritance and Evolution. Princeton: Princeton University Press.
Brink  RA . 1956. A genetic change associated with the R locus in maize which is directed and potentially reversible. Genetics. 41 (6 ):872–889. doi:10.1093/genetics/41.6.872.17247669
Capovilla  M, Robichon  A, Rassoulzadegan  M. 2017. A new paramutation-like example at the delta gene of Drosophila. PLoS One. 12 (3 ):e0172780. doi:10.1371/journal.pone.0172780.28355214
Charlesworth  D . 2006. Evolution of plant breeding systems. Curr Biol. 16 (17 ):R726–R735. doi:10.1016/j.cub.2006.07.068.16950099
Charlesworth  D, Barton  NH, Charlesworth  B. 2017. The sources of adaptive variation. Proc Biol Sci. 284 (1855 ):20162864. doi:10.1098/rspb.2016.2864.28566483
Choi  JY, Yuh  CGL. 2020. Double-edged sword: the evolutionary consequences of the epigenetic silencing of transposable elements. PLoS Genet. 16 (7 ):e1008872. doi:10.1371/journal.pgen.1008872.32673310
Clark  AG . 1998. Mutation-selection balance with multiple alleles. In: Woodruff  RC, Thompson  JN, editors. Mutation and Evolution. Dordrecht: Springer Netherlands (Contemporary Issues in Genetics and Evolution, vol. 7 ). p. 41–47. doi:10.1007/978-94-011-5210-5_4.
Coe  EH . 1966. The properties, origin, and mechanism of conversion-type inheritance at the B locus in maize. Genetics. 53 (6 ):1035–1063. doi:10.1093/genetics/53.6.1035.17248307
Crow  JF, Kimura  M. 2010. An Introduction to Population Genetics Theory. Indian ed. Jodhpur (NJ): Scientific Publisher (India); The Blackburn Press. p. 255–303.
Cubas  P, Vincent  C, Coen  E. 1999. An epigenetic mutation responsible for natural variation in floral symmetry. Nature. 401 (6749 ):157–161. doi:10.1038/43657.10490023
Day  T, Bonduriansky  R. 2011. A unified approach to the evolutionary consequences of genetic and nongenetic inheritance. Am Nat. 178 (2 ):E18–E36. doi:10.1086/660911.21750377
Dorador  AP, Dalikova  M, Cerbin  S, Stillman  CM, Zych  MG, Hawley  RS, Miller  DE, Ray  DA, Funikov  SY, Evgen'ev  MB, et al  2022. Paramutation-like epigenetic conversion by piRNA at the telomere of Drosophila virilis. Biology (Basel). 11 (10 ):1480. doi:10.3390/biology11101480.36290385
Furrow  RE . 2014. Epigenetic inheritance, epimutation, and the response to selection. PLoS One. 9 (7 ):e101559. doi:10.1371/journal.pone.0101559.25019291
Gao  Y, Zhao  Y. 2013. Epigenetic suppression of T-DNA insertion mutants in Arabidopsis. Mol Plant. 6 (2 ):539–545. doi:10.1093/mp/sss093.22973063
Geoghegan  JL, Spencer  HG. 2012. Population-epigenetic models of selection. Theor Popul Biol. 81 (3 ):232–242. doi:10.1016/j.tpb.2011.08.001.21855559
Geoghegan  JL, Spencer  HG. 2013a. Exploring epiallele stability in a population-epigenetic model. Theor Popul Biol. 83 (February ):136–144. doi:10.1016/j.tpb.2012.09.001.23044385
Geoghegan  JL, Spencer  HG. 2013b. The adaptive invasion of epialleles in a heterogeneous environment. Theor Popul Biol. 88 (September ):1–8. doi:10.1016/j.tpb.2013.05.001.23689020
Geoghegan  JL, Spencer  HG. 2013c. The evolutionary potential of paramutation: a population-epigenetic model. Theor Popul Biol. 88 (September ):9–19. doi:10.1016/j.tpb.2013.05.003.23751219
Glémin  S . 2010. Surprising fitness consequences of GC-biased gene conversion: i. Mutation load and inbreeding depression. Genetics. 185 (3 ):939–959. doi:10.1534/genetics.110.116368.20421602
Gouil  Q, Baulcombe  DC. 2018. Paramutation-like features of multiple natural epialleles in tomato. BMC Genomics. 19 (1 ):203. doi:10.1186/s12864-018-4590-4.29554868
van der Graaf  A, Wardenaar  R, Neumann  DA, Taudt  A, Shaw  RG, Jansen  RC, Schmitz  RJ, Colomé-Tatché  M, Johannes  F. 2015. Rate, spectrum, and evolutionary dynamics of spontaneous epimutations. Proc Natl Acad Sci U S A. 112 (21 ):6676–6681. doi:10.1073/pnas.1424254112.25964364
Greenspoon  PB, Spencer  HG. 2018. The evolution of epigenetically mediated adaptive transgenerational plasticity in a subdivided population. Evolution. 72 (12 ):2773–2780. doi:10.1111/evo.13619.30298912
Groszmann  M, Greaves  IK, Albertyn  ZI, Scofield  GN, Peacock  WJ, Dennis  ES. 2011. Changes in 24-Nt siRNA levels in Arabidopsis hybrids suggest an epigenetic contribution to hybrid vigor. Proc Natl Acad Sci U S A. 108 (6 ):2617–2622. doi:10.1073/pnas.1019217108.21266545
He  G, Zhu  X, Elling  AA, Chen  L, Wang  X, Guo  L, Liang  M, He  H, Zhang  H, Chen  F, et al  2010. Global epigenetic and transcriptional trends among two rice subspecies and their reciprocal hybrids. Plant Cell. 22 (1 ):17–33. doi:10.1105/tpc.109.072041.20086188
Hollick  JB . 2017. Paramutation and related phenomena in diverse species. Nat Rev Genet. 18 (1 ):5–23. doi:10.1038/nrg.2016.115.27748375
van Houwelingen  A, Souer  E, Mol  J, Koes  R. 1999. Epigenetic interactions among three dTph1 transposons in two homologous chromosomes activate a new excision-repair mechanism in petunia. Plant Cell. 11 (7 ):1319–1336. doi:10.1105/tpc.11.7.1319.10402432
Jacobsen  SE, Meyerowitz  EM. 1997. Hypermethylated SUPERMAN epigenetic alleles in Arabidopsis. Science. 277 (5329 ):1100–1103. doi:10.1126/science.277.5329.1100.9262479
Johannes  F, Porcher  E, Teixeira  FK, Saliba-Colombani  V, Simon  M, Agier  N, Bulski  A, Albuisson  J, Heredia  F, Audigier  P, et al  2009. Assessing the impact of transgenerational epigenetic variation on complex traits. PLoS Genet. 5 (6 ):e1000530. doi:10.1371/journal.pgen.1000530.19557164
Johannes  F, Schmitz  RJ. 2019. Spontaneous epimutations in plants. New Phytol. 221 (3 ):1253–1259. doi:10.1111/nph.15434.30216456
Khaitová  LC, Fojtová  M, Křížová  K, Lunerová  J, Fulneček  J, Depicker  A, Kovařík  A. 2011. Paramutation of tobacco transgenes by small RNA-mediated transcriptional gene silencing. Epigenetics. 6 (5 ):650–660. doi:10.4161/epi.6.5.15764.21521939
Klironomos  FD, Berg  J, Collins  S. 2013. How epigenetic mutations can affect genetic evolution: model and mechanism: problems & paradigms. BioEssays. 35 (6 ):571–578. doi:10.1002/bies.201200169.23580343
Kronholm  I, Collins  S. 2016. Epigenetic mutations can both help and hinder adaptive evolution. Mol Ecol. 25 (8 ):1856–1868. doi:10.1111/mec.13296.26139359
Li  L, Petsch  K, Shimizu  R, Liu  S, Xu  WW, Ying  K, Yu  J, Scanlon  MJ, Schnable  PS, Timmermans  MC, et al  2013. Mendelian and non-Mendelian regulation of gene expression in maize. PLoS Genet. 9 (1 ):e1003202. doi:10.1371/journal.pgen.1003202.23341782
Luff  B, Pawlowski  L, Bender  J. 1999. An inverted repeat triggers cytosine methylation of identical sequences in Arabidopsis. Mol Cell. 3 (4 ):505–511. doi:10.1016/S1097-2765(00)80478-5.10230403
Manning  K, Tör  M, Poole  M, Hong  Y, Thompson  AJ, King  GJ, Giovannoni  JJ, Seymour  GB. 2006. A naturally occurring epigenetic mutation in a gene encoding an SBP-box transcription factor inhibits tomato fruit ripening. Nat Genet. 38 (8 ):948–952. doi:10.1038/ng1841.16832354
Martin  A, Troadec  C, Boualem  A, Rajab  M, Fernandez  R, Morin  H, Pitrat  M, Dogimont  C, Bendahmane  A. 2009. A transposon-induced epigenetic change leads to sex determination in melon. Nature. 461 (7267 ):1135–1138. doi:10.1038/nature08498.19847267
Mittelsten Scheid  O, Afsar  K, Paszkowski  J. 2003. Formation of stable epialleles and their paramutation-like interaction in tetraploid Arabidopsis thaliana. Nat Genet. 34 (4 ):450–454. doi:10.1038/ng1210.12847525
Miura  K, Agetsuma  M, Kitano  H, Yoshimura  A, Matsuoka  M, Jacobsen  SE, Ashikari  M. 2009. A metastable DWARF1 epigenetic mutant affecting plant stature in rice. Proc Natl Acad Sci U S A. 106 (27 ):11218–11223. doi:10.1073/pnas.0901942106.19541604
Ong-Abdullah  M, Ordway  JM, Jiang  N, Ooi  SE, Kok  SY, Sarpan  N, Azimi  N, Hashim  AT, Ishak  Z, Rosli  SK, et al  2015. Loss of karma transposon methylation underlies the mantled somaclonal variant of oil palm. Nature. 525 (7570 ):533–537. doi:10.1038/nature15365.26352475
Otto  SP, Day  T. 2007. A Biologist's Guide to Mathematical Modelling in Ecology and Evolution. Princeton (NJ): Princeton University Press.
Pilu  R, Panzeri  D, Cassani  E, Cerino Badone  F, Landoni  M, Nielsen  E. 2009. A paramutation phenomenon is involved in the genetics of maize low phytic acid1-241 (Lpa1-241) trait. Heredity (Edinb). 102 (3 ):236–245. doi:10.1038/hdy.2008.96.18781168
Quadrana  L, Almeida  J, Asís  R, Duffy  T, Dominguez  PG, Bermúdez  L, Conti  G, Corrêa  da Silva JV, Peralta  IE, Colot  V, et al  2014. Natural occurring epialleles determine vitamin e accumulation in tomato fruits. Nat Commun. 5 (1 ):4027. doi:10.1038/ncomms5027.
Quadrana  L, Colot  V. 2016. Plant transgenerational epigenetics. Ann Rev Genet. 50 (1 ):467–491. doi:10.1146/annurev-genet-120215-035254.27732791
Rassoulzadegan  M, Grandjean  V, Gounon  P, Vincent  S, Gillot  I, Cuzin  F. 2006. RNA-mediated non-Mendelian inheritance of an epigenetic change in the mouse. Nature. 441 (7092 ):469–474. doi:10.1038/nature04674.16724059
Shen  H, He  H, Li  J, Chen  W, Wang  X, Guo  L, Peng  Z, He  G, Zhong  S, Qi  Y, et al  2012. Genome-wide analysis of DNA methylation and gene expression changes in two Arabidopsis ecotypes and their reciprocal hybrids. Plant Cell. 24 (3 ):875–892. doi:10.1105/tpc.111.094870.22438023
Sidorenko  LV, Chandler  VL, Wang  X, Peterson  T. 2024. Transcribed enhancer sequences are required for maize P1 paramutation. Genetics. 226 (1 ):iyad178. doi:10.1093/genetics/iyad178.38169343
Silveira  AB, Trontin  C, Cortijo  S, Barau  J, Del Bem  LE, Loudet  O, Colot  V, Vincentz  M. 2013. Extensive natural epigenetic variation at a de novo originated gene. PLoS Genet. 9 (4 ):e1003437. doi:10.1371/journal.pgen.1003437.23593031
Smithson  MW, Dybdahl  MF, Nuismer  SL. 2019. The adaptive value of epigenetic mutation: limited in large but high in small peripheral populations. J Evol Biol. 32 (12 ):1391–1405. doi:10.1111/jeb.13535.31529541
Soppe  WJJ, Jacobsen  SE, Alonso-Blanco  C, Jackson  JP, Kakutani  T, Koornneef  M, Peeters  AJM. 2000. The late flowering phenotype of fwa mutants is caused by gain-of-function epigenetic alleles of a homeodomain gene. Mol Cell. 6 (4 ):791–802. doi:10.1016/s1097-2765(05)00090-0.11090618
Springer  NM, McGinnis  KM. 2015. Paramutation in evolution, population genetics and breeding. Semin Cell Dev Biol. 44 (August ):33–38. doi:10.1016/j.semcdb.2015.08.010.26325077
Stenøien  HK, Pedersen  B. 2005. Mutation and epimutation load in haploid and diploid life forms. J Theor Biol. 233 (1 ):119–126. doi:10.1016/j.jtbi.2004.09.012.15615625
Uller  T, English  S, Pen  I. 2015. When is incomplete epigenetic resetting in germ cells favoured by natural selection?  Proc Biol Sci. 282 (1811 ):20150682. doi:10.1098/rspb.2015.0682.26136447
Vastenhouw  NL, Brunschwig  K, Okihara  KL, Müller  F, Tijsterman  M, Plasterk  RHA. 2006. Long-term gene silencing by RNAi. Nature. 442 (7105 ):882–882. doi:10.1038/442882a.16929289
Webster  AK, Phillips  PC. 2024. Heritable epigenetic variation facilitates long-term maintenance of epigenetic and genetic variation. G3 (Bethesda). 14 (2 ):jkad287. doi:10.1093/g3journal/jkad287.38113034
Xue  W, Ruprecht  C, Street  N, Hematy  K, Chang  C, Frommer  WB, Persson  S, Niittylä  T. 2012. Paramutation-like interaction of T-DNA loci in Arabidopsis. PLoS One. 7 (12 ):e51651. doi:10.1371/journal.pone.0051651.23272131
Yi  SV, Goodisman  MAD. 2021. The impact of epigenetic information on genome evolution. Philos Trans R Soc B: Biol Sci. 376 (1826 ):20200114. doi:10.1098/rstb.2020.0114.
