Skip to main content
  • Research article
  • Open access
  • Published:

Paramutation-like features of multiple natural epialleles in tomato



Freakish and rare or the tip of the iceberg? Both phrases have been used to refer to paramutation, an epigenetic drive that contravenes Mendel’s first law of segregation. Although its underlying mechanisms are beginning to unravel, its understanding relies only on a few examples that may involve transgenes or artificially generated epialleles.


By using DNA methylation of introgression lines as an indication of past paramutation, we reveal that the paramutation-like properties of the H06 locus in hybrids of Solanum lycopersicum and a range of tomato relatives and cultivars depend on the timing of sRNA production and conform to an RNA-directed mechanism. In addition, by scanning the methylomes of tomato introgression lines for shared regions of differential methylation that are absent in the S. lycopersicum parent, we identify thousands of candidate regions for paramutation-like behaviour. The methylation patterns for a subset of these regions segregate with non Mendelian ratios, consistent with secondary paramutation-like interactions to variable extents depending on the locus.


Together these results demonstrate that paramutation-like epigenetic interactions are common for natural epialleles in tomato, but vary in timing and penetrance.


Paramutation is an epigenetic process in plants (including pea, maize, tomato [1]) and animals (worm [2], fruit fly [3], mouse [4]) that is associated with gene silencing. It is unlike other epigenetic mechanisms, however, in that it involves transfer of the silent state from an allele with epigenetic modification to its active homologue. This paramutated allele then becomes silenced and it acquires the ability to silence other active alleles in subsequent generations so that inheritance patterns are non Mendelian [1, 5].

The best characterised examples of paramutation are from maize, at the b1, r1, and pl1 loci. Genetic screens have implicated NRPD1 (rmr6 [6]), the major subunit of Pol IV; the RNA-dependent RNA polymerase RDR2 (mop1 [7]) and NRPD2a (mop2 /rmr7 mutants [8, 9]) the shared subunit of Pol IV and Pol V. These proteins are all required for RNA-directed DNA methylation (RdDM) in which DNA methytransferases are guided to the target sequence in the genome by base pairing of small RNAs (sRNAs). RdDM is also associated with paramutation of the SULFUREA locus in tomato [10, 11] and trans-chromosomal DNA methylation in Arabidopsis thaliana hybrids [12, 13]. Based on these findings the dominant model of paramutation implicates RdDM in the establishment and/or maintenance of the epigenetic mark.

Our interest in paramutation follows from an earlier study of sRNA in tomato lines in which homozygous regions were introgressed from Solanum pennellii into Solanum lycopersicum (cv. M82) [14]. The resulting introgression lines (ILs) each carry many loci at which sRNAs are more abundant than in either parent: they are transgressive [15]. To explain these findings we invoked epigenetic mechanisms because, in some instances, there was hypermethylation of the genomic DNA corresponding to the sRNA locus.

In this present study we focus initially on one locus, H06, at which there is transgressive sRNA and DNA hypermethylation in multiple introgression lines [15]. We were prompted to consider the involvement of a paramutation-like process at H06 because presence of this epiallele in multiple ILs was not consistent with Mendelian inheritance. The ILs are produced by recurrent backcrossing of the F1 hybrid to the Solanum lycopersicum cv. M82 parent so that Mendelian features of the hybrid genome would co-segregate with specific regions of introgressed DNA. At H06, however, the segregation must have been independent of any introgressed regions. We could rule out that the anomalous behaviour of H06 was due to a spontaneous epigenetic change because we could reproducibly recapitulate the transgressive sRNA and DNA hypermethylation at this locus in F2 progeny of the S. lycopersicum × S. pennelli cross [15].

An alternative explanation of the H06 epiallele invoked non-Mendelian inheritance due to a hybrid-induced ‘paramutation’ that, once triggered, would be inherited in the recurrent backcross progeny independently of any one region in the S. pennelli genome. The genetic and molecular tests presented here are fully consistent with that hypothesis and they indicate further that the timing of H06 paramutation correlates with the production of sRNAs from the paramutagenic epiallele. We also identify other methylated DNA epialleles in the ILs with paramutation-like properties: they are absent from the parental lines, present in multiple ILs independently of a specific introgressed region from S. pennelli and they transfer their epigenetic mark to a non-methylated allele following backcrossing to M82. Based on the characterisation of these loci we propose that multiple paramutation-like events occur in the progeny of a S. lycopersicum cv. M82 × S. pennellii cross. These events illustrate how epigenetic effects, of which paramutation is an extreme example, may be induced by hybridisation of divergent genomes. Further characterisation of this epigenetic spectrum will reveal the defining features of paramutation loci and the extent to which hybrid-induced changes to the epigenome can influence transgressive segregation in crop plants and natural evolution.


Paramutation at the H06 locus in multiple lines

The H06 locus [15] is in the euchromatin of chromosome 8, nine kilobases upstream of two genes in divergent orientation (Fig. 1a). It is unmethylated and lacks small RNAs in M82 but, in IL8–3, it is methylated in all three contexts (CG, CHG and CHH, where H is any nucleotide but G) and produces abundant 24-nt sRNAs (Fig. 1ac). We refer to this as the H06IL epiallele. The homologous locus in Solanum lycopersicum cv. Micro-Tom and Solanum pimpinellifolium has an epigenetic profile similar to H06IL both in terms of DNA methylation and sRNA production (Fig. 1b and c). In Solanum pennellii, however, this locus was distinct from the other species in that there was full methylation at CG and CHG but no sRNAs and only partial CHH methylation in the corresponding region (Fig. 1b and c). This finding is different from our previous report in which the S. pennellii locus was described as hypomethylated at all three C contexts [15] (in the aerial part of two-week-old seedlings, whereas here we use only the leaf of a two-week-old seedling). The present data are, however, consistent with the previous finding in that the RdDM features of H06 - mCHH and sRNA - are low or absent in S. pennellii seedlings.

Fig. 1
figure 1

H06 epialleles. a Genomic location of H06, sRNAs in M82 and IL8–3 whole seedlings and methylation in seedling leaves. The H06IL epiallele in IL8–3 was methylated in the CG, CHG and CHH contexts and the source of abundant sRNAs. Neither DNA methylation nor sRNAs were detected in the H06M82 epiallele. b High sRNA production at H06 in seedlings of three introgression lines, the Micro-Tom cultivar and the wild relative S. pimpinellifolium. c H06 methylation in leaves of M82, S. pennellii, IL8–3, Micro-Tom and S. pimpinellifolium determined by Sanger bisulfite sequencing (at least 10 independent clones per genotype). Primer sequences are given in Additional file 10, and H06 sequences in Additional file 11

To find out whether the H06IL has properties consistent with paramutation, we crossed M82 × IL8–3 and monitored the DNA methylation in F1 (M82 × IL8–3), F2 and BC1 (M82 × F1) generations. The DNA was extracted from the leaves of 15-day-old plants and it was assayed with McrBC digestion. The first aim of these tests was to establish whether the epigenetic mark could be heritably transferred from H06IL to H06M82 to create an H06IL’ epiallele. The second aim was to find out whether H06IL’ could mediate secondary paramutation and transfer its epigenetic mark to H06M82.

