Salt Stress Induces Non-CG Methylation in Coding Regions of Barley Seedlings (Hordeum vulgare)

Abstract

Salinity can negatively impact crop growth and yield. Changes in DNA methylation are known to occur when plants are challenged by stress and have been associated with the regulation of stress-response genes. However, the role of DNA-methylation in moderating gene expression in response to salt stress has been relatively poorly studied among crops such as barley. Here, we assessed the extent of salt-induced alterations of DNA methylation in barley and their putative role in perturbed gene expression. Using Next Generation Sequencing, we screened the leaf and root methylomes of five divergent barley varieties grown under control and three salt concentrations, to seek genotype independent salt-induced changes in DNA methylation. Salt stress caused increased methylation in leaves but diminished methylation in roots with a higher number of changes in leaves than in roots, indicating that salt induced changes to global methylation are organ specific. Differentially Methylated Markers (DMMs) were mostly located in close proximity to repeat elements, but also in 1094 genes, of which many possessed gene ontology (GO) terms associated with plant responses to stress. Identified markers have potential value as sentinels of salt stress and provide a starting point to allow understanding of the functional role of DNA methylation in facilitating barley’s response to this stressor.

1. Introduction

Barley is an important crop for food, feed and brewing [1,2] and is used as a research model for temperate cereals [3,4]. Although considered relatively tolerant to salinity [5], barley grown under saline conditions often suffers substantial yield losses [6]. In recognition of a global increase in saline soils worldwide [5], there are continuing efforts to improve the salt tolerance of barley varieties to maintain current levels of production. As with other plant species, barley responds to salt stress through the coordination of processes that alleviate both osmotic stress and ion toxicity [7]. Acclimation to saline conditions requires the stimulation of multiple molecular networks, including stress sensing, signal transduction, and the expression of stress-specific genes and metabolites [3,7,8,9]. Modern genetic improvement strategies aimed at improving salt tolerance require the characterisation of genes activated in response to saline stress [10] and ideally, better understanding of their interactions and of any plasticity in their expression afforded by epigenetic regulation [11].

Epigenetic mechanisms that control gene regulation act independently of any changes to DNA sequences [12,13,14], although one, DNA methylation, does alter its chemistry. The term DNA methylation describes the addition of a methyl group to a specific cytosine base in DNA, and this change often plays a critical role in moderating gene expression [13,15,16]. Indeed, DNA methylation has been implicated in several critical aspects of plant development and in regulating a plant’s adaptation to stress [13,17,18,19,20]. Change in DNA methylation status can occur via de novo DNA methylation, which is generally associated with gene repression, or by demethylation, which usually enhances gene expression [16], although numerous exceptions to this rule are known [16,21,22]. There are several reasons for characterising changes to the global methylation status of the genome that occur in response to a stress such as excessively saline soil. At the simplest level, identifying salt-induced methylation changes to specific sites has the potential to diagnose the presence and level of salt stress experienced by roots, based solely on the methylation status of key epimarkers. Salt concentration in saline soils varies on both spatial and temporal scales [23,24] and so measuring the timing and extent of exposure of an individual plant can be difficult in field conditions. The possibility of being able to characterise fine scale salt exposure of individual (sentinel) plants based on changes to the methylome is therefore an attractive prospect, and one that may also facilitate ready identification of genotypes exhibiting differential responses (e.g., tolerance to osmotic stress, Na⁺ exclusion from leaf blades, tissue tolerance) [5]. At the same time, better identification of those genic regions most likely to be methylation-regulated in response to salt stress provides a useful starting point from which to identify candidate genes that may be implicated in plastic salt stress responses and to build broader understanding of the molecular mechanisms in play that confer plant resilience to saline stress, something that has the potential to open up new avenues for crop breeding [11].

Several studies have demonstrated that exposure to salt stress can significantly perturb plant methylation profiles [25,26,27,28]. Others have correlated stress-induced modifications to DNA methylation to changes in gene regulation across a range of species [14,18,29,30], although some controversy remains over the consistency of the DNA methylation sites described [28,31]. In general, most salt-induced changes to DNA methylation seem to occur within, or in proximity to, known stress response genes [7,25,32,33]. In maize, salinity induced de novo methylation to zmPP2C in roots but demethylation to zmGST in leaves, with both changes seemingly linked to altered expression levels [34]. De novo methylation significantly repressed the expression of zmPP2C in roots, whereas demethylation of zmGST enhanced its expression in leaves, implying that DNA methylation changes in response to salt stress might contribute to stress acclimation [34]. In barley, acute salt stress has been similarly shown to evoke methylation-modulated changes to the expression of several genes involved in metabolic and physiological processes implicated in the ability of plants to cope with stress [3,9,35]. However, to date, there has been a marked lack of reports linking salt-induced gene expression to global changes in DNA methylation or methylation-associated changes that apply across a representative sample of any crop species [36].

