Atlas of Age- and Tissue-Specific DNA Methylation during Early Development of Barley (Hordeum vulgare)

2.1 Plant material and growth conditions

Five spring barley varieties (Barque 73, Flagship, Hindmarsh, Schooner and Yarra), were selected based on their similarity in phenology in order to minimize epigenetic variability between varieties associated with developmental differences. Seeds from all varieties were provided by the Salt Focus Group at the Australian Centre for Plant Functional Genomics (ACPFG, Adelaide, South Australia), and planted at the same time in potting mix comprising 50% UC (University of California at Davis), 35% coco-peat and 15% clay/loam (v v⁻¹) in 3.3 L pots, 17.5 cm deep, free-draining and placed on saucers. The experiment was conducted from 30th January to 20th February 2015 in a greenhouse at the Waite Campus, University of Adelaide, South Australia (34°58′11″S, 138°38′19″E). The seedlings were grown under natural photoperiod, while temperatures were set at 22°C/15°C (day/night). The experiment consisted of five randomized blocks of five varieties (25 seedlings per block). Pots were watered to weight every 2 days to a gravimetric water content of 16.8% (w w⁻¹) (0.8 × field capacity) [26] until sampling 21 days after sowing, when seedlings were at three-leaf stage (Zadok stage 13 [27]). Blades and sheaths of leaves 1–3 were sampled separately. Leaves 1 and 2 were fully expanded prior to sampling, whilst leaf 3 had just completed growth. About 50 mg of plant material was cut from the middle section of each leaf blade and each leaf sheath and collected in 2 ml micro tubes. Roots were cut from the seedlings and washed using tap water to remove soil particles, then blotted dry with paper towels before sampling 50 mg of root tissue. All samples were snap frozen in liquid nitrogen, and then stored at −80°C until DNA extraction. In total, 175 tissue samples were collected, including 25 root samples (i.e. 5 plants per each of the five varieties used in the study), 75 leaf blade samples (i.e. from leaves 1, 2 and 3 from each of the 5 plants per variety used in the study) and 75 leaf sheath samples (i.e. from leaves 1, 2 and 3 from each of the 5 plants per variety used in the study).

2.2 DNA isolation

Prior to DNA extraction, frozen plant material was homogenized in a bead beater (2010-Geno/Grinder, SPEX SamplePrep^®, USA). DNA isolation was performed from pulverised plant samples using a Qiagen DNeasy kit and following the manufacturer’s instructions. DNA samples were quantified using a NanoDrop^® 1000 Spectrophotometer (V 3.8.1, ThermoFisher Scientific Inc., Australia) and concentrations were standardized to 10 ng/μl for subsequent library preparation.

2.3 Methylation-sensitive genotyping by sequencing (ms-GBS)

The ms-GBS was performed using a modified version [28] of the original GBS technique [56]. Genomic DNA was digested using the combination of a methylation-insensitive rare cutter, EcoRI (GAATTC), and a frequent and methylation-sensitive cutter, MspI (CCGG). Each sample of DNA was digested in a reaction volume of 20 μl containing 2 μl of NEB Smartcut buffer, 8 U of HF-EcoRI (High-Fidelity) and 8 U of MspI (New England BioLabs, Australia). The reaction was performed in a BioRad 100 thermocycler at 37°C for 2 h, followed by enzyme inactivation at 65°C for 10 min.

Then the ligation of adapters to individual samples was achieved in the same plates by adding 0.1 pmol of the respective barcoded adapters with an MspI cut site overhang, 15 pmol of the common Y adapter with an EcoRI cut site overhang, 200 U of T4 Ligase and T4 Ligase buffer (New England BioLabs, Australia) in a total volume of 40 μl. Ligation was carried out at 24°C for 2 h followed by an enzyme inactivation step at 65°C for 10 min.

DNA samples were allocated to plates, 81 samples each, including the negative control, water. Prior to pooling plate samples into a single 81-plex library, the ligation products were individually cleaned up to remove excess adapters using an Agencourt AMPure XP purification system (Beckman Coulter, Australia) at a ratio of 0.85 (AMPure magnetic beads/ligation product), following the manufacturer’s instructions. Individual GBS libraries were produced by pooling 25 ng of DNA from each sample. Each constructed library was then amplified in eight separate PCR (25 μl each) containing 10 μl of library DNA, 5 μl of 5× Q5 high fidelity buffer, 0.25 μl polymerase Q5 high fidelity, 1 μl each of Forward and Reverse common primers at 10 μM, 0.5 μl of 10 μM dNTP and 7.25 μl of sterile pure water. PCR amplification was performed in a BioRad T100 thermocycler, consisting of DNA denaturation at 98°C (30 s) and 10 cycles of 98°C (30 s), 62°C (20 s) and 72°C (30 s), followed by 72°C for 5 min. PCR products were next pooled to reconstitute libraries. DNA fragments between 200 and 350 bp in size were captured using AMPure XP magnetic beads following the manufacturer’s instructions. Bead-captured fragments were eluted in 35 μl of water, of which 30 μl were collected in a new labelled microtube. Libraries were next paired-end sequenced in an Illumina HiSeq 2500 (Illumina Inc., USA) at the Australian Genome Research Facility (AGRF, Melbourne Node, Australia). Sequencing results were deposited in the European Nucleotide Archive (ENA) (Study Accession Number: PRJEB27251).