Fourteen of the fifteen F1 plants had highly methylated H06 DNA in this assay (Fig. 2a) indicating that H06M82 had been converted to H06IL’. The level of methylation in the 30 F2 plants derived from highly methylated F1s was always high, showing that H06IL and the newly established H06IL’ are normally stable, thus fulfilling the first criterion of paramutation that the epigenetic mark could be transferred from H06IL to H06M82 to form a stable H06IL’ epiallele. There was, however, instability of H06IL in a single F1 plant and its F2 progeny in which the H06 DNA was completely unmethylated (Fig. 2a).

Fig. 2
figure 2

H06 paramutation assessed by McrBC. McrBC digests methylated DNA, so that comparing amplification of McrBC-treated and non-treated DNA by qPCR reveals the proportion of molecules that were methylated. a H06 McrBC in M82 × IL8–3 seedling leaves. Unmethylated H06M82 (n = 12), crossed with the methylated H06IL from IL8–3 (n = 7), became methylated in most F1s (14 out of n = 15), and remained methylated in the F2 (n = 30). Back-crossing the F1 with M82 demonstrated weak paramutagenicity (BC1, n = 26, p-value = 0.03776, one-sided binomial test). Paramutation also occurred with the IL as the female parent (IL8–3 × M82, n = 1). b H06 McrBC in M82 × Micro-Tom and S.pimpinellifolium. F1 seedling leaves are fully methylated (n = 6 for each cross), like Micro-Tom and S. pimpinellifolium themselves (n = 4 and 5 respectively). (C) H06 McrBC in M82 × S. pennellii. The F1 vegetative tissue was half-methylated (n = 4) whereas the F1 pollen and the F2 leaves were fully methylated (n = 5 and 11 respectively)

The secondary paramutation criterion was tested in the M82 × F1 backcrossed (BC) progeny. If H06IL but not H06IL’ is paramutagenic then half of the BC1 plants would have fully methylated H06 DNA and half would have at most 50% methylated DNA. In contrast, if both epialleles were paramutagenic, there would be more than half of the plants with full methylation of H06 DNA. The data (Fig. 2a) are consistent with secondary paramutation because there was an excess of highly methylated plants (18 highly methylated plants, 8 half-methylated or lower, p-value = 0.038, one-sided binomial test). It is likely, however, that some of the H06IL or H06IL’ alleles in the BC1 were unstable and reverted to H06M82 because five of these BC plants had less than 50% methylated DNA at H06.

To investigate the paramutagenic properties of the H06 allele in Micro-Tom, S. pimpinellifolium and S. pennellii we made a further series of crosses with M82 and analysed the F1 progeny. With Micro-Tom and S. pimpinellifolium the results were as with IL8–3: all of the H06 alleles were highly methylated in the F1, irrespective of the direction of the cross (Fig. 2b). In contrast, in M82 × S. pennellii, about half of the H06 alleles were methylated in the 15-day leaves (Fig. 2c). This level increased in flowers and, in pollen and 15 d leaves of F2 plants, it was close to 100% (Fig. 2c). From these data we conclude that our various Solanum genotypes carried three distinct epialleles at H06: H06M82, the S. pennellii allele referred to as H06pen and H06IL. The H06M82 allele had neither methylation nor sRNAs. H06pen had DNA methylation but without the sRNA characteristics of RdDM in leaves of seedlings. This allele triggered an epigenetic change to H06IL but only after reproductive development in the F1. The third epiallele - H06IL - was present in Micro-Tom and S. pimpinellifolium in addition to IL8–3, had both abundant sRNAs and DNA methylation in leaves of seedlings, and it could trigger a paramutation-like change to H06M82 without any evident lag. The H06IL’ alleles had the general characteristics of H06IL but with incomplete penetrance of its effect on H06M82.

H06 paramutation timing correlates with sRNA production

To further investigate the interaction of H06M82 with H06pen, we used allele-specific Sanger bisulfite sequencing of H06 DNA from leaves and pollen of M82, S. pennellii and their F1. Consistent with the McrBC results, the allelic methylation of the F1 in the leaves mirrored that of the parents: low methylation of the M82 allele and high methylation of the S. pennellii allele (Fig. 3a) at CG and CHG. In pollen, however, the M82 allele of the F1 became hypermethylated in CG, CHG and CHH contexts (Fig. 3b), whereas pollen in the M82 parent is hypomethylated.

Fig. 3
figure 3

Paramutation in the flowers of the M82 × S. pennellii F1. H06 methylation in leaf (a) and pollen (b) of M82, S. pennellii and the F1. At least 10 independent Sanger bisulfite clones were used for each plot. c H06 sRNAs are produced in S. pennellii and F1 flowers, although at low levels (normalised counts in libraries of about 10 million reads). Two M82, two S. pennellii, and four F1 libraries for each tissue

This gain of methylation by the M82 allele in the F1 correlated with an increase in sRNA production at H06 (Fig. 3c). The H06 sRNA levels were highest in flowers of S. pennelli and the F1, and very low in M82 leaves, flowers and pollen. These levels were, however, much lower than in IL8–3: in the flowers of S. pennellii there were 1–2 reads per million mapped whereas in IL8–3 seedlings there were more than 150.

From these data there is a clear correlation of sRNA with the paramutation-like properties of the various H06 loci. H06M82 lacks both sRNA and DNA methylation and is paramutable; H06pen with low levels of sRNAs triggers a delayed paramutation-like process and H06IL, at which the sRNA levels are high, mediates a rapid transfer of the epigenetic marks to H06M82. This correlation of paramutation and sRNA implicates RdDM in the process. Furthermore the changes in both mCHH and sRNA in pollen and flowers suggest that the reproductive phase may be a key developmental stage in the transfer of an epigenetic mark between alleles of H06.

De novo methylated H06 recapitulates paramutation

In principle the association of RdDM with the paramutation-like properties of H06 could be either a cause or consequence of the transition from H06M82 to H06IL. To test the possibility of a causal role of RNA we used virus-induced gene silencing (VIGS) in which an RNA virus is modified to carry a small host genomic sequence insert. Such RNA viruses may direct DNA methylation of the corresponding genomic DNA of the infected plant [16] and we have previously used this approach to recapitulate the silencing of the sulf locus in tomato [11].

To test this system at H06, we infected unmethylated M82 plants with a Tobacco Rattle Virus (TRV) containing a 500-bp H06 sequence (TRV-H06, Fig. 4a). Out of 11 successfully infected plants (V0 generation) spread across three replicate experiments, only one plant gave rise to methylated progeny (V1 generation) (Fig. 4b), with 30% of the plants having high methylation (20/59). None of 6 the control plants infected with unmodified TRV produced methylated progeny (56 tested V1 plants) and, with many tens of M82 plants tested over a period of 3 years, we have never observed spontaneous methylation of H06 DNA. It is likely therefore that the new epiallele (H06VIGS) was triggered by the TRV-H06 VIGS.

Fig. 4
figure 4

Virus-induced paramutation. a Tobacco Rattle Virus modified to contain a 500-bp H06 fragment. b DNA methylation in the progeny of the V0 infected plant that gave methylated offspring, assessed by McrBC (n = 59 V1 plants). V1 methylation patterns were stable in the V2 generation (n = 4, 3, 4 and 4 for V2_1 to V2_4), and exhibited weak paramutagenicity in backcrosses to M82 (V1BC1, n = 10 for both V1_3 × M82 and V1_4 × M82). c H06 methylation of a V1 plant assessed by Sanger bisulfite sequencing (9 independent clones)

