A pure-line rice cultivar, Jinongda18 (Ma et al., 2003) was used in this study. Uniform germinating seeds were pretreated with 355, 532, and 1,064 nm Nd³⁺ YAG pulsed laser for 6, 10, 15, 20 times, respectively (Table 1). Eleven treatment groups with different cumulative doses of laser irradiation were set, and the untreated wild-type (WT) plants were used as the mock control. More than 1000 plants were used for the mock and 11 stress treatments, and the treated seeds were allowed to continue germinating in petri-dishes for 5d in a growth cabinet (30°C during the day and 25°C during the night, 16/8 h photoperiod at 50 μmolm⁻²s⁻¹). Seedlings were then transferred to a homogeneous experimental paddy-field plot at Jilin Agricultural University, Changchun, China, in accord with season. Plant showing prominent phenotypic changes especially in plant height and heading date were tagged. The heading date was recorded as the first plant showing panicle emergence. Five plants were randomly chosen to measure their height. A single treated M0 individual with the shortest height and the latest heading date, designated as M0#11, was selected to generate the M1 progeny by self-pollination, then, a single M1 individual, designated as M1#1, was chosen to produce the M2 generation by self-pollination.

The flag-leaves of rice plants in WT group and the 11 low-dose laser-treated groups were sampled as pools, while the M0#11, M1, and M2 progeny were sampled as individuals for DNA isolation at the grain filling stage. DNA was isolated using the modified CTAB method and purified by phenol extractions. MSAP is a modified version of the amplified fragment length polymorphism (AFLP) to detect the stability and alteration in cytosine DNA methylation at the 5′-CCGG sites (Yaish et al., 2014). Genomic DNA was digested with EcoRI combined with HpaII or MspI (New England Biolabs, Beverly, Massachusetts), ligated with EcoRI and H/M adapters, and then amplified with one pair of pre-selective and 20 pairs of selective primers (Supplementary Table 1). The amplification products of MSAP were resolved by 5% denaturing polyacrylamide gel electrophoresis and visualized by silver staining. Only clear and completely reproducible bands in two independent experiments were scored.

MSAP is performed using HpaII and MspI, a pair of isoschizomers that recognize the same restriction site (5′-CCGG) but have different sensitivities to methylation of the cytosines. HpaII will not cut if either of the cytosines is fully (double-strand) methylated, whereas MspI will not cut if the external cytosine is fully- or hemi- (single-strand) methylated (McClelland et al., 1994). Thus, for a given DNA sample, H0M1 indicates CG methylation; H1M0 indicates CHG methylation; H1M1 indicates no methylation; and H0M0 indicates CG/CHG methylation. Based on the principles, the changes of MSAP patterns were divided into four categories: (a) CG hyper: H1M1 to H0M1, H1M0 to H0M0; (b) CHG hyper: H1M1 to H1M0, H0M1 to H0M0; (c) CG hypo: H0M1 to H1M1, H0M0 to H1M0; and (d) CHG hypo: H1M0 to H1M1, H0M0 to H0M1.