2.4 Analysis of global differences in DNA methylation between samples

Differences in ms-GBS profiles between samples were explored by performing principal component-linear discriminant analysis (PC-LDA) (a supervised clustering approach for high dimensional data), using different levels of hierarchy between samples as the putative drivers in DNA methylation differences (i.e. grouping samples by organ (root vs. leaf), tissue (root vs. blade vs. sheath, and tissue) and age (root vs. leaf 1 vs. leaf 2 vs. leaf 3 vs. sheath 1 vs. sheath 2 vs. sheath 3)). PC-LDA was implemented using the R package FIEmspro 1.1-0 [29] on the standardized coverage, the count per million reads (CPM) of the 913,697 ms-GBS markers generated. PC-LDA results were visualized by a scatter plot of the first two discriminant factors (DFs), and a 3D plot using the first three DFs. Finally we used an unsupervised hierarchical cluster analysis to generate a dissimilarity tree based on Mahalanobis distance [30] generated also based on the standardized coverage (CPM) of the 913,697 ms-GBS markers.

2.5 Detection of DMMs in barley

Differentially methylated DNA was assessed in ^mCCGG motifs (recognised by MspI), between barley leaf parts (blade and sheath) and roots. To do so, samples were grouped according to organ type (root, blade and sheath) regardless of the genotype of origin, making 25 samples per organ. This approach aimed to minimise genotype-dependent methylation markers. DMMs were identified using the package, msgbsR, developed by Mayne et al. [31]. DMMs were selected based on FDR adjusted P-values with a threshold of 0.05 [32, 33]. The significance of the marker also fulfilled the condition that the read counts reached at least 1 CPM and was present in at least 20 samples per organ type (maximum sample per group = 25). The logFC (logarithm 2 of fold-change) was computed to estimate the intensity and directionality of differential DNA methylation between tissues. Determining the directionality of DNA methylation uses the fold change as an inverse proxy for change in the methylation level. That is, higher methylation levels on a specific locus will reduce the number of MspI restriction products and therefore reduce the number of sequences generated for that locus [34].

2.6 Distribution of DMMs around genomic features