This new epiallele (H06VIGS, Fig. 4c) was distinct from H06pen and H06IL (Fig. 1c) in that the hypermethylation was in a restricted region of 200 bp rather than 500 bp. The H06VIGS epialleles were stable in the V2 generation (Fig. 4b). We tested whether this new epiallele was paramutagenic by backcrossing V1 plants to M82: if H06VIGS is inherited as a standard Mendelian locus without paramutation, the BC1 plants would have no more than 50% of methylated H06 DNA in an McrBC assay; if however H06VIGS is paramutagenic, some backcrossed plants would exhibit methylation above 50%. Of 21 BC1 plants, 4 had substantially more than 50% of methylated H06 DNA (Fig. 4b), and we conclude that H06VIGS has weak paramutagenic activity. There were also many plants with less than 50% of methylated H06 in these BC1 plants that is due, presumably, to instability of H06VIGS in the heterozygous condition. The ability to epigenetically modify and confer paramutation properties to H06 by VIGS indicates that sRNAs could be causal in paramutation. However the low efficiency of VIGS in methylating H06 precludes definitive conclusions about the mechanism of paramutation.

Genome-wide ‘paramutation’ in introgression lines?

Is H06 paramutation an isolated example, or could there be other similar loci in Solanum genomes? To address this question we screened the DNA methylome of M82 and three different introgression lines, IL1–1, IL2–5 and IL8–3, in duplicates (summarised in Additional file 1). We reasoned that paramutation loci would be differentially methylated compared to M82 in multiple introgression lines. To identify DMRs with this characteristic we compared methylation counts in CG, CHG and CHH of each introgression lines to M82 in 300-bp sliding windows. Each introgression line had 10,000–30,000 DMRs depending on context (CG, CHG and CHH) as compared to M82 (Fig. 5). Only a subset of these DMRs were shared in all three ILs: around 4600 CG hypermethylated DMRs, 1800 CG hypomethyated DMRs, 3600 CHG hypermethylated DMRs, 1000 CHG hypomethylated DMRs, 262 CHH hypermethylated DMRs and 900 CHH hypomethylated DMRs (Fig. 5 and Additional file 2). Of these DMRs there were 25, including H06, that were hypermethylated in all three contexts (Additional file 3). There were also 10 hypomethylated DMRs in all three contexts (Additional file 4).

Fig. 5
figure 5

Introgression lines DMRs. a Number of DMRs in each introgression line compared to M82 for CG, CHG and CHH contexts. Two plants per genotype. b DMR overlap across introgression lines. Hypermethylated DMRs have higher methylation in the introgression lines compared to M82. Area-proportional Venn diagrams

The RdDM model of paramutation suggests that there would be 24-nt sRNAs at the paramutated loci, as is the case for H06. Using seedling sRNA libraries from [15], we identified 134 loci upregulated for these sRNAs in the ILs, and 68 down-regulated (Fig. 6a and b). Of the 134 upregulated sRNA loci, 9 overlapped hypermethylated DMRs and 3 overlapped hypomethylated DMRs. In addition to H06, one other locus (ch12:62124801–62,125,100, referred to as Hyper1) had increased sRNAs and hypermethylation in all three contexts in the introgression lines (Additional file 5: Figure S1). This locus and the other hypermethylated DMRs that overlap differential sRNAs are prime candidates for paramutation that would follow the same pattern as H06.

Fig. 6
figure 6

Differential sRNA and gene expression in the introgression lines. a Differential analysis of sRNAs between M82 (two plants) and IL1–1/IL2–5/IL8–3 seedlings (one plant per genotype). Data from [15]. b Classification of differential sRNA loci by ShortStack. c Differential gene expression between M82 (n = 40) and IL1–1/IL2–5/IL8–3 seedlings (n = 30). Data from [17]. sRNA loci and gene that are differentially expressed (p-adj < 0.05) are coloured in red

DMRs can be associated with differences in gene expression, so we performed differential expression analysis between M82 and the three ILs under scrutiny with data from [17]. In this data set for seedlings grown in sun and shade, we found 124 genes that were upregulated in ILs compared to M82, and 108 that were downregulated (Fig. 6c). The reported log2 fold changes of the differentially expressed genes were modest, in part because a large proportion of differentially expressed genes had low average expression (and hence low counts, leading to a shrinkage of the log2 fold change estimates, Fig. 6c and Additional file 6: Figure S2). There was no significant enrichment in biological process, molecular function or cellular compartments of gene ontology terms associated with the differentially expressed genes (upregulated, dowrnegulated or the full list). 10 and 5 of the upregulated genes were overlapped (over their 2 kb promoter or transcribed sequence) by shared hypermethylated and hypomethylated DMRs, respectively, while 4 of the down-regulated genes overlapped hypermethylated DMRs and 7 hypomethylated ones. These genes may be of particular interest to investigate transcriptional and physiological effects of paramutation following a cross between M82 and S. pennellii. Although this data set encompasses two environmental conditions (sun and shade), DMRs may affect the transcription of more genes in a tissue- or development-specific manner.

Paramutation-like activity of DMRs

Hypermethylated DMRs

To evaluate the paramutagenic properties of the hypermethylated DMRs, we selected 9 regions based on hypermethylation in multiple cytosine contexts in the three ILs, amenability to McrBC assays, presence of differential sRNAs, and proximity to differentially expressed genes (Hyper1–9, Table 1, Additional file 5: Figure S1). We tested their inheritance in the progeny of M82 × IL8–3 by McrBC-qPCR. Mendelian segregation would give 50% methylation in the F1, and in the F2 either 100% (1/4 progeny), 50% (1/2 progeny) or 0% (1/4 progeny) methylation. We tested these DMRs in 10 different F1s and in 33 F2s from three different F1 plants. Any evidence for more than 50% methylation of the F1 alleles or more than 3/4 F2 progeny with 50% or higher methylation would indicate non-Mendelian inheritance and evidence for a paramutation-like mechanism.

Table 1 Summary of paramutation candidates and their validation

According to these criteria, there was evidence of a paramutation-like process at four loci (Hyper1–4, Fig. 7a and Table 1). For three of those (Hyper1–3), some F1s had more than 50% of methylated alleles. At the fourth locus (Hyper4, Fig. 7a) the F1 progeny had approximately 50% of methylated alleles consistent with Mendelian inheritance but more than 3/4 of the F2 progeny had more than 50% methylated alleles (only 3/33 F2s had low methylation, p-value = 0.021, one-sided binomial test). This pattern is compatible with a partial gain of methylation by the M82 allele later in the development of the F1 (as was the case for H06 in M82 × S. pennellii). Of note, Hyper1 corresponds to a region of highly transgressive sRNAs (log2 fold change of 7) that are 24-nt in length (Additional file 5: Figure S1). It is therefore very similar in behaviour to H06. As for Hyper3, it is located in the immediate promoter of the Solyc09g064640 gene, which is expressed at lower levels in the introgression lines compared to M82 (Additional file 6: Figure S2).

Fig. 7
figure 7

Validation of paramutation by McrBC. a Hypermethylated DMRs. b Hypomethylated DMRs. A solid frame indicates strong evidence for paramutation, and a dashed frame partial evidence. For each region, the left hand panel shows the results of methylation analysis by McrBC for individual plants (10 M82, 10 IL8–3, 10 F1, 3 × 10 to 11 F2). This information is collated in the right hand panel with the splitting of the F2s according to their F1 parent. Low methylation: ≤ 33%. Intermediate: > 33% and ≤ 66%. High: > 66%