Bands showing alteration in methylation patterns in the MSAP profiles in the laser-irradiated M0 plants and/or their progenies relative to WT were isolated, boiled with ddH₂O, and amplified with the pre-selective primers used in the original MASP analysis. The PCR products were purified with Wizard PCR Preps Purification System (Promega). Then the purified PCR products were ligated into pM18-T vector (Promega), and then sequenced. Homology analysis was performed by BlastX at the NCBI website (http://www.ncbi.nlm.nih.gov/). Sequence alignment was done by the CLUSTALW program using the Genedoc software.

Total RNA was isolated from the flag-leaves at the same developmental stage as that used for DNA methylation analysis by the Trizol Reagent (Invitrogen). Then RNA was treated with DNaseI (Invitrogen), reverse-transcribed by the SuperScript TMRNase H-Reverse Transcriptase (Invitrogen), and subjected to qRT-PCR analysis using gene-specific primers (Roche LightCycler 480). Genes encoding DNA methytransferase1 or MET1 (LOC_Os07g08500, DMT707), CMT3 (LOC_Os05g13790, DMT703), DRM2-1 (LOC_Os03g02010, DMT706), and DRM2-2 (LOC_Os05g04330, DMT710), DME1 (LOC_Os01g011900, DNG702), DME2 (LOC_Os05g37350, DNG701), DDM1 (LOC_Os03g51230, CHR741), Argonaute AGO1-1 (LOC_Os02g45070, AG0711), AGO1-2 (LOC_Os04g47870, AG0708), AGO4-1 (LOC_Os04g06770, AG0705), AGO4-2 (LOC_Os01g16870, AG0703) were analyzed in WT, M0, M1, and M2 progeny to interrogate the impacts of laser irradiation on their staedy-state transcript abundance. These gene-specific primers were designed by the Primer Premier 5 software (Supplementary Table 2). A β-actin gene of rice (LOC_Os05g0438800) was used as an internal control. DNA contamination was tested by inclusion of RNAs without RT. Three batches of independently isolated RNAs were used as biological replications. The relative amounts of the gene transcripts were determined using the Ct (threshold cycle) method and the fold-change data were analyzed by the 2^−ΔΔCt method.

Based on the principle that presence vs. absence of mPing produces 433 bp difference in length in given locus, excision of mPing from a known locus can be easily detected by locus-specific PCR amplification. Thus, a set of 53 pairs of locus-specific primers was designed by the Primer Premier 5 software (Supplementary Table 3), and each pair of primers bracketed an intact mPing based on the rice reference genome sequence of the standard laboratory genotype Nipponbare (http://rgp.dna.affrc.go.jp). The amplicons were visualized by ethidium bromide staining after electrophoresis through 2% agarose gels. A set of empty donor sites for mPing excision were identified, isolated and sequenced, together with their corresponding mPing-containing loci.

Transposon display (TD) was performed by combining the MseI-adaptor-specific primers either with a mPing subterminal-specific primer or with a set of inter-simple sequence repeat (ISSR) exactly as previously reported (Jiang et al., 2003). Major procedures, i.e., digestion, ligation, pre-selective amplication, selective amplication were the same as for MSAP. And the restriction enzyme MseI (New England Biolabs) was used, the pre-selective primers were MseI +0 (5′- GATGAGTCCTGAGTAA) and mPing internal amplification primer (5′- GCTGACGAGTTTCACCAGGATG), the selective primers were 8 MseI +2 (Supplementary Table 4) and mPing proximal ends amplification primers (5′- TGTGCATGACACACCAGTG). The novel bands in TD gels of laser irradiated plants (M0, M1, and M2) relative to WT were considered as putative mPing de novo insertions.

Genomic DNA (3 μg per lane) of the various plants was digested by HindIII or XbaI (New England Biolabs). Digested DNA was run through 1% agarose gel and transferred onto Hybond-N⁺ nylon membranes (Amersham Pharmacia Biotech, Piscataway, New Jersey) by the alkaline transfer recommended by the supplier. And the mPing (positions: 6–430), Pong (BK000586.1), and Ping (AB087616.1) probes were prepared as previously described (Lin et al., 2006). The probe-fragments were gel-purified and labeled with fluorescein-11-dUTP by the Gene Images random prime-labeling module (Amersham Pharmacia Biotech). Hybridization signal was detected by the Gene Images CDP-Star detection module (Amersham Pharmacia Biotech). The filters were exposed to X-ray films.

Genomic DNA was modified using an EZ DNA Methylation-Gold kit (Zymo Research) according to the manufacturer's recommendations. Modified DNA was purified using a Zymo-Spin IC column (Zymo Research). The primers of retrotransposon Tos17, an established locus in the rice genome for assessing DNA methylation liability, for bisulfite sequencing were designed using the MethPrimer program (http://www.urogene.org/methprimer/) (Supplementary Table 5). For each PCR, 1.0–3.0 μl of bisulfite-treated DNA was used, and the PCR products were cloned into the pMD18-T vector and sequenced. More than 20 clones were sequenced for each sample. Analyses of the bisulfate sequencing results were conducted at the Kismeth website (http://katahdin.mssm.edu/kismeth). The methylation levels per site (CG, CHG, and asymmetric CHH) were calculated by dividing the number of non-converted (methylated) cytosines by the total number of cytosines within the assay (Ngezahayo et al., 2009).

Statistical significance was determined using SPSS 11.5 for Windows (http://www.spss.com/statistics/). Continuous variables e.g., plants height, fold-change data in qRT-PCR were presented as the mean ± standard deviation and statistically tested by an unpaired, two-tailed t-test. A value of P < 0.05 was considered significant and P < 0.01 was very significant. The Pearson correlation analysis between the DNA methylation pattern variations (detected by MASP) and the expression levels of methylation-related genes (detected by qRT-PCR) was calculated by the range method using SPSS. −1 ≦ Pearson's r ≦ 1. Pearson's r > 0 indicates the two variables have a positive correlation, Pearson's r < 0 indicates the two variables have a negative correlation, Pearson's r = 0 indicates the two variables have no correlation. The color depth symbolizes the values of Pearson's r, and the red and blue blocks denote positive and negative correlations, respectively.

Phenotyping was conducted at the reproduction stage after the laser treatment in rice cultivar Jinongda 18. Phenotypic variations especially plant height and heading dates were observed in all the 11 laser treatment groups with variable extents (Table 1). Of 1,000 laser irradiated germinating seeds, ca. 990 plants grew to maturity in the experimental paddy-field. The heading date was recorded as the first plant showing panicle emergence in one group, and it was 84 days in WT, and which was delayed in most laser-treated groups. For example, heading date was 88 days in groups 7 and 10, and 89 days in group 11 (Figure 1). These observations indicated that long-wave laser treatment was more effective than short-wave laser, and the treatment frequency was a key factor which influences the degree of phenotypic variations in the same wavelength.

Five plants were randomly chosen to quantitatively measure the plant height, and the mean plant height was 101 ± 1.58 cm in WT (Table 1), and the common feature was dwarfing in all the 11 laser-treated groups. For example, it was 95 ± 2.12 cm in group 7 (P < 0.01), 95 ± 1.58 cm in group 10 (P < 0.01), the shortest height was 93 ± 3.00 in group 11 (P < 0.01). It indicated that laser irradiation negatively influenced plant growth and development, and the degrees of impact were largely dose-dependent.

A total of 1,140 clear and reproducible bands (between two technical replicates, starting from the first step, i.e., DNA isolation) were scored using 20 pairs of selective primers in the MSAP analysis (Supplementary Table 1) for the 11 selected M0 individuals that showed clear phenotypic variations, along with WT (Figure 2A). The tabulated methylation level was 21.2% (^mCG + ^mCHG) in WT, and all but one (M0-1) of the 11 M0 plants showed similar levels of methylation to that of WT, which collectively ranged from 19.4 to 21.1% (Figure 2B). This suggests that the laser irradiation treatments did not cause large changes in total DNA methylation level at the 5′-CCGG sites. Notably, however, for the single M0 plant (M0-1) that did show clear loss of methylation, it mainly occurred at the CG sites (Figure 2B). In contrast to the general stability in methylation level, we found that all four methylation patterns (CG hyper, CG hypo, CHG hyper, and CHG hypo) that are discernible based on the MSAP data (Materials and methods) showed significant changes in all 11 laser irradiated M0 plants relative to WT. The frequencies of these methylation pattern alterations ranged from 1 to > 4% (Figure 2C). The similar frequencies of hyper vs. hypo methylation changes at both CG and CHG sites clearly have offset each other and render the total methylation levels broadly constant (Figure 2B). Notably, in all the MSAP profiles we did not find parallel band changes between the two isoschizomers, HpaII and MspI, indicating nucleotide sequence changes at the 5′-CCGG sites should be rare if occurred at all in these low-dose laser treated rice plants (Figure 2A).

Next, we asked whether the methylation changes were inheritable to the next organismal generation. For this purpose, we chose plant M0-11 which showed the most obvious phenotypic changes (Table 1) to produce the M1 generation via self-pollination. A total of 543 reproducible bands (between two technical replicates) were scored using 10 informative selective MSAP primers in 18 randomly chosen M1 plants (Figure 3A). Results showed that the total methylation level of M0#11 was 20.3%, and which ranged from 19.7 to 22.8% in the M1 plants (Figure 3B), suggesting a moderate fluctuation in methylation level occurred in the M1 plants, which contrasted with the broad methylation level stability among the M0 plants (Figure 2B). Analysis of the four methylation patterns revealed the important observation that a great majority of the altered methylation patterns in the M0 plant (M0-11) was inherited to all 18 plants of the M1 generation (Figure 3A). Nevertheless, it is important to note that heritability of the alterations was clearly much less than 100% and variable across the 18 M1 individuals (Figure 3C). Moreover, gaining of additional methylation pattern alterations, especially CG hypo and both CG hyper and CHG hyper patterns, were apparent in most of the M1 plants (Figure 3C). The less than 100% heritability, together with unbalanced gain vs. loss of additional alteration, is consistent with the discernible fluctuations in total methylation level in the M1 plants (Figure 3B).

The still dynamic nature of methylation patterns in the M1 plants promoted us to address the issue further. We thus chose a single M1 plant (M1-1) to produce the M2-generation plants by self-pollination. A total of 543 reproducible bands (between two technical replicates) were scored by using seven pairs of selective MSAP primers across 16 randomly selected M2 individuals (Figure 4A). Results indicated that the methylation levels of all 16 M2 plants remained similar to that of the M1 plant (Figure 4B). Importantly, except for one M2 plant (M2-2), the two hyper-methylation patterns (CG hyper and CHG hyper) were uniform and the same as those of M0-11 and M1-1 (Figure 4C), pointing to their stabilization. By contrast, the two hypomethylation patterns (CG hypo and CHG hypo) were still as fluctuating among the M2 plants (Figure 4C) as they were among the M1 plants (Figure 3C). Together, our results suggest that although some of the newly acquired methylation pattern alterations due to the treatment were inherited and stabilized by the M2 generation, DNA methylation patterns in general are still dynamic in the M2 generation.

Tos17 is an endogenous long terminal repeat (LTR) retrotransposon in the rice genome, and which is usually heavily methylated over it entire length under normal conditions. The methylation state of Tos17 however is highly labile and loss of methylated cytosines (^mCs) may occur under various conditions, which may (but not necessarily) lead to its transpostional activation (Liu et al., 2004). Thus, Tos17 has been widely accepted as a sensitive mark to assess methylation liability in rice (Liu et al., 2004; Ding et al., 2007). To verify that laser irradiation has caused changes in DNA methylation patterns and their transgenerational heritability by an independent approach, we analyzed methylation state of the 5′ portion of Tos17 together with its immediate flank (Table 2) by locus-specific bisulfite sequencing in several M1 and M2 individuals that showed clear methylation changes in the MSAP profiles along with their WT. First, for the entire sequenced segment, there are two CG sites located in the flank region (Table 2), which are slightly methylated (ca. 10%) in WT, but the methylation was completely lost in all analyzed progeny individuals derived from the low-dose laser-treated mother plant (Figure 5A). There are nine CHG sites in the sequenced segment, four in the flank and five in the 5′ portion of Tos17 (Table 2). In addition, there are 50 CHH sites in the sequenced segment, 27 in the flank and 23 in the body region of Tos17 (Table 2). Analysis of these sequences found that (i) in WT, CG, CHG, and CHH sites were low-methylated, while only a small fraction of CHG sites had a high frequency of methylation. (ii) In the M1 plant, CG sites were slightly or moderately methylated compared with WT, most CHG sites were maintained the methylation patterns of M0, while CHH sites showed inheritance of M0 or reversal of the methylation patterns to those of WT. (iii) In the M2 plant, nearly all three kinds of methylation at CG, CHG and CHH were identical with those of M1 (Figure 5), denoting largely stable epigenetic inheritance thereafter.

Establishment, erasure and maintenance of DNA methylation patterns are accomplished by the interlaced action of multiple enzymes encoded by a large set of genes collected termed chromatin genes (McCue et al., 2015; Rigal et al., 2016; Tan et al., 2016). To test whether the laser irradiation-induced heritable changes of DNA methylation patterns might be related to altered expression of the chromatin genes, especially those directly related to DNA methylation, we analyzed expression levels of a set of genes encoding DNA methyltransferases (MET1-1, CMT3-1, DRM2-1, and DRM2-2), 5-methylcytosine DNA glycosylases (DME1, DME2), the SWI/SNF chromatin remodeler (DDM1), and four siRNA pathway-related AGO proteins (AGO1-1, AGO1-2, AGO4-1, AGO4-2) in all the 11 M0 plants (Table 1). Real-time PCR results indicated that all the analyzed genes showed significantly perturbed expression in these M0 plants (Supplementary Table 6). Specifically, genes coding for DNA methyltransferases CMT3-1, DRM2-1, and DRM2-2 were all significantly down-regulated in the M0 plants. Similarly, genes coding for 5-methylcytosine DNA glycosylases DME1 and DME2 were also down-regulated in the M0 plants, and so were genes coding for the siRNA pathway-related AGO proteins AGO1-1, AGO1-2, AGO4-1, and AGO4-2 (Supplementary Table 6). In contrast, the gene coding for SWI/SNF chromatin remodeler DDM1 was significantly up-regulated in all the M0 plants (Supplementary Table 6).

To test whether the significantly perturbed expression states of these chromatin genes in the M0 plants were also transgenerationally persistent, all the 18 M1 plants and 16 M2 plants were subjected to the same qRT-PCR analysis (Supplementary Table 7). Results indicated that DDM1, which was significantly up-regulated in M0, displayed a further up-regulation in virtually all the M1 plants, whereas MET1-1, CMT3-1, DRM2-1, DME1, DME2, AGO1-1, AGO1-2, AGO4-1, AGO4-2, which were down-regulated in M0 (Supplementary Table 6), showed up-regulated expression in nearly all the M1 plants (Supplementary Table 7). DRM2-2, which was also down-regulated in M0 (Supplementary Table 6), showed largely reversion to the expression state of WT (Supplementary Table 7). Notably, DRM2-1 and DRM2-2 showed a dramatic upregulation in two (M1-11 and M1-12) of the M1 plants (Supplementary Table 7). Further analysis of the M2 generation plants indicated that CMT3, DRM2-1, DME1, and AGO1-1 showed stable inheritance of the expression states of their M1 mother plant (Figure 6A), while MET1, DRM2-2, DME2, AGO4-1, and AGO4-2 showed largely reversion of expression states to those of the M0 plant, which nevertheless were still significantly different from those of WT (Figure 6B). DDM1 which was significantly up-regulated in the M0 and M1 plants, showed a tendency of down-regulation in all but two (M2-5 and M2-6) of the studied M2 plants, and in these two plants it showed a highly outlier upregulation (Supplementary Table 8).

In addition, A Pearson correlation analysis was used to analyze the possible correlations between the DNA methylation patterns (detected by MSAP) and the perturbed expression of chromatin genes (detected by qRT-PCR). The value of Pearson's r symbolizes the correlation tendency between the two variables. We found that, in M0 plants, CG hypermethylation was positively correlated with expression AGO4-2 (r = 0.71, P < 0.05) (Figure 7A). In M1 plants, CG hypomethylation was positively correlated with expression of AGO4-1 (r = 0.60, P < 0.01), AGO4-2 (r = 0.69, P < 0.01) and MET1-1 (r = 0.51, P < 0.05), while CHG hypermethylation was negatively correlated with expression AGO4-1 (r = −0.52, P < 0.05) and AGO4-2 (r = −0.55, P < 0.01) (Figure 7B). However, no correlation was detected in M2 plants, consistent with the MSAP profiles showing that the DNA methylation patterns were stabilized by large in the M2 generation (Figure 4C). Together, the above data of transgenerational expression dynamics suggested that the immediate changes, inheritance, additional changes, reversion to parental state as well as stabilization of the altered CG and CHG methylation patterns in the laser irradiated plants and their offspring are closely linked to the transgenerationally perturbed expression states of at least some of the chromatin genes analyzed.

It should be cautioned, however, that the expression dynamics of these chromatin-related genes both among the laser-irradiated individuals and between generations clearly did not follow a persistent trend. Therefore, we cannot rule out the possibility that stochastic factors like expression noise may have also contributed to the expression dynamics.

To get some insights into the functional relevance of the genomic loci underwent methylation pattern changes as a result of laser-irradiation, 36 reproducible variant bands (shared among at least three M0 individuals) were eluted from the MSAP gel profile, cloned and sequenced. A BlastX query indicated that 16 loci were homologous to known or predicted protein-coding genes, three loci were TEs, and the other 15 loci showed no similarity with known sequences in the database (Supplementary Table 9), suggesting that at least certain variant methylation loci are likely functionally consequential.

mPing (430 bp in length) is the most active MITE-type TE in the rice genome, and can be induced to transpose under various stress conditions, and its mobility can be associated with DNA methylation of itself and/or its flanks (Ngezahayo et al., 2009). Given the extensive DNA methylation pattern changes induced by the laser irradiation, we sought to test the possibility that mPing had been mobilized by the treatment. For this purpose, 53 mPing-containing loci identified in the WT were analyzed by locus-specific PCR, and the results showed that mPing excision occurred in 21 (ca. 40%) loci in one or more of the 11 selected laser-irradiated M0 plants (Table 3), and mPing mobility was most active in M0#11 (Figure 8A). It has been documented previously that excision of mPing is not necessarily accompanied by reinsertions of the mobilized copies (Shan et al., 2005). To test if mPing reinsertions occurred in these plants, transposon-display (TD) was performed. Although exhaustive TD was conducted, only six reinsertion events in the M0#11 plant or some of its 19 M1 progenies (bands appeared de novo relative to WT), and all six loci were isolated and confirmed by newly designed PCR primers (e.g., Figure 8B). This suggests that although excisions of original mPing copies occurred abundantly (ca. 40%, 21 out of 53 loci), most mobilized copies failed to reinsert back to the rice genome. To confirm this possibility further, Southern blotting was conducted using the full-length mPing as a probe. Indeed, only loss of bands were detected in M0#11 and its progenies at the resolution of Southern blotting (Figure 8C), confirming the possibility, which is consistent with a previous study (Shan et al., 2005).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.