To test whether there was a relationship between tissue-specific DMMs and particular genomic features (e.g., genes and repeat regions as defined in Ensembl database (http://plants.ensembl.org/biomart/martview/)), DMM distribution was assessed in the barley genome. Therefore, DMMs stable between tissues were mapped to the barley reference genome. Then, the number of DMMs within genomic features (repeats, genes, exons, UTRs and tRNA genes) and per 1 kb bins within 5 kb flanking regions [24, 28] was tallied using the shell module, bedtools/2.22.0 [35].

3.1 Methylation-sensitive genotyping by sequencing

The sequencing of the 170 samples of barley tissue which met DNA quality requirements yielded over 900 million raw reads, with more than 91% bases above Q30 (99.9% accuracy of base call [36]) across all samples (Table 1). Of these reads, 99.27% contained the barcode and EcoRI/MspI adapters ligated during library construction. Further filtering was performed to retain reads that strictly aligned with the barley reference genome. In this way, we obtained nearly 450 million reads (50.10%), with a mean of 2,637,916 high quality reads per sample. These high-quality reads accounted for 913,697 sequence tags, representing 32.30% of the 2,828,642 CCGG sites in the barley genome. Of these sequence tags, 748,594 (80.62%) showed some form of polymorphism for methylation between samples.

3.2 Estimation of tissue- and tissue rank-dependent epigenetic differentiation

The PC-LDA plots revealed clear evidence of structuring of methylation between samples (Figure 1a). A 3D plot using the first three discriminant factors (DF1, DF2 and DF3) revealed that blades and sheaths were further grouped according to the rank of the leaf from which they were harvested. The distance between blades and sheaths seems to shrink with leaf rank (Figure 1b). This leaf rank-dependent grouping was also supported by hierarchical cluster analysis (HCA) of the distances between sample group centres (Figure 1c), based on the Mahalanobis distance [29, 30], and sample clusters matched the leaf developmental age (Figure 1c). Leaf rank-dependent DNA methylation differences were further assessed between tissues by comparing the methylation profiles of blades and sheaths for each rank of leaf appearance. No DMMs were observed between the three leaf blades, whereas sheaths 1 and 3 presented 18 DMMs (Table 2).

3.3 Differentially methylated DNA markers between roots and leaves

DMMs between barley roots and leaves were obtained through comparison of the read count per million of tissue types, independently of genotypes. This comparison revealed substantial DMMs between both roots vs. blades and roots vs. sheaths (Figure 2a), and there were more DMMs between roots and blades (6510 DMMs Figure 2b) than between roots and sheaths (4116 DMMs Figure 2c). Of these markers, 3266 DMMs were present in both blades and sheaths when compared to roots, and their methylation changed consistently in the same direction in each comparison (Figure 3a). The number of DMMs between roots and leaf blades increased with leaf-rank, whereas DMMs between roots and leaf sheaths did not show any relationship with rank (Figure 2a). Tissue-specific DMMs were predominantly hypomethylated (95–98%) in leaf parts (sheath or blade) compared to roots (Figure 2a). This result was in line with the median of the fold-changes of DMMs, which indicated an overall DNA hypomethylation in leaves (Figure 4a and b). From here on, DMMs consistently present in roots vs. sheaths and roots vs. blades will be designated as stable markers between roots and leaves.

3.4 Differentially methylated DNA markers between the leaf blade and sheath

There was only a small number of DMMs between leaf blades and sheaths (0–73 DMMs, Table 2 and Figure 2d). These DMMs were basically between leaf blades and sheaths 1 and 2; and there was none between blade 1 and sheath 3. There was only 1 DMM between sheath 3 and blades 2 and 3 (Table 2 and Figure 2d). Pairwise comparisons between blades 1–2 and sheaths 1–2 revealed 20 common DMMs, which were all hypermethylated in sheaths compared to blades (Figures 2e and 4b). Half of the 20 common DMMs between blades and sheaths were located on chromosome 5H. Furthermore, there were no DMMs in pairwise comparisons among blades 1–3 and among sheaths 1–3, except between sheath 1 and sheath 3 which had 18 DMMs (Table 2). However, comparing blades and sheaths of the same leaf rank showed 32 DMMs between blade 1 and sheath 1, 36 DMMs between blade 2 and sheath 2 and 1 DMM between blade 3 and sheath 3.

3.5 Distribution of tissue-specific DMMs around genes

Relatively few of the tissue-specific DMMs were located around gene exons. Indeed, of the 3266 stable DMMs between root and leaf samples, only 60 (1.8%) were located within 5 kb of a gene, including 21 overlaps with genes and 39 DMMs that were spread within 5 kb upstream and downstream of genes (Figure 5a). Apart from the absence of DMMs within 1 kb upstream of transcription start sites, there was no obvious tissue-specific DMM distribution pattern around the genes (Figure 5a). The same assessment process showed that, as with common DMMs, only a small proportion of blade-specific DMMs (44 of 3246, 1.3%) was positioned close to a gene (Figure 5b). Of these, 15 DMMs overlapped with a gene transcript, whereas the remaining 29 DMMs were distributed within 5 kb of the gene without any clear pattern (Figure 5b), except that the number of DMMs located between 2 and 3 kb bins was higher both upstream and downstream, than any other 1 kb bin within the 5 kb flanking regions (Figure 5b). There were fewer sheath-specific methylation markers within 5 kb from genes than blade-specific markers (13 of 2391 DMMs, 0.5%) (Figure 5c). The majority of these (10 out of 13 DMMs) were sited within 3 kb of a gene, and no DMMs were present 3–5 kb from transcription margins (Figure 5c). Of 37 gene-body DMMs detected across all comparisons (Figure 5a–c), 27 overlapped with an exon and the remaining 10 markers were in intronic regions, 70–604 bp upstream of exons, except 1 DMM, which was 62 bp downstream an exon (Appendix A).

3.6 Distribution of tissue-specific DMMs near repeat regions

Many more tissue-specific DMMs were detected near repeats than near genes. The DMMs around repeat regions (as defined in the Ensembl database (http://plants.ensembl.org/biomart/martview/)) were concentrated either within the repeats or within 1 kb of their margins (Figure 6a). A similar distribution pattern was obtained with both blade-specific and sheath-specific DMMs when contrasted with roots, with more DMMs overlapping with the repeats themselves than in the 1 kb stretches flanking their margins (Figure 6b and c). The few markers that were differentially methylated between blades and sheaths (20 DMMs in total) were all located within 1 kb of a repeat (Figure 6d). Therefore, stable tissue-specific DMMs appeared to occur preferentially within repeats and 1 kb flanking regions, with higher frequency within 1 kb downstream than within 1 kb upstream, regardless of tissue types (Figure 6a–d).

3.7 Distribution of genes around differentially methylated (DM) repeats

To investigate a possible interaction between differentially methylated (DM) repeats and genes, the distance of genes from DM repeats between root and leaf samples was evaluated. In this way, we found 105 genes near repeats (up to 5 kb either side), of which 37 overlapped with a repeat and the remaining genes were scattered up- and downstream from the repeat (Figure 7). The number of DM repeats surrounded by genes thus represented only a tiny proportion of the total repeats that were differentially methylated between roots and leaves (105 out of 3266 DM repeats, 3.21%). About half of genes near DM repeats (52 of 105 genes) were also differentially methylated, whereas the remainder (53 genes) were not.

4.1 Extensive epigenetic differentiation between roots and leaves

In this study, we detected large numbers of DDMs between roots and leaves that were conserved across a diverse array of barley genotypes, and so were deemed far more likely to be organ-specific than genotype-dependent. Of these, hypomethylation of the ^mCCGG motif predominated in leaves (Figures 2b and c, 3b and 4a). More surprisingly, we also detected similarly conserved DMMs between leaf-blades and leaf-sheaths (Figures 2e and 4b). The number of conserved DMMs between blades and sheaths (20 DMMs), all hypermethylated in sheaths, was relatively consistent with the closeness of these structures in position and function. These findings are broadly congruent with previous studies, which reported differential DNA methylation between variable tissues (e.g. endosperm, pollen, leaves, and roots) in diverse plant species [7, 8, 9, 10], but additionally hint that the developmental closeness of structures being compared may also be reflected in the distinctiveness of their methylation profiles. However, controversy over the extent and validity of organ-specific DMMs [9, 10, 21, 22, 23] could cast doubt over their utility for organ diagnosis or as a tool to gain greater insight into the genes responsible for organ development/identity. Here, we sought to mitigate against the possibility of type I errors in DMM assignment through the unprecedented use of five diverse varieties and five biological replicates of each variety in the identification of these marks. In contrast to our findings, previous workers have reported little difference in the methylation levels of both ^mCG and ^mCHG motifs between roots and leaves in Arabidopsis [9] and sorghum [10]. Further, no significant difference was detected at all for ^mCG and ^mCHG methylation levels between tissues in cotton [37]. These divergences may simply reflect genuine biological differences between taxonomic groups. However, it is also important to recognise that such differences may also arise from the approach used to identify organ-specific DMMs. Variability in the techniques used to assess plant methylation profiles may introduce different forms of bias and preclude or complicate comparison among studies. DMM detection can be influenced by factors such as (1) the genome coverage of the methylation profiling method (low coverage methods such as MSAP are likely to miss many markers) [7], and (2) the data analysis approach used, which can compare either global methylation levels (e.g. percent methylation) [9] or methylated loci (e.g. DMMs) [28]. We contend that relying solely on global methylation levels can be misleading in comparing tissue profiles, because similar methylation levels may show completely different patterns and so vital information content is lost.

The current study revealed that tissue-specific DNA methylation occurred abundantly in the ^mCHG context (in particular ^mCCGGs) (Figure 2a and c). This concurs with reports of the CHG context similarly dominating differential DNA methylation between organs in Brachypodium distachyon [8] and sorghum [10]. Although tissue-specific methylation also occurs in other cytosine contexts [10], our results and other studies [10, 22] suggest that ^mCCGG is a primary motif of epigenetic distinctiveness of plant organs. However, since MspI activity is affected by the presence of cytosine hydroxymethylation on its recognition sequence [38], some of the markers identified here as being cytosine methylation induced, could be due to the presence of (de)hydroxymethylation events instead. Additionally, while tissue-specific DMMs were mostly hypomethylated in leaves compared to roots in the present study (Figure 3b), in Arabidopsis, Widman et al. [9] found that hypermethylation prevailed in leaves compared with roots. This apparent contradiction in the directionality of methylation in DMMs between roots and leaves may be a reflection a difference in the polarity of early divisions in the monocotyledonous barley and the dicotyledonous Arabidopsis embryos or else the methylation profiling method implemented.

4.2 DNA methylation flux is tissue-specific during barley seedlings development

In addition to tissue-specificity of methylation profiles, one notable finding in the current study was that leaf cohorts exhibited a strong tendency to co-cluster. This suggests that the nature of methylation divergence between organs is not absolutely fixed and instead appears to change with developmental progression. This observation accords with previous reports that genome-wide methylation patterns are not static during plant development [39]. Additionally, a considerable portion of DMMs between roots and leaves was also specific to the leaf rank, due to the steady decrease in the number of DMMs between roots and leaf blades with the rank of the latter (Figure 2a–c). In this case, therefore, the slow but progressive accumulation of additional methylation marks in the leaves increases their divergence from root profiles and enables the separation of leaf cohorts. However, the small number of DMMs distinguishing between leaf blades and leaf sheaths ran counter to this trend such that there were no DMMs capable of discrimination between these leaf parts among the oldest cohort studied (leaf 1) (Figure 2d and Table 2). It seems intuitively improbable that older cohorts of leaves would simply lose differentiation between structurally distinct parts, especially if these marks had a functional role in defining function. Perhaps the most plausible biological explanation for the apparent erosion of divergence lies in the different chronological ages of the leaf cohorts that were sampled. Put simply, the third leaves were the least mature of the three cohorts collected and so it is entirely possible that the blade-sheath differential marks had yet to appear in these samples. Thus, it is important to consider the developmental and ageing progression chronology when assigning DMMs and that some organ- or structure-specific marks may only become organ-specific late in their development. Such late-emerging developmental DMMs should mean that the cumulative number of tissue-specific markers increases and so the organs or structures become more distinct, through leaf growth stages [40], each of which may carry a specific epigenetic profile. Certainly, others have noted that methylation profiles vary progressively as the organ develops [3, 41, 42] before reaching, at maturity, a “default” methylome, which may be conserved across varieties [24]. These results suggest that, once leaves are differentiated and mature, they do not show significant differences in DNA methylation profiles, regardless of their rank of appearance. Additionally, the location of half of the 20 common DMMs between blades and sheaths on chromosome 5H implies that this chromosome carries loci important for blade and sheath identities.

4.3 Tissue-specific DNA methylation preferably occurs in repeat regions of the barley genome

Organ-specific DMMs identified here were primarily associated with repeat regions. No significant difference was observed between the frequency of CCGG sites in and around genes and repeats. However, 84% of the barley genome is comprised of mobile elements or other repeat structures [25, 43], indicating that the fact that the majority of detected DMMs are located within or in the proximity of a repeat is due to the intrinsic repetitive nature of the studied genome. Nevertheless, the fact that 27 DMMs overlapped with exons and 10 were located in introns (Appendix A) contradicts previous claims that CHG methylation marks are exclusively restricted to repeat regions and intergenic regions [20, 21, 44, 45]. The possible regulatory significance of such gene body CHG methylation marks requires further investigation [46]. However, it is already well-established that tissue-specific DMMs can influence gene expression by enhancing gene transcription [9] and alternative splicing [47] or through repression due to immediate proximity to transcription start site [48].

The predominance of DMMs around and within repeats leads us to speculate that they could play an important role in defining organ identity in barley, and accords with previous findings in Brachypodium distachyon [8]. This flux of DNA methylation patterns in repeats [8, 42, 49] has been proposed to regulate [44] developmental shifts during plant growth and development [11, 39]. Nevertheless, the association between DMMs in/around repeat regions and organ identity described here does not establish a causal link between the two. However, there are grounds for reasoning that this may be the case and that the possibility warrants further study. First, repeat regions were previously proposed to be involved in alternative promoters, a substantial proportion of which (>40%) was reported to shape tissue differentiation [16]. Therefore, tissue-specific DMMs in repeats may contribute to alternative promoters, and thus determine organ identity. Second, differential gene expression between roots and leaves [25, 50] implies a firm regulatory system, including epigenetic mechanisms to guarantee tissue-specific cell development. Tissue-specific DMMs in repeats show that repeats are not the so-called “selfish parasites” of the genome [51], but can directly or indirectly affect tissue-specific gene expression [42, 52, 53]. Finally, it has been suggested that transposons coordinate splice variants, a genomic event that occurs in more than 60% of plant genes [54, 55], thus generating multiple mRNA transcripts from a single gene [56, 57]. Many splice variants are tissue-specific [58], suggesting that it is entirely possible that tissue-specific DMMs in repeats affect alternative splicing and subsequent gene expression. Also, some DM genes might potentially be regulated simultaneously by their own methylation and that of repeats [53, 59], due to proximity with DM repeats.