It is clear, however, that the DNA methylation marks at these loci are not completely stable. At Hyper1 and Hyper2, for example, there were F1 s with very little methylation, suggesting that there must have been loss of methylation even from the IL8–3 allele. It is likely therefore that there are two oppositely acting mechanisms influencing the methylation status of these loci: a paramutation-like event leading to de novo establishment of DNA methylation and a second process leading to removal of DNA methylation. The net effect of these processes is likely to account for the extent to which the pattern of DNA methylation deviates from Mendelian ratios in these F1 and F2 plants. However once established in the F1, methylation may be fairly stable: F2 progeny from two out of three F1s were consistently highly methylated for Hyper1 and Hyper2.

In addition to these four loci (Hyper1–4, Fig. 7a) we also investigated inheritance at five other loci that were hypermethylated in the IL1–1, IL2–5 and IL8–3 (Additional file 7: Figure S3 and Table 1). At Hyper5 and Hyper6 there was no conclusive evidence for paramutation-like behaviour (Fig. 7a) and we conclude either that these are conventional Mendelian loci or that the instability of the methylation mark in the progeny of M82 × IL8–3 offset the effects of the paramutation-like mechanism. At Hyper7 and Hyper8, high methylation in some M82 plants as assessed by McrBC does not allow us to conclude that the lack of F2s with low methylation reflects non-Mendelian segregation (Additional file 7: Figure S3). Methylation in IL8–3 at Hyper9 was variable, and there was no evidence for paramutation-like interactions (Additional file 7: Figure S3).

Hypomethylated DMRs

We also tested inheritance of 8 hypomethylated DMRs in IL1–1, IL2–5 and IL8–3 (Hypo1–8, Fig. 7b and Additional file 7: Figure S3). For one locus (Hypo1), there was a high proportion of the hypermethylated allele in the F2 (from two F1s in particular) although in the F1 the trend was towards medium-to-low methylation. These proportions indicate that methylation of this locus is subject to instability but that the hypermethylated M82 epiallele has some paramutation-like activity.

Similarly at Hypo2 and Hypo3, there was evidence for contrasting dynamics of methylation: 4 and 5 (out of 10) F1s were highly methylated, while 3 in each case had very low methylation (Fig. 7b). Overall the proportions of lowly methylated F2s did not significantly depart from Mendelian ratios, but the F1 of origin had clear effects, with progenies of some F1s homogeneously hypermethylated. These results argue for paramutation-like activity as well as instability of the hypermethylated M82 epialleles.

At Hypo4 and Hypo5 there was no evidence for non-Mendelian segregation (Fig. 7b), while for the remaining three loci (Hyper6–8, Additional file 7: Figure S3), McrBC estimates of methylation in the parents were too variable to interpret the methylation levels of their progeny.

In summary, of the 17 DMRs tested in detail, 7 showed at least partial evidence of a heritable gain of methylation upon crossing (Table 1). Because the methylation differences were observed in three independent introgression lines, it is unlikely that they are due to a trans effect from the introgressed S. pennellii region and a more likely explanation is that there is an epigenetic effect that may be partially related to paramutation-like mechanisms when methylation is gained in the M82 × IL8–3 F1s and F2s.


Origins of DNA methylation differences: Paramutation-like mechanisms versus spontaneous epigenetic variation

The main aim of this study was to explore the epigenetic changes associated with wide cross hybridisation of M82 × S. pennellii through the characterisation of H06 and other DMRs in at least three different ILs relative to the parental lines. In principle these shared DMRs, as with DMRs in hybrid progeny of maize [18], rice [19] and Arabidopsis thaliana [12, 20], could be due to spontaneous epimutation [21, 22] or genetic differences in seedstocks (e.g. a transposon insertion) [23, 24] as well as to genetic or epigenetic interactions between the parental genomes. For many of the DMRs shared by the three ILs under study, we cannot say which of these explanations is likely to apply. For H06, however, we could reproduce its epimutation whenever we made the M82 × S. pennellii cross and paramutation in backcrosses of IL8–3 with M82 occurred at at high frequency, suggesting that the epigenetic state of the ILs at H06 is due to the hybridisation process. To explain hybrid-induced paramutation it could be that there are genetic or epigenetic differences between the parental genomes. A genetic mechanism might apply if the two genomes differ by insertion or deletion mutation at the affected locus. In such a scenario it could be, by analogy with epigenetic marks in Neurospora crassa that the transition is triggered by unpaired DNA during meiosis of the F1 [25]. The difference between the plant and fungal systems in the properties of the epigenetic mark: in N. crassa it would be a standard heritable mark with Mendelian inheritance whereas in tomato it would be a heritable mark that can transmit between alleles.

This requirement for unpaired DNA would account for the delayed onset of the paramutation-like process of H06 until the meiosis in the reproductive phase of the F1 (Fig. 3). An alternative epigenetic explanation could invoke sRNAs from one of the interacting alleles that, in the parental line, are at too low an abundance to trigger paramutation. If the second allele has characteristics that favour production of secondary sRNA, for example repeats [5], then the formation of the hybrid would trigger a high level of sRNA and RdDM. The well established feedback loops of the RdDM process [26] would ensure that, once established at one allele, the amplified RdDM is stable and has the epigenetic drive characteristic of paramutation in subsequent generations.

These hypotheses are not mutually exclusive because both RdDM and meiotic silencing of unpaired DNA are dependent on an RNA-dependent RNA polymerase [27, 25]. The involvement of RdDM is supported by several lines of evidence that link small RNAs and in particular 24-nt sRNAs with paramutation-like events in tomato in addition to the findings from other species [11, 5, 12, 13]. At H06, for example, the timing of the paramutation-like process correlates with the timing or sRNA production (Fig. 3) and in VIGS the viral sRNAs produced could mediate (although inefficiently) the stable methylation of H06 and the onset of paramutation-like properties (Fig. 4). In addition to H06, out of four loci where the hypermethylated epiallele was associated with upregulated sRNAs, three showed evidence of paramutation-like behaviour (Hyper1, Hyper2 and Hypo3, Fig. 7 and Table 1).

There are, however, epialleles with paramutation-like characteristics (Hyper3, Hyper4, Hypo1 and Hypo2, Fig. 7) that are not associated with differential small RNAs. It remains possible that RNA-independent mechanisms are involved although some of these examples may be due to sRNA production and gain of methylation in tissues or developmental stages that are not well represented in the samples used for this analysis. Conversely, not all epialleles with differential sRNA have distinct epigenetic marks or, even if they do, are paramutagenic. Clearly paramutation requires additional, as yet unknown, factors beyond those associated with heritable epigenetic marks.

Stability of epialleles and penetrance of paramutation-like effects

For one locus with two epialleles (unmethylated and methylated), there are four possible states in a non-mosaic diploid organism: ‘unmethylated/unmethylated’, ‘unmethylated/methylated’, ‘methylated/unmethylated’, ‘methylated/methylated’ (Additional file 8: Figure S4). Local positive feedback loops in the maintenance of methylation [26, 27] or non-methylation [28] contribute to the stability of each state, while spontaneous epimutations and epiallele interactions in epi-heterozygotes trigger state transitions and distortions in the expected Mendelian ratios (Additional file 8: Figure S4).