For food crops with large genomes, the use of genome-wide bisulfite sequencing to characterise genome-wide flux in methylation from a representative range of genotypes is effectively precluded by cost and the complexity of bioinformatics [37]. For this reason, most works on stress-induced methylome change have elected to either target particular loci [29,38] or else to survey only a subset of the genome. Of the many methods available, Methylation Sensitive Amplification Polymorphism (MSAP) analysis has proved particularly popular to study stress-induced changes to genome-wide methylation patterns [25,27,32], in part because of the reproducible reputation of the technique [39,40,41]. However, the MSAP method only generates relatively small numbers of anonymous markers [42,43] and so has limited utility for studies aiming to establish links between changes in methylation and altered gene expression. While some works have sought to overcome this limitation by targeted sequencing of MSAP amplicons [7,25,32,33], others have argued that this amendment of the method is still cumbersome, costly and time-consuming [44]. The ability of Next Generation Sequencing to analyse large numbers of loci in multiple methylomes in parallel provides the opportunity to overcome these limitations. The use of methylation-sensitive Genotyping-By-Sequencing (ms-GBS) provides workers with the possibility of identifying differentially methylated markers (DMMs) with a better depth and coverage of the genome [44,45]. By using methylation-sensitive restriction enzymes to reduce genome complexity during library preparation, differentially methylated fragments are produced that are appropriate for high throughput sequencing [44,45]. This approach presents the advantage of detecting methylated sites that are dispersed across the genome and is particularly appealing for species with large genomes such as barley [44].

In this study, we used ms-GBS to assess the level of salt-induced changes to methylation site distribution patterns in the roots and leaves of five diverse barley genotypes (Barque 73, Flagship, Hindmarsh, Schooner and Yarra) with a range of salinity tolerance levels, and matched these against the reference barley genome to characterize the genomic locations of changed loci. We then combined these results with publicly available data about the gene expression of barley roots under salt to postulate the possible functional implications of DNA methylation flux on gene regulation in barley under salt stress.

2. Results

2.1. Methylation-Sensitive Genotyping-By-Sequencing (ms-GBS)

Overall, we generated in excess of 1 billion raw reads (1,015,703,602) from ms-GBS libraries, sequenced on a HiSeq 2500 (Illumina, San Diego, CA, USA). A high proportion of the raw reads passed the filter for the presence of the barcoded adapter, the MspI restriction product site and the EcoRI adapter (1,004,318,258; 98.87%). However, when these reads were filtered further to identify those that were uniquely mapped to the draft barley reference genome [4], the numbers fell substantially to 496,960,365 reads (i.e., 49.48% of raw reads). This yielded an average of 2,484,801 high quality reads per library and represented 892,859 unique sequence tags. Tags represented in this set amounted to 31.56% of the MspI recognition sites (5′-CCGG-3′) estimated for the barley reference genome (2,828,642;

Table 1. Data yields of the methylation-sensitive Genotyping-By-Sequencing (ms-GBS), generated using the Illumina HiSeq 2500 platform.

Sequencing Results	Yield
Raw reads	1,015,703,602
Reads that matched barcodes	1,004,318,258
Reads aligned to barley reference genome	496,960,365
Samples	200
Average reads per sample	2,484,801
Total unique tags	892,859
Polymorphic tags	645,297

).

2.2. Salt-Induced DNA Methylation Changes Are Organ and Concentration Specific

In total, 24,395 and 3777 unique sequence tags were deemed “significant differentially methylated markers” (false discovery rate (FDR) < 0.01) in leaf and root samples respectively across all five varieties and salt treatments (Figure 1 and Figure 2a). Curiously, the overall number of leaf DMMs increased progressively with salt concentration (75, 150 and 200 mM NaCl), whereas there was no such pattern seen in the roots (Figure 1). Soil salinity was found to induce more hypomethylated DMMs than hypermethylated DMMs in both leaves and roots, regardless of concentration (Figure 1). Although the number of salt-induced DMMs was higher in leaves (24,395 DMMs) than roots (3777), the scale of the change evoked by salt stress was far higher in roots when measured by p-values (Figure 2a) and the fold-changes in read counts (Figure 2b,c). A comparison of the median fold-change in methylation across all markers in the two organs revealed that salt induced net hypomethylation in roots and hypermethylation in leaves (Figure 2a–c), even though the number of salt induced hypomethylated sites exceeded those of hypermethylated sites in both organs (Figure 1).

2.3. Stability of Salt-Induced DMMs across Treatments

We next surveyed the appearance of DMMs across treatments and organs. Only a small proportion of DMMs appeared across all salt concentrations (Figure 3a,b). Moreover, of the 24,395 salt-induced DMMs detected in leaf samples, 52% were specific to 75, 150 or 200 mM NaCl (2390, 4070 and 6202, respectively) (Figure 3a), implying a positive association between the salt concentration and the number of loci affected by methylation changes. In roots, there was no obvious relationship with salt concentration, with 633, 1642 and 88 salt concentration-specific DMMs for 75, 150 and 200 mM NaCl, respectively (Figure 3b) (62% of the total).

There were, nevertheless, many stable DMMs that appeared in all salt concentrations. These dose-insensitive DMMs accounted for 22.9% (5593 of 24,395) of all salt-induced DMMs recovered from leaves and 14% (528 of 3777) of those recovered from roots (Figure 3a,b, Supplemental Data Set S1). These dose-insensitive DMMs invariably presented the same directionality of methylation change across all concentrations (i.e., always hyper- or hypomethylated) (Figure 4a,b). The dose-insensitive DDMs followed the global trend (see above) and so, mostly became hypomethylated following salt exposure in both leaves (4744, 84.82%) and roots (329, 62.31%). Of these, just 22 were shared between leaf and root samples, most of which, again, became hypomethylated following salt exposure (Figure 4c). Of these, 20 markers shared the same directionality of methylation change following salt exposure between organs, but two markers (“2:1:467135271” and “6:1:259709553”) became hypermethylated in leaves and hypomethylated in roots following exposure to salt (Figure 4c).

2.4. Distribution of Salt-Induced DMMs around Annotated Genomic Features

