Data yields from ms-GBS, generated using the Illumina HiSeq 2500 platform.
The barley (Hordeum vulgare) genome comprises over 32,000 genes, with differentiated cells expressing only a subset of genes; the remainder being silent. Mechanisms by which tissue-specific genes are regulated are not entirely understood, although DNA methylation is likely to be involved. To shed light on the dynamic of DNA methylation during development and its variation between organs, methylation-sensitive genotyping by sequencing (ms-GBS) was used to generate methylation profiles for roots, leaf-blades and leaf-sheaths from five barley varieties, using seedlings at the three-leaf stage. Robust differentially methylated markers (DMMs) were characterised by pairwise comparisons of roots, leaf-blades and leaf-sheaths of three different ages. While very many DMMs were found between roots and leaf parts, only a few existed between leaf-blades and leaf-sheaths, with differences decreasing with leaf rank. Organ-specific DMMs appeared to target mainly repeat regions, implying that organ differentiation partially relies on the spreading of DNA methylation from repeats to promoters of adjacent genes. Identified DMMs indicate that different organs do possess diagnostic methylation profiles and suggest that DNA methylation is important for both tissue differentiation and organ function and will provide the basis to the understanding of the role of DNA methylation in plant organ differentiation and development.
- Hordeum vulgare
- tissue-specific methylation
- developmental epigenomics
DNA methylation is an important characteristic of plant genomes [1, 2], and can occur in all cytosine contexts (CG, CHG and CHH, where H = A, C or T) . The effect of DNA methylation variants on plant development has been demonstrated through methylation alteration tests, which led to plant abnormalities [4, 5]. Furthermore, DNA methylation has been reported to vary from tissue to tissue in many species [6, 7, 8, 9, 10], and these methylation changes seemed to be essential for normal plant development [11, 12].
Additionally, tissue-specific methylation was proposed to have a strong correlation with the differential expression of some tissue-specific genes. Examples include tissue-specific pigmentation in maize, which is reported to be under epigenetic control , and differential gene expression between organs attributed to differentially methylated regions in soybean  and sorghum . These studies extended our understanding of the functional importance of tissue-specific DNA methylation, including its role in setting developmental trajectories [9, 13, 15].
A substantial proportion of developmentally expressed genes have alternative promoters (multiple promoters that regulate the same gene) which are under different regulatory programmes . Maunakea et al.  proposed that alternative promoters are, at least sometimes, controlled by intragenic DNA methylation. This form of developmental gene regulation is reasoned to be dependent on transposon activity  and by implication would mean that silencing of transposons due to DNA methylation may be central to tissue-specific gene expression. Also, tissue-specific gene expression has been associated with methylation changes in promoter regions [2, 18, 19], especially CG islands within promoters . These studies indicate that tissue-specific gene expression does not rely on a single methylation pattern in the genome but, probably, on a combination of variable DNA methylation features.
The magnitude of differential methylation between tissues has been the subject of controversy. It was believed that significant distinctive DNA methylation existed only between specialised tissues such as endosperm, pollen, leaves and roots [9, 10, 21, 22]. Nevertheless, many of these studies also showed that differential DNA methylation between organs, such as roots and leaves, was minor in rice , maize , sorghum  and Arabidopsis . DNA methylation differences between roots and leaves were small in both mCG and mCHG contexts [9, 10], with about 1% and 5% divergence, respectively, reported in Arabidopsis . While these studies of differential DNA methylation between tissues generally compared the overall methylation levels [9, 10, 24], these results differ from comparisons made with differentially methylated markers (DMMs) between the same tissues , probably due to differences in methylation profiling methods, making it difficult to compare results from different studies. Therefore, it is difficult to know whether differences in the results concerning tissue-specific DNA methylation are due to the plant species or to the approach taken. The study of DNA methylation patterns in plant tissues is important for a better understanding of how these epigenetic markers determine tissue differentiation. Thus, further investigation is warranted to clarify organ specificity of cytosine methylation and the distribution patterns of tissue-specific DNA methylation markers in the plant genome.
To undertake such an investigation, we used barley, a globally important cereal crop, the genome of which has been sequenced recently . The availability of a reference genome made barley a model for the study of cereal crops such as wheat, oats or rye. In this study, we assessed differential DNA methylation between two barley (
2. Materials and methods
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−1) 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−1) (0.8 × field capacity)  until sampling 21 days after sowing, when seedlings were at three-leaf stage (Zadok stage 13 ). 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  of the original GBS technique . Genomic DNA was digested using the combination of a methylation-insensitive rare cutter,
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
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  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  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
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,
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 ) across all samples (Table 1). Of these reads, 99.27% contained the barcode and
|Reads that matched barcodes||895,013,295|
|Reads aligned to barley reference genome||448,445,748|
|Average reads per sample||2,637,916|
|Total unique tags||913,697|
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).
|Blade 1||Blade 2||Blade 3||Sheath 1||Sheath 2||Sheath 3|
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  and sorghum . Further, no significant difference was detected at all for mCG and mCHG methylation levels between tissues in cotton . 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) , and (2) the data analysis approach used, which can compare either global methylation levels (e.g. percent methylation)  or methylated loci (e.g. DMMs) . 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
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 . 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 , 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 . 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 . However, it is already well-established that tissue-specific DMMs can influence gene expression by enhancing gene transcription  and alternative splicing  or through repression due to immediate proximity to transcription start site .
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
This study provides a comprehensive set of robust tissue specific epimarkers which were conserved in all barley genotypes tested and can therefore be considered genotype independent. Such markers have potential to be converted into locus-specific methylation sensitive cleaved amplified polymorphic sequence markers (ms-CAPS) to be used as diagnostic of sample origin. Moreover, these markers provide a basis for the understanding of the role of DNA methylation in plant organ differentiation and development. Our data illustrates that during tissue development, DNA methylation evolves to reach a default profile once the tissue is completely differentiated at maturity. It is possible that the plant organ formation and maturation is under at least partial control of DNA methylation changes. In addition, repeats could play an important role in tissue definition. The existence of tissue-specific mCCGG sites suggests that this context carries important factors of tissue differentiation. Expression analysis of tissue samples would conclusively demonstrate the role of tissue-specific DMMs in gene regulation. These markers will provide a basis for future studies of the role of DNA methylation in plant organ differentiation and development.
The authors are grateful to the Australian Agency for International Development (AusAID) for supporting MK with an Australian Awards Scholarship. MJW was paid by the BBSRC (BBS/E/W/0012843C). This work is supported by the National Institute of Food and Agriculture, U.S. Department of Agriculture, Hatch Program number 2352987000, and award number 2019-67013-29168.
Conflict of interest
“The authors declare no conflict of interest.”
|Start||End||ID||Rank||Start||End||bp to exon|
MK conceived and performed the experiments, analysed the data and wrote the manuscript; BJM performed ms-GBS data alignments; MJW, ESS, BB, CMRL conceived the experiments and supervised the work. All authors read and commented on the manuscript.