Epialleles associated with paramutation in maize can be highly stable as with the b1 locus, or unstable and revert to the active state, as with the pl1 locus [29]. Similarly, with the paramutation-like loci in tomato there are varying degrees of stability. For some loci there were F1 plants without methylation (H06, Hyper1, Hyper2, Hypo1, Hypo2, Hypo3, Fig. 7), indicating a degree of instability of the methylated epiallele. Currently there are no clear mechanisms to explain the loss of methylation in an epi-heterozygote, apart from a dilution by a factor two of the sRNAs that could weaken a RdDM positive feedback loop acting in trans. The penetrance of the paramutation-like interaction (transition to the ‘methylated/methylated’ state) also varied: while 9/10 F1s had high methylation at Hyper1, only 5/10 were highly methylated at Hyper2 (Fig. 7). It is therefore crucial to garner information about stability and paramutagenicity of a large number of epialleles in order to identify their determinants. However we note that a binary depiction of the epigenetic state (e.g. methylated/unmethylated) is not complete: while H06IL is very stably methylated, the paramutated H06IL’ can revert when backcrossed to M82 (Fig. 2).

Frequency of paramutation in tomato hybrids

Paramutation with perfect efficiency will be fixed very rapidly in a population and will therefore be difficult to detect. Using a wide cross such as Solanum lycopersicum cv. M82 × Solanum pennellii LA0716 increases the likelihood of detecting highly penetrant paramutation events. We based our genome-wide search for paramutation on the characteristics of H06 paramutation, which is associated with a gain of methylation in all sequence contexts and abundant sRNAs in multiple ILs at the seedling stage. We only found one other locus, Hyper1, that shared all of these characteristics, and this locus also displayed paramutation-like properties. However segregation analyses (Fig. 7) revealed that loci without differential sRNAs or without hypermethylation in all contexts could also engage in paramutation-like behaviours. Unexpectedly, loci that were hypomethylated in the ILs could also gain methylation upon crossing to M82.

These findings suggest that paramutation-like interactions are in fact common; occurring at hundreds of loci in a M82 × IL8–3 cross. However fully penetrant paramutation is rare: primary paramutation at H06 was the most penetrant, but it displayed weak secondary paramutation (Fig. 2); and the penetrance of epigenetic changes and their inheritance in the seven validated loci (from seventeen tested) was incomplete. Therefore we conclude that full penetrance is rare and the IL crossing scheme may be too stringent (three consecutive backcrosses) to identify many traces of paramutation from the initial M82 × S. pennellii cross. Further studies aiming at assessing the prevalence of paramutation will need to be powered to detect incomplete penetrance.

Effects on gene expression

DNA methylation is primarily associated with the silencing of repeated sequences. Consistently earlier examples of paramutation are often linked with transposable elements [5]. Thus, as a mechanism of identity-based silencing, paramutation is thought to contribute to silencing multiple copies of transposable elements throughout the genome, and differential methylation may only occasionally also affect gene expression. We found no evidence for differential expression of the genes surrounding H06 in the seedling expression data set. However H06 is a target of the RIN (RIPENING INHIBITOR) transcription factor [30] and its surrounding genes are most highly expressed in the ripening fruit [31] (Additional file 9: Figure S5), so its differential methylation may have an effect in this particular tissue. For three of the validated DMRs that engaged in paramutation-like interactions (Hyper3, Hyper4 and Hypo3), the neighbouring genes were differentially expressed in the introgression lines compared to M82 (Additional file 6: Figure S2B).

Generating novel epialleles can be a source of phenotypic diversity [32, 33]. If the epialleles are stable enough and their epigenetic state can be efficiently controlled, epigenetic drives via paramutation may be used to improve crops.


In this study we show that interactions between natural tomato epialleles are relatively common, can lead to heritable differences in methylation, and occur mostly in one direction (gain of methylation) while stochastic loss of methylation may act as a counterbalance. There were however variations between pairs of epialleles in timing, penetrance and heritability of the gain of methylation. We propose that interactions between epigenomes should be viewed as a continuum, of which fully penetrant and fully heritable paramutation is an extreme manifestation. Taking into account the partial penetrance and stochasticity of epigenetic interactions will be crucial to properly investigate transgenerational inheritance.


Plant material

Tomato (Solanum lycopersicum) cv. M82, S. lycopersicum Micro-Tom, Solanum pimpinellifolium and Solanum pennellii plants were raised from seeds in compost (Levington M3) and maintained in a growth room at 23 °C with 16/8 h light/dark periods with 60% relative humidity, at a light intensity of 150 μmol photons m− 2.s− 1. Unless otherwise indicated, S. pennellii refers to accession LA0716 that was used to generate the introgression lines [14].

DNA extraction

DNA from leaf (2-week-old seedlings), mature flowers or pollen was extracted with the Puregene kit (QIAgen) following manufacturer’s instructions.


IL8–3 was genotyped based on the TG510 marker: the sequence was amplified with primers TG510 fw and TG510 reverse, and digested with AluI. The M82 sequence is cleaved in a 133-bp and a 245-bp fragment, while the introgressed sequence from S. pennellii remains intact (378 bp). The H06 sequence from M82 was differentiated from S. pennellii and S. pimpinellifolium by digesting H06 amplicons with HpaI. The S. pennellii and S. pimpinellifolium are digested into two fragments, 260 and 210 bp in length, while the M82 sequence remains uncut. Even brightness of the H06 genotyping bands for F1 s on an agarose gel indicate that McrBC-qPCR results should not be distorted by differences in DNA amplification efficiency between alleles.

RNA extraction

Total RNA samples were prepared from 100 mg of tissue using TRIzol (LifeTechnologies).


H06 genomic insert was cloned into the binary TRV RNA2 vector using the KpnI and XhoI restriction sites of the multiple cloning site as described previously [16, 34]. Cotyledons of tomato seedlings were agro-infiltrated 10 d after sowing with a 1:1mixture of Agrobacterium tumefaciens (strain GV3101:pMP90 + pSOUP) carrying TRV RNA1 and RNA2 at OD600 = 1.5. 11 plants were succesfully infected with TRV-H06 and 6 with wild-type TRV.

Viral load quantification

Quantification of viral load was performed with Precision One-Step qRT-PCR (Primerdesign) on 15 ng RNA per well and normalised with TIP41.


sRNAs from leaf (2-week-old seedlings), flower and pollen of M82 (two individual plants), S. pennellii LA0716 (two plants) and their F1 (four plants) were cloned from 10 μg total RNA using the Illumina TruSeq Small RNA cloning kit and libraries were indexed during the PCR step (12 cycles) according to the manufacturer’s protocol. Gel size-selected, pooled libraries were sequenced on a HiSeq 2000 50SE. Padubiri Shivaprasad prepared Micro-Tom (one plant) and S. pimpinellifolium (one plant) sRNA libraries from the aerial part of two-week-old seedlings according to [15]. Sequences were trimmed and filtered with Trim Galore! (with the adapter parameter ‘-a TGGAATTCTCGGGTGCCAAGG’). Mapping to H06 (S. pennellii and M82 sequences) was performed with Bowtie 1.1.1 [35] without mismatches (options ‘-v 0 -a’).