We assessed the distribution of concentration-independent, salt-induced DMMs relative to annotated features of the barley genome (e.g., protein coding genes, repeats, tRNAs, etc.). Proximity to a repeat sequence appeared to be a strong predictor of the location of DMMs induced by salt. Indeed, 96.5% of DMMs induced by salt in leaves and 99.8% in roots occurred either within the repeats themselves or within 1 Kb of them (Figure 5a,b).

We next sought to identify genes positioned within the proximity of the dose-insensitive, salt-induced DMMs. The expression of these genes was considered most likely to be consistently influenced by salt-induced methylation flux. In leaves, 19.1% (1070/5593) of dose-insensitive DMMs were located within 5 Kb of genes (Figure 5c; Supplemental Data Set S2), with the majority located within the gene body itself (56.4%, 603 DMMs; Figure 5c). In roots, just 24 (i.e., 4.5%) of the dose-insensitive DMMs lay within 5 Kb of a gene, five of which were located within the gene body, 14 were upstream and five were downstream (Figure 5d; Supplemental Data Set S2). Additionally, it is worth mentioning that of the 22 dose-insensitive DMMs shared in leaves and roots (Figure 4c), only one was positioned within 5 Kb of a gene (3994 bp upstream MLOC_63677 on chromosome 2H).

Given that the effect of DNA methylation on gene expression may depend on the position of the change relative to the transcribed sequences [16,46], we further investigated DMM distance to 5′UTRs, 3′UTRs, and exons of differentially methylated genes in leaves and roots. In leaves, it appeared that salt-induced DMMs near 5′UTRs were most abundant within 1 Kb (277 DMMs) of the 5′UTR in the downstream direction, with those falling between 1 and 2 Kb being the second most common (120 DMMs; Figure 6a). Outside these windows, DMMs occurred in the range of 40–65 DMMs per Kb (Figure 6a). DMMs were more common in the upstream direction of 3′UTRs, with the 1 Kb bin immediately upstream containing the highest number of DMMs (197 DMMs), decreasing gradually to reach background levels (50–70 DMMs per KB) after 4 Kb (Figure 6b). In comparison, there were insufficient gene-associated DMMs from root samples to provide strong evidence of clustering around either the 5′UTRs or 3′UTRs.

The majority of DMMs within gene bodies from leaf samples lay within exons (81.4%, 498 of 612; Figure 6e). The remaining DMMs were generally within 1 Kb of an exon (Figure 6e). Three out of the five gene body DMMs from roots were similarly exonic or within 1 Kb of the nearest exonic region (Figure 6f). Considered collectively, gene body DMMs were most commonly associated with the first exons (57.5%; 355/617), and included 296 overlaps, 45 downstream and 14 upstream (Figure 6e,f). Additionally, there were 41 DMMs from leaves and two DMMs from roots that clustered around tRNA genes (Figure 6g,h). While only one DMM overlapped with a tRNA in leaves, 14 out of the 41 DMMs were within 1 Kb upstream (nine DMMs) and downstream (five DMMs) of a tRNA gene (Figure 6g). The two DMMs nearest tRNA genes in roots were located within 1 and 4 Kb downstream of the gene respectively (Figure 6h).

2.5. Gene Ontology Analysis of Salt-Induced DMMs