sRNA counts on H06 in Fig. 1 were obtained by mapping sRNA libraries from [15] to Heinz genome SL2.50 with Bowtie 1.1.1 without mismatches (options ‘-v 0 -a’), and extracting the counts for the interval “SL2.50ch08:54487325–54,487,842”. We found that there was a mistake in the labelling of the libraries deposited in [15]: the libraries labelled IL1–1, IL2–5, IL8–1-1, IL8–1-D, IL8–1-3, IL8–2, IL8–2-1, IL8–3 and IL8–3-1 actually correspond to IL8–3-1, IL8–3, IL8–2-1, ILH8–2, IL8–1-5, IL8–1- D, IL8–1-1, IL2–5 and IL1–1. For convenience the correctly labelled sRNA libraries for the introgression lines used in this study (IL1–1, IL2–5 and IL8–3) are included in the GEO data set.

To find differential sRNA loci between M82 and ILs, sRNA libraries from IL1–1, IL2–5, IL8–3 and two M82 seedlings (from [15]) were mapped and clustered on Heinz genome SL2.50 using ShortStack v3.3.3 [36] with default parameters. sRNA counts on the defined loci were analysed with DESeq2 v1.8.1 [37].

Gene expression analysis

Transcript counts for two-week-old IL1–1, IL2–5, IL8–3 and M82 seedlings from [17] were analysed with DESeq2 [37] with the design factors ‘condition’ (sun or shade), ‘experimental batch’, and ‘genotype’ (either IL or M82). In each of the five experimental replications, there were 4 M82 and 1 IL1–1, 1 IL2–5, and 1 IL8–3 in each treatment (sun or shade), totalling 40 M82 and 30 IL libraries on which to run the differential expression analysis. Genes were considered differentially expressed between genotypes when the adjusted p-value was for the ‘genotype’ factor was < 0.05.

Gene ontology enrichment analysis

Differentially expressed gene IDs were tested for gene ontology term enrichment against the full set of S. lycopersicum genes with the PANTHER Overrepresentation Test (Released 20,171,205) (, GO Ontology database Released 2017–12-27, Fisher’s exact test with FDR multiple test correction).

Analysis of DNA methylation


Analysis of methylation by McrBC was performed as previously described in [16].

Sanger bisulfite

For Sanger bisulfite sequencing, 450 ng of DNA was bisulfite-converted with EZ DNA Methylation-Gold Kit (Zymo Research) and amplified with primers specific to the region and the Kapa Uracil+ HotStart DNA polymerase (Kapa Bioscience). Amplification products were size selected on a 1.5% agarose gel and gel extracted with the QIAquick gel extraction kit (QIAgen), A-tailed following the same protocol as for the library preparation (below), and cloned into pGEM-T easy (Promega). Sequences aligned with MUSCLE were then analysed with CyMATE [38] and PCR-duplicates were removed. Each figure represents at least nine independent clones.


Bisulfite library preparation was performed in duplicates for each genotype with a custom protocol similar to ref. [39]. 1.2 μg DNA was sonicated on a Covaris E220 to a target size of 400 bp and purified on XP beads (Ampure, ratio 1.8). DNA was end-repaired and A-tailed using T4 DNA polymerase and Klenow Fragment (NEB) and purified again using XP beads (ratio 1.8×). Methylated Illumina Y-shaped adapters for paired-end sequencing were ligated using Quick-Stick Ligase (Bioline). 450 ng of purified (ratio 1.8×), adapter-ligated DNA was bisulfite-converted using EZ DNA Methylation-Gold Kit (Zymo Research) according to manufacturer’s instructions. DNA was barcoded using 12 cycles of PCR amplification with KAPA HiFi HotStart Uracil+ Ready Mix (Kapabiosystems) with PE1.0 and custom index primers (courtesy of the Sanger Institute). Duplicate libraries for M82, IL1–1, IL2–5 and IL8–3 (leaves of two-week-old seedlings) were sequenced in paired-end mode. Reads were trimmed and filtered with Trim Galore! v0.4.2 (default parameters), then mapped on Heinz genome SL2.50 using Bismark v0.17.0 [40] (first in paired-end mode with options ‘–score-min L,0,-0.2 -p 4 –reorder –ignore quals –no-mixed –no-discordant –unmapped’, then unmapped read1 was mapped in single-end mode with the same quality parameter ‘-N 1’). Reads were deduplicated with ‘bismark deduplicate’ and methylation calls were extracted using Bismark ‘methylation extractor’ (with options ‘-r2 2’ for paired-end reads). Based on Bismark’s cytosine report, methylated and unmethylated counts for cytosines of both strands were pooled into 300-bp bins sliding by 200 bp and separated by context (CG, CHG, and CHH). Bins with fewer than 10 counts or with coverage exceeding the 99 th percentile in one or more libraries were excluded from the analysis. We also excluded the introgressed regions: ch01 1–86 Mb, ch02 44–54 Mb and ch08 61 Mb – 65,866,657 bp (end of chromosome). Each introgression line was then compared to M82 by fitting a logistic regression on the methylated and unmethylated counts (R glm with ‘family = binomial’). The p-values were subjected to a Benjamini-Hochberg correction to control the false discovery rate at 5%. Differentially methylated bins were further filtered by imposing thresholds on the absolute difference in methylation (average of replicates) between M82 and the introgression line: a difference of 25% in the CG context, 20% in CHG and 10% in CHH. DMRs were constructed by merging overlapping differentially methylated bins. We included the region Hyper2 into the McrBC validation analysis despite it being inside the introgressed region in IL2–5, because it was hypermethylated in CG, CHG and CHH contexts and produced abundant sRNAs in IL1–1 and IL8–3, thus making it a likely candidate for paramutation. Chloroplast DNA was used as a control for bisulfite conversion efficiency, and sequencing statistics are collected in Additional file 1.


Please refer to Additional file 10.



Differentially methylated region


Introgression line


RNA-directed DNA methylation


Small RNA


Virus-induced gene silencing


  1. Chandler VL, Stam M. Chromatin conversations: mechanisms and implications of paramutation. Nat Rev Genet. 2004;5(7):532–44.

    Article  CAS  PubMed  Google Scholar 

  2. Sapetschnig A, Sarkies P, Lehrbach NJ, Miska EA. Tertiary siRNAs Mediate Paramutation in C. elegans. PLoS Genetics. 2015;11(3):1005078.

    Article  Google Scholar 

  3. de Vanssay A, Bougé A-L, Boivin A, Hermant C, Teysset L, Delmarre V, Antoniewski C, Ronsseray S. Paramutation in Drosophila linked to emergence of a piRNA-producing locus. Nature. 2012;490(7418):112–5.

    Article  CAS  PubMed  Google Scholar 

  4. Rassoulzadegan M, Grandjean V, Gounon P, Vincent S, Gillot I, Cuzin F. RNA-mediated non-mendelian inheritance of an epigenetic change in the mouse. Nature. 2006;441(7092):469–74.

    Article  CAS  PubMed  Google Scholar 

  5. Hollick JB. Paramutation: a trans-homolog interaction affecting heritable gene regulation. Curr Opin Plant Biol. 2012;15(5):536–43.

    Article  CAS  PubMed  Google Scholar 

  6. Hollick JB. Rmr6 maintains meiotic inheritance of Paramutant states in Zea mays. Genetics. 2005;171(2):725–40.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Dorweiler JE, Carey CC, Kubo KM, Hollick JB, Kermicle JL, Chandler VL. Mediator of paramutation1 is required for establishment and maintenance of paramutation at multiple maize loci. Plant Cell. 2000;12(11):2101–18.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Stonaker JL, Lim JP, Erhard KF, Hollick JB. Diversity of pol IV function is defined by mutations at the maize rmr7 locus. PLoS Genet. 2009;5(11):1000706.

    Article  Google Scholar 

  9. Sidorenko L, Dorweiler JE, Cigan AM, Arteaga-Vazquez M, Vyas M, Kermicle J, Jurcin D, Brzeski J, Cai Y, Chandler VL. A dominant mutation in mediator of paramutation2, one of three second-largest subunits of a plant-specific RNA polymerase, disrupts multiple siRNA silencing processes. PLoS Genet. 2009;5(11):1000725.

    Article  Google Scholar 

  10. Hagemann R. Somatische Konversion bei Lycopersicon esculentum mill. Zeitschrift für Vererbungslehre. 1958;89(4):587–613.

    CAS  PubMed  Google Scholar 

  11. Gouil Q, Novák O, Baulcombe DC. SLTAB2 is the paramutated SULFUREA locus in tomato. J Exp Bot. 2016;67(9):2655–64.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Zhang Q, Wang D, Lang Z, He L, Yang L, Zeng L, Li Y, Zhao C, Huang H, Zhang H, Zhang H, Zhu J-K. Methylation interactions in Arabidopsis hybrids require RNA-directed DNA methylation and are influenced by genetic variation. Proc Natl Acad Sci. 2016;113(29):4248–56.

    Article  Google Scholar 

  13. Greaves IK, Eichten SR, Groszmann M, Wang A, Ying H, Peacock WJ, Dennis ES. Twenty-four-nucleotide siRNAs produce heritable trans-chromosomal methylation in F1 Arabidopsis hybrids. Proc Natl Acad Sci. 2016;113(44):6895–902.

    Article  Google Scholar 

  14. Eshed Y, Zamir D. An introgression line population of Lycopersicon pennellii in the cultivated tomato enables the identification and fine mapping of yield- associated QTL. Genetics. 1995;141(3):1147–62.

    CAS  PubMed  PubMed Central  Google Scholar 

  15. Shivaprasad PV, Dunn RM, Santos BA, Bassett A, Baulcombe DC. Extraordinary transgressive phenotypes of hybrid tomato are influenced by epigenetics and small silencing RNAs. EMBO J. 2011;31(2):257–66.

    Article  PubMed  PubMed Central  Google Scholar 

  16. Bond DM, Baulcombe DC. Epigenetic transitions leading to heritable, RNA-mediated de novo silencing in Arabidopsis thaliana. Proc Natl Acad Sci. 2015;112(3):917–22.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Chitwood DH, Kumar R, Headland LR, Ranjan A, Covington MF, Ichihashi Y, Fulop D, Jimenez-Gomez JM, Peng J, Maloof JN, Sinha NR. A quantitative genetic basis for leaf morphology in a set of precisely defined tomato introgression lines. Plant Cell. 2013;25(7):2465–81.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Regulski M, Lu Z, Kendall J, Donoghue MTA, Reinders J, Llaca V, Deschamps S, Smith A, Levy D, McCombie WR, Tingey S, Rafalski A, Hicks J, Ware D, Martienssen RA. The maize methylome influences mRNA splice sites and reveals widespread paramutation-like switches guided by small RNA. Genome Res. 2013;23(10):1651–62.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Chodavarapu RK, Feng S, Ding B, Simon SA, Lopez D, Jia Y, Wang G-L, Meyers BC, Jacobsen SE, Pellegrini M. Transcriptome and methylome interactions in rice hybrids. Proc Natl Acad Sci. 2012;109(30):12040–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Greaves IK, Groszmann M, Ying H, Taylor JM, Peacock WJ, Dennis ES. Trans chromosomal methylation in Arabidopsis hybrids. Proc Natl Acad Sci. 2012;109(9):3570–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Becker C, Hagmann J, Müller J, Koenig D, Stegle O, Borgwardt K, Weigel D. Spontaneous epigenetic variation in the Arabidopsis thaliana methylome. Nature. 2011;480(7376):245–9.

    Article  CAS  PubMed  Google Scholar 

  22. Schmitz RJ, Schultz MD, Lewsey MG, O’Malley RC, Urich MA, Libiger O, Schork NJ, Ecker JR. Transgenerational epigenetic instability is a source of novel methylation variants. Science. 2011;334(6054):369–73.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Eichten SR, Briskine R, Song J, Li Q, Swanson-Wagner R, Hermanson PJ, Waters AJ, Starr E, West PT, Tiffin P, Myers CL, Vaughn MW, Springer NM. Epigenetic and genetic influences on DNA methylation variation in maize populations. Plant Cell. 2013;25(8):2783–97.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Stuart T, Eichten SR, Cahn J, Karpievitch YV, Borevitz JO, Lister R. Population scale mapping of transposable element diversity reveals links to gene regulation and epigenomic variation. elife. 2016;5:2–15.

    Article  Google Scholar 

  25. Shiu PKT, Raju NB, Zickler D, Metzenberg RL. Meiotic silencing by unpaired DNA. Cell. 2001;107(7):905–16.

    Article  CAS  PubMed  Google Scholar 

  26. Matzke MA, Mosher RA. RNA-directed DNA methylation: an epigenetic pathway of increasing complexity. Nat Rev Genet. 2014;15(6):394–408.

    Article  CAS  PubMed  Google Scholar 

  27. Law JA, Jacobsen SE. Establishing, maintaining and modifying DNA methylation patterns in plants and animals. Nat Rev Genet. 2010;11(3):204–20.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Inagaki S, Miura-Kamio A, Nakamura Y, Lu F, Cui X, Cao X, Kimura H, Saze H, Kakutani T. Autocatalytic differentiation of epigenetic modifications within the Arabidopsis genome. EMBO J. 2010;29(20):3496–506.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Stam M, Scheid OM. Paramutation: an encounter leaving a lasting impression. Trends Plant Sci. 2005;10(6):283–90.

    Article  CAS  PubMed  Google Scholar 

  30. Zhong S, Fei Z, Chen Y-R, Zheng Y, Huang M, Vrebalov J, McQuinn R, Gapper N, Liu B, Xiang J, Shao Y, Giovannoni JJ. Single-base resolution methylomes of tomato fruit development reveal epigenome modifications associated with ripening. Nat Biotechnol. 2013;31(2):154–9.

    Article  CAS  PubMed  Google Scholar 

  31. The Tomato Genome Consortium. The tomato genome sequence provides insights into fleshy fruit evolution. Nature. 2012;485(7400):635–41.

    Article  Google Scholar 

  32. Cortijo S, Wardenaar R, Colomé-Tatché M, Gilly A, Etcheverry M, Labadie K, Caillieux E, Hospital F, Aury J-M, Wincker P, Roudier F, Jansen RC, Colot V, Johannes F. Mapping the epigenetic basis of complex traits. Science. 2014;343(6175):1145–8.

    Article  CAS  PubMed  Google Scholar 

  33. Yang X, Kundariya H, Xu Y-Z, Sandhu A, Yu J, Hutton SF, Zhang M, Mackenzie Sa. MutS HOMOLOG1-derived epigenetic breeding potential in tomato. Plant Physiol. 2015;168(1):222–32.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Liu YL, Schiff M, Dinesh-Kumar SP. Virus-induced gene silencing in tomato. Plant J. 2002;31(6):777–86.

    Article  CAS  PubMed  Google Scholar 

  35. Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009;10(3):25.

    Article  Google Scholar 

  36. Axtell MJ. ShortStack: comprehensive annotation and quantification of small RNA genes. RNA. 2013;19(6):740–51.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.

    Article  PubMed  PubMed Central  Google Scholar 

  38. Hetzl J, Foerster AM, Raidl G, Mittelsten Scheid O. CyMATE: a new tool for methylation analysis of plant genomic DNA after bisulphite sequencing. Plant J. 2007;51(3):526–36.

    Article  CAS  PubMed  Google Scholar 

  39. Urich MA, Nery JR, Lister R, Schmitz RJ, Ecker JR. MethylC-seq library preparation for base-resolution whole-genome bisulfite sequencing. Nat Protoc. 2015;10(3):475–83.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Krueger F, Andrews SR. Bismark: a flexible aligner and methylation caller for bisulfite-Seq applications. Bioinformatics. 2011;27(11):1571–2.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Winter D, Vinegar B, Nahal H, Ammar R, Wilson GV, Provart NJ. An “electronic fluorescent pictograph” browser for exploring and analyzing large-scale biological data sets. PLoS One. 2007;2(8):718.

    Article  Google Scholar 