A gene ontology (GO) analysis was performed for all differentially methylated genes (i.e., within 5 Kb of a salt-induced DMM) from both leaves and roots. The 1070 differentially methylated (DM) genes identified from leaves included 1017 that were hypomethylated and 53 that were hypermethylated following salt exposure. These genes yielded 433 and 99 high level GO terms, for the hypomethylated and hypermethylated groups, respectively (

Table 2. Number of differentially methylated genes (DMGs) and associated GO terms in barley leaves and roots. GO, gene ontology; hypo, hypomethylated genes; hyper, hypermethylated genes. GO groups were determined using REVIGO (Available online: http://revigo.irb.hr/).

Methylation Change	DMGs	GO Terms per GO Group			Total GO Terms
		Biological Process	Cellular Component	Molecular Function
Leaf hypo	1017	315	40	73	433
Leaf hyper	53	64	21	14	99
Root hypo	15	19	10	0	29
Root hyper	9	13	11	0	24

). The top five functional groups retrieved from the hypomethylated genes in leaves were the “protein modification process”, “cellular amide metabolism”, “cell cycle” and “negative regulation of signal transduction” (Figure 7a, Supplemental Data Set S3). Hypermethylated genes were enriched with GO terms that were associated with “organophosphate biosynthesis”, “peptide metabolism”, “peptide metabolism transport chain”, “generation of precursor metabolites and energy”, and “photosynthesis” (Figure 7b, Supplemental Data Set S3).

In roots, salt-induced hypomethylated markers were associated with 15 genes, whereas hypermethylated DMMs were in, or were proximal to, nine genes. These genes were significantly enriched for 29 (hypomethylated) and 24 (hypermethylated) GO terms (Table 2). The GO terms derived from hypomethylated genes in roots fell into three main function groups, in this order: “generation of precursor metabolites and energy”, “peptide metabolism” and “carbohydrate derivative metabolism” (Figure 8a, Supplemental Data Set S3). Hypermethylated genes enriched GO terms that were related to one main biological function: “peptide biosynthesis”. The details concerning all GO terms enriched by differentially methylated genes in roots are listed in Supplemental Data Set S3.

These GO terms, enriched from differentially methylated genes, give an indication of the biological pathways in which activity might be modified in response to salinity. Some GO terms, although not dominant, are related to functions essential for plant responses to salt stress, such as “ion transmembrane transport”, “potassium ion transport”, “cation transmembrane transporter activity”, “response to osmotic stress, “response to chemical stimulus”, “oxidation–reduction process”, “regulation of innate immune response”, “cellular response to stress”, and “defence response”, among others (Supplemental Data Set S3).

2.6. Differentially Expressed Genes in Barley Roots

To investigate whether observed changes in DNA methylation could be associated with changes in gene expression, salt-induced DMMs were compared to publicly available gene expression responses to salt exposure. These datasets were related to two genotypes (Sahara and Clipper) and included four biological replicates of each variety (see Section 4). Differential gene expression between salt treatments revealed 124 upregulated and 34 downregulated transcripts (

Table 3. Number of differentially expressed (DE) genes and associated GO terms in barley roots. GO, gene ontology; GO groups were determined using REVIGO (Available online: http://revigo.irb.hr/).

Expression Change	DE Genes		GO Terms per GO Group			Total GO Terms
	Total Transcripts	Annotated	Biological Process	Cellular Component	Molecular Function
Upregulated	124	76	94	22	29	145
Downregulated	34	18	23	12	0	53

, Supplemental Data Set S4) among which, 76 and 18 transcripts, respectively, matched barley reference genes in the public database “Ensembl” (Available online: http://plants.ensembl.org/biomart/martview). The ontology of these annotated genes revealed many pathways that were regulated by salinity in barley roots. The top five gene representatives of significantly enriched GO terms in upregulated genes were “organophosphate biosynthesis”, “peptide metabolism”, “protein modification process”, “electron transport chain”, “monovalent inorganic cation transport” and “photosynthesis” (Figure 9, Supplemental Data Set S4). Downregulated genes enriched GO terms which clustered around the functional pathway “peptide metabolism” and to a small extent, around “generation of precursor metabolites and energy”.

We then cross-referenced the differentially expressed (DE) genes against the DMMs identified in the current study. This was achieved by searching for DE genes within 5 Kb of DMMs. There were no differentially methylated genes amongst DE genes, with a false discovery rate (FDR) below 5%, and so, we extended the gene list by reducing the stringency of the FDR cut-off to 10%. With this setting, seven DE genes were found to be differentially methylated, one of which contained two DMMs (MSTRG.43260, one hypo- and one hypermethylated) (

Table 4. List of differentially methylated DE genes in barley roots. DE, differentially expressed gene; DMM, differentially methylated markers, Chrom, chromosome; FDR, false discovery rate; dist2Gene, DMM position relative to gene.

DE Genes		DMMs			Statistics				Annotation
Gene ID	Range	Chrom	Position	Methylation	logFC	p-Value	FDR	dist2Gene
MSTRG.4246	1:435681474–435731845	1H	435,689,351	hyper	−1.76	0.000	0.053	0	-
MSTRG.31525	5:507135444–507397451	5H	507,332,872	hypo	−1.47	0.002	0.083	0	MLOC_2917
MSTRG.43260	7:427906474–427974581	7H	427,925,930	hyper	−1.05	0.006	0.093	0	MLOC_73155
MSTRG.43261	7:427906474–427974581	7H	427,948,871	hypo	−1.05	0.006	0.093	0	MLOC_73155
MSTRG.10572	2:543673444–543674117	2H	543,678,039	hypo	−1.05	0.006	0.095	3922	-
MSTRG.6485	2:17425326–17624569	2H	17,517,122	hypo	1.39	0.007	0.095	0	-
MSTRG.6418	2:15418194–15419914	2H	15,414,469	hypo	1.53	0.004	0.089	3725	MLOC_24124
MSTRG.10644	2:545135370–545135958	2H	545,131,372	hyper	3.43	0.003	0.086	3998	MLOC_48766

). However, there was no correlation between gene methylation status and the direction of gene expression. Some hypomethylated genes were downregulated, whereas others were upregulated, and vice versa for hypermethylated genes (Table 4). Only four of these differentially methylated transcripts matched with annotated barley genes in public databases. The gene ontology analysis of these genes revealed that hypomethylated and hypermethylated genes enriched functionally close GO terms, which were all related to cellular components: plastid, cytoplasmic part and intracellular membrane-bounded (Supplemental Data Set S4).

3. Discussion

A growing number of studies are highlighting the role of DNA methylation in the coordination of the adaptive responses of plants to stress [14,30,31]. The primary challenge, particularly for crops with large genomes, resides in assembling a genome-wide picture of the role of methylation in orchestrating the molecular response to the stressor. As with the study of other stresses, works on methylation-based responses to salt stress have, therefore, largely relied on low throughput-targeted approaches, low genome coverage or anonymous markers applied to a low number of genotypes [7,25,26,27,33,47]. However, our use here of methylation sensitive Genotyping-By-Sequencing (ms-GBS) to study salt-induced changes in DNA methylation in ^mCCGG contexts has allowed us to survey methylome flux across a reasonably representative portion of the genome (Figure 2). The application of this approach allowed us to characterize genotype-independent, salinity-induced methylation flux in both leaf and root samples, and then to relate the pattern of differentially methylated markers to specific genomic features.

3.1. Consistency of Salt-Induced DMMs

The five barley varieties included in this study (Barque 73, Flagship, Hindmarsh, Schooner and Yarra) were selected based on the similarity in phenology. The lines varied however in their salinity tolerance, in terms of growth rates and sodium accumulation [48]. By comparing varieties that grow at similar pace we aimed to minimise epigenetic variability between samples associated with developmental differences. The high prevalence of concentration-specific DMMs in leaves (52%) and roots (62%) could be taken to imply that many or most are merely stochastic, statistical outliers, as suggested in previous studies [25,28]. However, it should be remembered that all DMMs described here were conserved across all five diverse barley genotypes and five biological replicates (per variety), with the direction of salt-induced methylation flux being conserved in all cases. This element of the experimental pipeline was introduced to minimize the effect of stochastic noise and was intended to strongly enrich conserved responses to salt stress. There are several aspects of the resulting data that suggest that this action did uncover at least some consistent and robust epimarks of salt exposure. At a simplistic level, the progressive increase in the number of leaf DMMs as the salt concentration increased could be taken as being suggestive of an incremental response as the salt concentration rose. This pattern would appear to be in accordance with the established theory that a large number of DMMs only become activated above a threshold concentration of salt, as hypothesized by Soen and co-workers [49]. Following this reasoning, as the salt concentration increases, more thresholds are exceeded and so, more DMMs become recruited into the global methylation flux. In this way, DMM abundance increases proportionally to the salt concentration. However, it is important to note that the use of only three concentration points in the titration series provides somewhat limited scope for confidence in this explanation. It is also pertinent to note that this pattern was not repeated among the root materials, and that DMMs were less frequent in roots, despite being exposed to higher levels of salt than leaf tissues [50]. However, the 22% (6121/28,172) of DMMs conserved across all salt concentrations are likely to be far more robust. The fact that these dose-insensitive DMMs invariably exhibited the same directionality of methylation change across all concentrations (i.e., always hyper- or hypomethylated) despite differing between DMMs is compelling evidence that most are truly salt-induced DMMs that appear across diverse genotypes and respond to salt exposure in a consistent manner.

3.2. Salt Induces Different Changes to the DNA Methylation of Leaves and Roots

The observation that the vast majority of salt-induced, dose-insensitive DMMs were also organ-specific warrants consideration. Roots and shoots respond differently to salinity, with roots playing a role in Na⁺ sequestration to reduce the amount of Na⁺ reaching the shoot where it can inhibit photosynthesis. As such, it is likely that different cellular processes are activated in the different organs. Another plausible explanation for the relative paucity of DMMs shared between roots and leaves can be taken from the fact that, although most research on crop tolerance to salt stress has focused on the effect of Na⁺ [51,52,53], it is now well stablished that the toxic effects and exclusion mechanisms of Na⁺ and Cl⁻ in barley are different and independent [54,55]. A parallel assessment of the barley epigenetic response to each independent ion (i.e., Na⁺ and Cl⁻) is crucial in order to fully understand the contribution of DNA methylation to salt stress. Barley can effectively compartmentalize excess Na⁺ in the roots and inhibit its transport to the leaves [56]. This response gives rise to very different Na⁺ environments in the two organs, with Na⁺ concentrations being generally higher in the roots than the shoots [50,56]. We reason that the dose-insensitive DMMs identified in the present study are all likely to respond to low salt thresholds, and so, appear in all salt treatments across the titration. Even though Na⁺ levels differ between the roots and shoots [50], the levels in both organs are generally within the same order of magnitude, and a low induction threshold would be reached in both. However, the divergence between leaf and root DMMs is in accordance with the known physiological response of the species to saline exposure of the roots and implies that both sets of DMMs would, therefore, provide a robust indication of exposure to salt.

Certainly, it has been widely reported that salinity imposes extensive, genome-wide modification of DNA methylation patterns, with more methylation changes being reported in leaves when compared with roots [25,26,27,33,57,58,59,60], a trend that is clearly in accordance with our findings. Taken at face value, the greater abundance of salt-induced methylation changes in leaves than in roots appears counterintuitive, since roots are in direct contact with salt stress. That said, we found the scale of change in methylation was greater in roots than in leaves, suggesting that although salt evokes changes in fewer loci, the effects on these sites are greater. Provided these changes are associated with concurrent changes to expressions of key genes involved with responses to salt stress, these observations can accommodate root-specific epigenetic responses to saline environments while plants are undergoing osmotic stress and salt toxicity [8,61]. Should at least one of these processes be focused on repressed transport of Na⁺, then milder ion accumulation in the leaves will slowly increase stress [5,62] at a lower level because of ‘leakiness’ of the system, evoking a widespread, but more measured, response in the leaves.

Our results agree with many previous studies that have reported that the overall level and direction of methylation flux in response to salinity varies according to organ type, with a tendency towards hypomethylation in roots and hypermethylation in leaves [7,25,26,27,33,47]. However, we also noted that the proportion of de novo methylation and demethylation events varied in the same manner in both roots and leaves, with a prevalence of hypomethylated events in both organs, albeit at different frequencies. It is possible that divergence between our findings and those of previous studies [7,25,26,27,33,47] may simply be a feature of barley. However, it is also possible that the trend towards hypomethylation is a more general one and that our findings diverge because of methodological differences in the present work such as (1) the high-throughput sequencing used to generate methylation profiles; (2) the level of stringency in selecting DMMs (FDR < 0.01); and 3) the diversity of barley varieties used in this study to account for genotype-dependent DNA methylation [25,26,27]. Most studies of salt-induced DNA methylation have relied on MSAP analysis of a single variety to assess flux in DNA methylation [7,25,26,27,33,47]. However, MSAP generates anonymous markers, and it is difficult to interpret the gain or loss of markers as providing clear evidence of hypo- or hypermethylation [42].

3.3. Salt-Induced Changes in DNA Methylation May Influence Gene Regulation

DNA methylation is modulated in the genome in three ways: de novo methylation (hypermethylation), methylation maintenance, and methylation removal (hypomethylation) [63]. Modification of DNA methylation in response to stress is hypothesised to be at least partially directed to specific genomic regions where the DNA methylation status influences the expression of stress-response genes [18,30,64,65]. This is in accordance with our finding that dose-insensitive salt-induced DMMs appeared more commonly in sites that could facilitate perturbation in the expression of genes that hypothetically could be part of a molecular response to salinity, as reported elsewhere [65]. There is evidence from previous studies to suggest that salt-induced DMMs can play important roles in evoking metabolic differences between seedlings growing under control and saline conditions [25,26,29,34,64,66], although the markers found in these works were either few in number or else identified from a single genotype. The provision here of a robust list of consistent, salt-responsive DMMs therefore provides a useful starting point from which to gather a more holistic picture of DNA methylation-mediated regulation of molecular responses to salt exposure.

Proof of a functional link between the change in methylation status in these DMMs and associated alterations in the expression of proximal stress response genes is beyond the scope of the current study. Nevertheless, there are several grounds for reasoning that at least some of the markers identified here are indeed functionally important. Certainly, others have argued that the close proximity of DMMs relative to the target genes is at least one requirement for such a relationship [19,67,68,69]. Viewed in this context, our observed clustering of DMMs around untranslated regions (UTRs) and exons is at least consistent with salt-induced DMMs mediating a functional response to the stress. Others have shown that a high frequency of salt-induced DMMs in gene extremities (towards 5′UTR and 3′UTR) can influence gene regulation through 5′UTR’ and 3′UTR’ closed-loop regulation systems which generate inactive transcripts [70,71], or through independent gene regulation by each UTR type [72]. Karan et al. [25] similarly observed that salt-induced DNA methylation changes generally occur in exon and UTR regions and could affect diverse biological functions in plants. There is also a strong body of evidence to suggest that gene body methylation, in general, can affect gene expression [19,55,68] by enhancing or inhibiting transcription and translation processes [71,72].

It has been claimed that, of all cytosine contexts, only ^mCG methylation occurs within gene bodies [68,73,74,75]. Our findings and those of others [76] do not support this stance, with non-CG types of methylation, such as ^mCCGG, being found frequently in transcribed regions from DNA isolated from both leaves and roots of barley. It is, however, still open to question whether these markers, like ^mCG, play roles in regulating gene expression [77]. Our observation of salt-induced DMMs associated with tRNA genes is more surprising and perhaps in accordance with the suggestion of a role for methylation-dependent regulation to support the RNA quality control system and protein synthesis [78,79,80]. More work is clearly required to investigate this possibility.

3.4. Salt-Induced DMMs Correlate with Stress Related Genes

There is circumstantial support to argue that at least some of the DMMs identified here may play functional roles in the expression of salt-response genes. Salt stress in barley alters the expression patterns of genes involved in a diverse range of physiological and regulatory pathways [3,9]. Given that salt-induced DMMs have the potential to regulate gene expression, the functions of differentially methylated genes were explored for possible correlations with stress responsive genes. The correlation of DM genes with GO terms that are related to plant responses to stress, such as “negative regulation of signal transduction”, “photosynthesis”, “response to osmotic stress” and “ion transmembrane transport”, is at least consistent with the possibility that salt-induced DMMs targeted genes which could play active functions in a plant’s response to salt, in broad accordance with previously expression studies of salt response [29,81,82,83]. In addition, some of the DM genes identified here were enriched for GO terms, such as “hydrolase activity”, “oxidoreductase activity”, “nucleic acid binding”, and “translation factor activity”. This agrees with similar reports of differentially methylated genes associated with salt stress in rice [25,26].

The present study also revealed differential DNA methylation of genes implicated in “organophosphate biosynthesis”. Should the change in methylation state alter the expression of these genes, it would be in accordance with previous studies which have shown that salt stress induces an increase in the amount of intra-cellular organophosphate solutes, such as di-myo-inositol-phosphate, inositol (1,4,5) trisphosphate, b-mannosylglycerate, b-mannosylglycerate and glutamate [84,85]. Furthermore, it was reported that salinity induces inorganic phosphate toxicity when Pi exceeds 0.10 mM in the substrate [86,87]. This salt-induced phosphate toxicity may arise from an excess of phosphate, not only due to P uptake, but also due to salt-induced increases of intracellular organophosphate solutes [85,86]. However, it is important to recognise that although the presence of DMMs near a gene may be an indication of responsiveness to salt stress, it does not provide sufficient evidence of a functional role of DNA methylation in the regulation of that gene [21,68]. Gene expression analysis is required to assess the link between DNA methylation and gene activity under salt stress. Finally, our current understanding of plant epigenetic responses to the environment has predominantly arisen from the analysis of mixed populations of cells contained within plant organs. Recent works have identified tissue [77] and cell type-specific [88] patterns of DNA methylation associated with tissue specific gene expression [77]. In light of these results, it is clear that in order to achieve a better understanding of the role that DNA methylation plays in plant responses to environmental challenges, future studies should include approaches that independently interrogate the methylome and transcriptomes of different tissues and cell types.

The difficulty in extrapolating a functional link is highlighted by the fact that only seven of the genes previously reported to be differentially expressed in roots [3] were also differentially methylated under salt stress. This result may imply that few of the markers found here are functionally important or, alternatively, this may be attributable to the use of different biological samples for methylation profiling and gene expression analyses, different growing conditions [3], partial coverage of the barley reference genome used here [4] or to the possibility of bias due to salt-induced DNA degradation [33,54,89,90].

4. Materials and Methods

4.1. Plant Material and Stress Treatment

Five diverse spring barley varieties were used in this investigation: Barque 73, Flagship, Hindmarsh, Schooner and Yarra. Seeds were kindly provided by the Salt Focus Group at the Australian Centre for Plant Functional Genomics (ACPFG, Adelaide, Australia). The lines were chosen due to their similar phenology, but they had a range of salinity tolerance levels [48]. The experiment was designed in randomized blocks of five replicates and four salinity treatments: control (0), 75, 150 and 200 mM NaCl.

Seeds were germinated, and seedlings were grown in 3.3 L free-draining pots, placed on saucers containing 2915 g of growth substrate (50% UC (University of California at Davis, Davis, CA, USA) potting mix, 35% coco-peat, and 15% clay/loam (v/v)). The five barley varieties were sown in each pot and positions were randomized in each pot to minimize a block effect. Two seeds were sown per variety and thinned to one seedling at 8 days after sowing. Salinity treatments were applied 10 days after sowing in four increments over 4 consecutive days, to minimise osmotic shock [91]. The required amount of NaCl for each salt concentration was calculated based on the substrate soil’s dry weight and the target gravimetric water content of 16.8% (g/g) [91]. At the time of salt application, the water content reached 26.4% and dropped down to the final concentration through evapotranspiration. Pots were watered to weight every 2 days to maintain the target gravimetric water content (16.8% (g/g)) [91] until sampling.

This experiment was conducted from 30 January to 20 February 2015, in a greenhouse at the Waite Campus, University of Adelaide (Adelaide, Australia), (34°58′11″ S, 138°38′19″ E). The seedlings were grown under a natural photoperiod and the temperature was set at 22 °C/15 °C (day/night).

4.2. DNA Extraction

At day 11 after the first salt stress imposition to barley seedlings (21 days after sowing, three leaves stage), 50 mg samples were collected from the middle sections of the 3rd leaf blades and roots. In total, 200 samples were collected (five varieties, four treatments and two organs), and were snap-frozen in liquid nitrogen, and then stored in a −80 °C freezer until needed for DNA extraction. Prior to DNA extraction, frozen plant material was disrupted in a bead beater (2010-Geno/Grinder, SPEX SamplePrep^®, Metuchen, NJ, USA). Genomic DNA was isolated using a Qiagen DNeasy kit (Qiagen, Dusseldorf, Germany), following the manufacturer’s instructions. DNA samples were then quantified in a NanoDrop^® 1000 Spectrophotometer (V 3.8.1, ThermoFisher Scientific Inc., Waltham, MA, USA) and concentrations were standardised to 10 ng/µL for subsequent ms-GBS library preparation.

4.3. Methylation Sensitive Genotyping by Sequencing (ms-GBS)

ms-GBS was performed using a modified version [44,45] of the original GBS technique [92,93]. Genomic DNA was digested using the combination of a rare cutter, EcoRI (GAATTC), and a frequent, methylation sensitive cutter MspI (CCGG). Each sample of DNA was digested in a reaction volume of 20 µL, containing 2 µL of New England BioLabs Smartcut buffer, 8 U of HF-EcoRI (High-Fidelity) and 8 U of MspI (New England BioLabs Inc., Ipswich, MA, USA). 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 (NEB T4 DNA Ligase #M0202) 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 (#A63880, Beckman Coulter, Brea, CA, USA) at a ratio of 0.85, 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 reactions (25 μL each), containing 10 μL of library DNA, 5 μL of 5× Q5 high fidelity buffer, 0.25 µL of polymerase Q5 high fidelity, and 1 μL of each Forward and Reverse common primer at 10 μM, 0.5 µL of 10 μM dNTP and 7.25 µL of pure sterile water. PCR amplification was performed in a BioRad T100 thermocycler (BioRad, Hercules, CA, USA), consisting of DNA denaturation at 98 °C (30 s) and ten 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 Agencourt AMPure XP magnetic beads following the manufacturer’s instructions. Bead-captured fragments were eluted in 35 μL of water, and 30 μL of the elution was collected in a new labelled microtube. Next, libraries were 125 bp paired-end sequenced in an Illumina HiSeq 2500 platform (Illumina Inc., San Diego, CA, 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) (See Supplemental Data Set S5 for sample sequencing information).

4.4. Data Analysis

The ms-GBS data was analysed following a workflow requiring bioinformatics tools in both Linux bash shell and R environments. Fastq files from the Illumina sequencing platform were first de-multiplexed and checked for read quality by the sequencing service provider, reporting read quality encoded in symbolic ASCII format as Phred-like quality score +33. Only fragments with at least 95% of the reads having Phred > 25 were retained. Reads that did not have a barcode were put into undetermined files and removed from any downstream analyses. Prior to demultiplexing, Illumina adaptor sequences, used for library construction, were also removed. The second step consisted of preparing the reads for alignment with the barley reference genome. As this was paired-end read sequencing data, both strands were merged together in a single read, using the module bbmap in bash. Merged reads were next aligned to the barley reference genome downloaded from the Ensembl database (Available online: http://plants.ensembl.org/Hordeum_vulgare/). This required the module bowtie/2-2.2.3 to build a bowtie2 index for the barley genome, and the module samtools/1.2 to perform alignments. As paired reads were merged into single reads, only those that overlapped were retained, to allow proper mapping. This alignment step yielded bam files containing only reads that matched with the reference genome. Next, a read count matrix was generated using only marker sequence tags that matched with MspI cut sites on known chromosomes (1H to 7H), and those on contigs were discarded. This count matrix was then used as data source for the performance of subsequent analyses using R packages.

4.5. Salinity Induced Differentially Methylated Markers in Barley

The alteration of DNA methylation in barley seedlings exposed to salinity was assessed in ^mCCGG contexts by the use of MspI during sample preparation. Differentially methylated markers (DMMs) were identified using the package msgbsR developed by Mayne et al. [94] (Available online: https://github.com/BenjaminAdelaide/msgbsR, accessed on 26 August 2016), in which a generalised linear model was fitted to the design with the trimmed mean of M-values normalisation option (TMM). Then, the Benjamini–Hochberg method was used to determine p-values. DMMs were selected based on FDR < 0.01 for differences in read counts per million between the salt-free control and salt treatments (75, 150 or 200 mM NaCl), with at least 1 count per million (CPM) reads. To obtain robust salt-induced markers, we selected DMMs that were conserved in all barley genotypes and present in at least 20 samples per treatment. The logFC (logarithm 2 of fold-change in CPM reads) was computed to evaluate the intensity of salt treatment-induced alteration of DNA methylation and to infer whether the change was a de novo methylation or demethylation event. This approach for determining the directionality of DNA methylation uses the fold change as an inverse proxy for changes in the methylation level. That is, higher methylation levels at a specific locus will reduce the number of restriction products for that locus [39] and therefore, reduce its number of CPM reads.

4.6. Distribution of Salt-Induced DMMs around Genomic Features

To determine whether there was a correlation between salt-induced DNA methylation and genomic features in barley, the distribution of DMMs was assessed around genes and repeat regions, as defined in the Ensembl database (Available online: http://plants.ensembl.org/biomart/martview/). This was done by mapping stable, salt-induced DMMs with repeats and genes in the barley reference genome. Then, we tallied the number of DMMs within genomic features (repeats, genes, exons) and per 1 Kb bins within 5 Kb flanking regions, both up- and downstream [47,95], using the shell module bedtools/2.22.0 [96]. The same procedure was repeated to estimate the number of DMMs around exons and UTRs of differentially methylated genes, and tRNA genes.

4.7. Gene Ontology of Differentially Methylated Genes

Genes within 5 Kb of a DMM were referred to as differentially methylated genes (DMGs). These genes were used for the gene ontology analysis, to investigate whether salt-induced changes in DNA methylation correlated with salt responsive genes. DM genes were grouped in hypermethylated and hypomethylated genes per organ (leaf or root), which were next used separately for GO term enrichment, using two R packages: GO.db and annotate [97,98]. Significant GO terms were selected based on Bonferroni adjusted p-values [99] at a significance threshold of 0.01 and a total GO enrichment of DM and non-DM genes equal to at least 10. The results of the GO analysis were visualized in treemaps generated in REVIGO [100].

4.8. Gene Expression and Ontology Analysis of Root Transcriptome

We further investigated whether differentially methylated genes were known to be differentially expressed in the plant. To do so, we used, as an example, a dataset of root transcriptomes of two barley varieties (Clipper and Sahara-3771) grown under salt stress (100 mM NaCl) and control conditions [3]. The raw data was downloaded from Available online: https://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-4634/, and samples from the root maturation zone, as defined by the authors [3], were used. The data contained four biological replicates of two varieties and two salt treatments (control and 100 mM NaCl), to give a total library size of over 390 million reads. Quality control was performed on these reads, which were then merged to form a single, large fastq file for each sample. Merged read pairs were trimmed using AdapterRemoval [101], followed by a second round of quality control.

After alignment using hisat2-2.0.4 in bash [102], a salt-induced differential gene expression analysis was performed, using a custom GTF file from Ensembl, created by the tool StringTie 1.3.1c [103]. This GFF file was restricted to transcripts on known chromosomes (1H to 7H). Read counts were assigned to genes in the GTF file using featureCounts v1.5.1 [104] and loaded as DGEList objects in R. As the data contained paired-end reads, the parameters were set to only count fragments (i.e., template molecules), instead of individual reads. This dataset was next filtered to keep only genes with CPM > 0.5 in at least four samples. Gene transcripts passing these conditions and present on chromosomes 1H to 7H were retained for differential expression analysis.

Before comparing treatments, the dataset was explored for sample variability using the MDS plot. Differential gene expression was then estimated using the lmFit function in limma::voom, a gene-wise linear model [105], and differentially expressed genes were defined as having an absolute fold-change > 2, with an FDR adjusted p-value < 0.05. Differentially expressed genes were first used “as is” for the gene ontology analysis as described above (previous section). Differentially expressed genes were then assessed for their proximity to salt-induced DMMs within 5 Kb in both directions. Genes found in this proximity with DMMs and referred to as differentially methylated DE genes, were used for another GO analysis. The results of these GO enrichments were visualized in treemaps produced in REVIGO [100] to show the main GO representatives.

5. Conclusions

To our knowledge, this study has provided the most comprehensive set of robust leaf and root epimarkers to indicate the exposure of barley to salt stress. These markers were conserved in both identity and direction across five diverse genotypes, biological replicates and all salt concentrations used. The leaf markers have potential value as epigenetic sentinels of the exposure of individual plants to soil salt stress. Viewed collectively, the root and leaf markers provide a useful starting point from which to assemble a more comprehensive picture of the functional role of DNA methylation in facilitating the plastic molecular responses of barley to this important stressor.