Download references


The authors would like to thank Shuoya Tang and Mel Steer for horticultural assistance, Aleksandar Ivanov for his work on sRNA libraries, Padubiri Shivaprasad for the Micro-Tom and S. pimpinellifolium sRNA libraries, Bruno Santos, Felix Krueger, and Sebastian Müller for bioinformatics advice.


This work was supported by the European Research Council Advanced Investigator Grant ERC-2013-AdG 340642, the Balzan Foundation, the Biotechnology and Biological Sciences Research Council and the Frank Smart Studentship. D.C.B. is the Royal Society Edward Penley Abraham Research Professor. The funding bodies had no roles in the design of the study, in the collection, analysis, interpretation of data, or in the writing of the manuscript.

Availability of data and materials

The datasets generated and analysed during the current study are available in the GEO repository under accession GSE97247 (

Author information

Authors and Affiliations



QG and DCB designed the research and wrote the manuscript, QG performed the research. All authors read and approved the final version of the manuscript.

Corresponding authors

Correspondence to Quentin Gouil or David C. Baulcombe.

Ethics declarations

Ethics approval and consent to participate

S. lycopersicum cv. M82, S. lycopersicum cv. Micro-Tom, S. pimpinellifolium, and S. pennellii LA0716 seeds were kindly provided by Tamas Dalmay. Seeds for the introgression lines were kindly provided by Dani Zamir. All plants were grown, propagated, crossed and collected in growth rooms at the Plant Growth Facility of the Department of Plant Sciences. Plant work was done under licence from the Department of Food and Rural Affairs (DEFRA) n° 50,973/198463 and VIGS experiments were conducted under DEFRA Plant Health licence 50,973/235589.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1:

Summary of MethylC-seq data. Tsv format. (TXT 358 bytes)

Additional file 2:

DMRs shared between introgression lines. Tsv format. (TXT 453 kb)

Additional file 3:

DMRs shared between introgression lines and hypermethylated in all cytosine contexts. Average proportion of methylated cytosines between two replicates, for each context and genotype. Tsv format. (TXT 4 kb)

Additional file 4:

DMRs shared between introgression lines and hypomethylated in all cytosine contexts. Average proportion of methylated cytosines between two replicates, for each context and genotype. Tsv format. (TXT 1 kb)

Additional file 5:

Figure S1. Genomic position, small RNAs and DNA methylation of the IL DMRs Hyper1–9 and Hypo1–8. The selected DMR is highlighted in red, the plotted region includes 2 kb upstream and downstream. Genes: ITAG2.4 gene models. Repeats: RepeatMasker annotation. Hypomethylated regions shared between the three introgression lines are annotated in blue, and hypermethylated regions in red. sRNAs: coverage of sRNAs in seedlings. Methylation: percentage of methylated cytosines in 100-bp regions in each context (two replicates per genotype). (PDF 853 kb)

Additional file 6:

Figure S2. A) Distribution of zero-count data for the genes differentially expressed between M82 and ILs (70 samples: 40 M82 and 30 ILs, data from [17]). B) Differential expression of genes whose promoters overlap shared IL DMRs selected for segregation analysis. Hypermethylation of Solyc09g064640’s promoter (Hyper3 and Hyper5) is associated with decreased expression in introgression lines. Solyc09g064640 codes for a protein related to the family with sequence similarity 65 member A (FAM65A). Solyc04g025030, encoding a MADS-box transcription factor similar to Agamous-like MADS-box protein AGL62, is associated with Hyper4, and Solyc06g075840, encoding a putative zinc-binding protein of the Yippee family, with Hypo3. Normalised counts in seedlings (40 M82 and 30 ILs, data from [17]). (PDF 64 kb)

Additional file 7:

Figure S3. Segregation of DNA methylation patterns by McrBC. These regions show variable methylation in the parental controls M82 and IL8–3. (A) Hypermethylated DMRs Hyper7–9. (B) Hypomethylated DMRs Hypo6–8. For each region, the left hand panel shows the results of methylation analysis by McrBC for individual plants. This information is collated in the right hand panel with the splitting of the F2 s according to their F1 parent. Low methylation: ≤ 33%. Intermediate: > 33% and ≤ 66%. High: > 66%. (PDF 135 kb)

Additional file 8:

Figure S4. Schematic of possible epiallelic transitions for one locus in a diploid state with unmethylated (white) and methylated (black) epialleles. There are four possible epigenetic configurations. Recurrent (curved) arrows represent the propensity to conserve the current state, from the action of local positive feedback loops maintaining methylation/unmethylation. Transition (straight) arrows between states represent spontaneous epimutations and, when starting from epigenetically heterozygous states, paramutation-like interactions. The epigenetic state of a genome may be seen as the result of such a stochastic process at each cell/organism generation. (PDF 2 kb)

Additional file 9:

Figure S5. Potential role of H06 in fruit ripening. Expression profiles across tissues of the two genes immediately downstream of H06, Solyc08g066020, encoding a Serine C-palmitoyltransferase like protein (A) and Solyc08g066030, encoding an unknown protein (B). Their expression peaks in mature green fruit. Data from [31] visualised on the eFP Browser [41] at Copyright permissions for the re-use of the tomato drawings kindly granted by Prof. Provart. (C) Demethylation of the 3′ end of H06 during fruit ripening correlates with the binding of the RIN transcription factor at this locus. RIN binding, demethylation and ripening are compromised in the rin (ripening inhibitor) and cnr (colorless non-ripening) mutants. Data from [30]. (PDF 971 kb)

Additional file 10:

Oligonucleotides used in this study. Tsv format. (TXT 1 kb)

Additional file 11:

H06 sequences for S. lycopersicum cv. M82, S. pimpinellifolium and S. pennellii LA0716. Fasta format. (FASTA 1 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Gouil, Q., Baulcombe, D.C. Paramutation-like features of multiple natural epialleles in tomato. BMC Genomics 19, 203 (2018).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: