Greenhouse Spatial Effects Detected in the Barley (Hordeum vulgare L.) Epigenome Underlie Stochasticity of DNA Methylation

Environmental cues are known to alter the methylation profile of genomic DNA, and thereby change the expression of some genes. A proportion of such modifications may become adaptive by adjusting expression of stress response genes but others have been shown to be highly stochastic, even under controlled conditions. The influence of environmental flux on plants adds an additional layer of complexity that has potential to confound attempts to interpret interactions between environment, methylome, and plant form. We therefore adopt a positional and longitudinal approach to study progressive changes to barley DNA methylation patterns in response to salt exposure during development under greenhouse conditions. Methylation-sensitive amplified polymorphism (MSAP) and phenotypic analyses of nine diverse barley varieties were grown in a randomized plot design, under two salt treatments (0 and 75 mM NaCl). Combining environmental, phenotypic and epigenetic data analyses, we show that at least part of the epigenetic variability, previously described as stochastic, is linked to environmental micro-variations during plant growth. Additionally, we show that differences in methylation increase with time of exposure to micro-variations in environment. We propose that subsequent epigenetic studies take into account microclimate-induced epigenetic variability.

DNA methylation has emerged as the prominent epigenetic signature for past or contemporary exposure of a plant to environmental insults (e.g. Xie et al. (2017) and has been implicated in the moderation of stress response (Bird and Jaenisch, 2003;Zilberman and Henikoff, 2007;Boyko and Kovalchuk, 2008). For instance, Tricker et al. (2012) reported that Arabidopsis thaliana responded to high relative humidity stress by suppressing the expression of two genes that control stomatal development through DNA methylation. DNA methylation has been similarly implicated in the response of various plant species to a range of stresses, including excess salt (Karan et al., 2012;Konate et al., 2018), temperature extremes (Steward et al., 2002;Bastow et al., 2004;Hashida et al., 2006;Pecinka et al., 2010;Song et al., 2012), herbivory (Herrera and Bazaga, 2011;Herrera and Bazaga, 2013), and heterogeneous environmental pressure (Wang et al., 2016). However, the relationship between DNA methylation and the stress effect is imprecise. Many of the methylation changes observed under stress fail to occur consistently across all genotypes or populations studied, and many others are not obviously associated with exonic regions. Fewer still can be directly tied to a particular stress response gene. Such observations have been described as stochastic (Karan et al., 2012;Tricker et al., 2012), spontaneous (Raj and Van Oudenaarden, 2008;Becker et al., 2011;Van Der Graaf et al., 2015), and without clear triggering factors (i.e. occurring randomly in the genome independently of stress). Many have considered the random and spontaneous alteration of DNA methylation is an adaptive biological process in its own right; one that drives diversity and evolution in a Lamarckian-like fashion (Feinberg and Irizarry, 2010;Meyer and Roeder, 2014;Soen et al., 2015;Van Der Graaf et al., 2015;Vogt, 2015) and with the clear potential to alter fitness (Consuegra and Rodrıǵuez Loṕez, 2016). Additionally, Soen et al. (2015) proposed a conceptual framework of random variations in the genome, instigated in response to environmental cues. They hypothesized that imposition of diverse types of stress upon individual organisms during development gives rise to an adaptive improvisation which deploys random phenotypic variations that allows some individuals to cope with unstable ambient conditions. However, the authors did not suggest an epigenetic mechanism that might be involved in the regulation of such adaptive phenotypic variation.
In a pivotal piece, Vogt (2015) provided insight into the concept of random variability. The author linked "stochastic developmental variation" to stochastic occurrence of DNA methylation (Bird and Jaenisch, 2003;Field and Blackman, 2003). However, Vogt did not consider in depth the possible role that microclimatic variation may play in this apparent stochasticity. Herrera and Bazaga (2010) suspected a role for mesoclimate in driving the epigenetic variability of natural populations but did not anticipate marked environmental differences to occur under controlled experimental conditions (greenhouse, growth room).
Moreover, since genome-by-environment interactions have been shown to be at least partially regulated by DNA methylation (Verhoeven et al., 2010), even minor perturbations of growing conditions attributable to positional effects within a controlled growing environment has the potential to introduce confounding variation in methylation patterning. One way of dealing with spatial variation, if it cannot be prevented, is to deploy an appropriate experimental design in order to distinguish treatment from positional effects (Brien et al., 2013;Cabrera-Bosquet et al., 2016). Experimental design normally accounts for such variability by combining blocking and randomization, along with appropriate statistical analyses (Addelman, 1970;Ruxton and Colegrave, 2011). Despite the usefulness of this approach, experimental design cannot entirely remove environmental variability (microclimate). This presents a potential challenge when attempting to link changes in DNA methylation to environmental stimuli. It is difficult to discriminate between the so-called stochastic methylation and position-dependent methylation due to the capacity of plants to promptly sense and epigenetically respond to subtle variation in ambient conditions (Gutzat and Mittelsten Scheid, 2012;Meyer, 2015).
In the present study, we combine methylation-sensitive amplified polymorphism (MSAP) and phenotypic analyses to assess the effect of microclimate on DNA methylation of barley plants growing under greenhouse conditions. To provide an indication of scale, we sought to compare the glasshouse positional effects on MSAP profiles and those generated after imposing mild salt stress to a replicate sample of plants grown in the same conditions. For this, nine spring barley varieties were grown in a randomized plot design under mild soil salt stress or control conditions. Environmental, phenotypic, and DNA methylation data collected at two time points are used to explore whether stochastic epigenetic may be linked to trivial environmental fluctuations. We also explore how phenotypic variability observed in these experiments correlates with differences in DNA methylation patterns.

Experimental Strategy
The central aim of this study was to assess the impact of microclimate (caused by differing plant positions within a glasshouse) on above ground biomass of barley plants, and of any associated change in global leaf-blade epigenome (as detected by MSAP). In this context, the MSAP profiles are being considered as a component of the DNA methylome (phenotype of the genome) and biomass is considered a component of their physical phenotype. We also sought to assess the scale of any changes to the epigenome through comparison with responses seen to mild salt stress among biological replicates in the same experiment. We sought to control possible sources of confounding variation (independent variables) by standardizing the source of material used for DNA extraction (tissue-to-tissue variation), use of a panel of varieties with similar growth rates (genetic variation), and collecting samples at two life stages (developmental variation).
Nine varieties of spring barley ( Table 1) were grown in a controlled temperature greenhouse at the Plant Accelerator ® (Australian Plant Phenomics Facility (APPF), Waite Campus, University of Adelaide, Australia) from June 26 to October 12, 2013. Varieties with similar flowering times (Menz, 2010) were selected to minimize discrepancies in sampling times between varieties. The experiment comprised eight randomized blocks with two plants of the same variety per plot ( Figure 1). Three seeds were sown in white pots (20 cm height × 15 cm diameter, Berry Plastics Corporation, Evansville, USA) containing 2400 g potting mixture (composed of 50% UC (University of California, Davis) potting mix, 35% coco-peat and 15% clay/loam (v/v)). Seedlings were thinned to one seedling per pot 2 weeks after sowing. Two soil salt treatments (0 and 75 mM NaCl ("control" and "salt stress," respectively, hereafter) were applied to three-leaf stage seedlings (25 days after sowing (DAS)), using the protocol described by Berger et al. (2012). Pots were watered every 2 days for up to 60 days after sowing to 16.8% (g/g) gravimetric water content, corresponding to 0.8 × field capacity. From day 61 after sowing, plants were watered daily to 16.8% (g/g) until seed set. For all samples, (50-100 mg) of leaf material was taken from the midpoint of the selected leaf blades at two time points. These comprised the 4th leaf blade after full emergence (15 days after salt treatment and 40 DAS) and from the flag leaf blade from the primary tiller at anthesis (62 days after salt treatment and 87 DAS). Samples were taken from plants growing in blocks 1, 3, 4, 6, and 8 ( Figure 1). This sampling strategy covered all varieties in all blocks. The nine barley varieties used exhibit very similar growth and development rates and so all reached both stages over the same time period. Thus, any epigenetic variation attributable to developmental or organ-to-organ variation was minimized. All leaf samples were immediately snap frozen in liquid nitrogen and stored at −80°C until DNA extraction. Whole plants were harvested at maturity and above-ground biomass was dried and weighed.

Greenhouse Environmental Conditions
The experiment was conducted in a 24-m 2 greenhouse (~8 m × 3 m), with a gable roof 4.5 m above the floor at the lowest and 6 m at the highest point. The greenhouse (34°58′16 S, 138°38′23 E) was oriented West-East ( Figure 1). To investigate the possible causes of position dependent variability of barley response across the greenhouse, environmental factors (temperature, relative humidity, and photosynthetic active rate) were recorded during the same period of the year (June 26, to October 12, 2015), using four sensor nodes located along the benches (Figure 1). Based on this period of the year, we deemed daytime to be between 7:00 AM and 6:00 PM. The sensor nodes were positioned 2 m apart and 1 m from the east and west walls ( Figure 1). Each node had a combination of sensors for photosynthetic active radiance (PAR) (model Quantum, LI-COR, Lincoln, Nebraska, USA) and for humidity/ temperature (Probe HMP60, Vaisala INTERCAP ® , Helsinki, Finland). Environmental data were recorded every minute for the duration of the experiment using wireless data loggers (National Instruments, Sydney, New South Wales, Australia). Before use for further analyses, recorded data were quality controlled to remove time slots when data were not present for all four nodes. To show the overall daily fluctuation of environmental factors between sensor nodes during the experiment, the average measure of each factor per hour was plotted for each node. Then, the vapor pressure deficit (VPD) for each time point was calculated according to Murray (1967): * 610:7 * 10 7:5T (237: 3+T) where RH = relative humidity, T = temperature, and the factor 610.7 × 10 7.5T/(237.3+T) = saturated vapor pressure (SVP).
Pairwise comparisons of each environmental factor at sensornode positions were performed using the Wilcoxon signed-rank test (Wilcoxon, 1945), on the R package "ggpubr" (Kassambara, 2019). These comparisons were performed independently for day and night periods.

DNA Extraction
Frozen plant material was homogenized in a bead beater (2010-Geno/Grinder, SPEX SamplePrep ® , USA) prior to DNA extraction using a Qiagen DNeasy kit according to the manufacturer's instructions. DNA samples were then quantified in a NanoDrop ® 1000 Spectrophotometer (V 3.8.1, ThermoFisher Scientific Inc., Australia) and concentrations were standardized to 10 ng/µl for subsequent MSAP analyses. To ensure marker reproducibility, DNA samples were analyzed in two technical replicates. Thus, samples were digested using a methylation insensitive restriction enzyme EcoRI in combination with either HpaII or MspI (isoschizomers), which show differential sensitivity to cytosine methylation at CCGG positions. Digested DNA fragments were ligated to adapters ( Table 1) with one end cohesive with restriction products generated by EcoRI or HpaII/MspI. Digestion and ligation reactions were performed in a single solution of 11 µl comprising: 1.1 µl T4 ligase buffer; 0.1 µl HpaII; 0.05 µl MspI; 0.25 µl EcoRI; 0.05 µl T4 ligase; 0.55 µl BSA ; 1.1 µl NaCl ; 1 µl Adapter EcoRI; 1 µl Adapter HpaII/MspI; 5.5 µl DNA sample and 0.3 µl pure water. Enzymes and buffer were acquired from New England Biolabs, Australia (NEB) and oligos were produced at Sigma-Aldrich, Australia. The solution was incubated for 2 h at 37°C, then enzymes were inactivated at 65°C for 10 min.

PCR
Two PCR amplifications were performed using products of the restriction/ligation reaction. First, a pre-amplification PCR was performed, in which primers complementary to adaptors but with 3' overhangs for a unique nucleotide (HpaII/MspI primer +C and EcoRI primer +A, Table 2) were used in a pre-optimized PCR master mix (BioMix ™ , Bioline, Meridian Bioscience; Australia) following the manufacturer's instructions. DNA digestion/ligation product (0.5 µl) was used for PCR amplification, with the following profile as per Rois et al. (2013): 72°C for 2 min, 29 cycles of 30 s denaturing at 94°C, 30 s annealing at 56°C and 2 min extension at 72°C, ending with 10 min at 72°C to ensure completion of the extension.
Pre-amplification products were quality assessed by 1% w/v agarose electrophoresis (80 V for 2 h), before performing the selective amplification using two selective primer combinations, EcoRI_AAG vs. HpaII/MspI_CCA and EcoRI-ATG vs. HpaII/ MspI_CAA. Amplified fragment detection through capillary electrophoresis was facilitated by labeling HpaII/MspI selective primers with the 6-FAM reporter molecule (6-carboxyfluorescein). Just 0.3 µl of pre-amplification product was used in the pre-optimized PCR master mix and the PCR was performed as follows (Rois et al., 2013); 94°C for 2 min, 12 cycles of 94°C for 30 s, 65°C (and decreasing by 0.7°C each cycle) for 30 s, and 72°C for 2 min, followed by 24 cycles of 94°C for 30 s, 56°C for 30 s, and 72°C for 2 min, ending with 72°C for 10 min.

Capillary Electrophoresis
The products of the selective PCR were fractionated by capillary electrophoresis on an ABI PRISM 3730 (Applied Biosystems, Foster City, California, USA) at the Australian Genome Research Facility Ltd (Adelaide, Australia). For this, 2 µl of selective PCR products were first combined with 15 µl of HiDi formamide (Applied Biosystems) and 0.5 µl of GeneScan ™ 500 ROX ™ Size Standard (Applied Biosystems). The mixture was then denatured at 95°C for 5 min and snap-cooled on ice for 5 min before sample fractionation at 15 kV for 6 s and at 15 kV for 33 min at 66°C.

MSAP Data Analysis
MSAP profiles obtained using HpaII and MspI were used to generate; 1) a qualitative binary matrix of allelic presence/ absence scores, and 2) a quantitative matrix of allelic peak height using GeneMapper Software v4 (Applied Biosystems). Qualitative epigenetic changes associated with greenhouse positional effect were analyzed using fragment sizes between 100 and 550 base pairs, which were selected to estimate epigenetic distance between individual plants (EpiGD) and subpopulations of plants (PhiPT) and perform Principal Coordinate Analyses (PCoA), using GenAlex 6.501 (Peakall and Smouse, 2012). Quantitative analysis of peak height was used to examine the effect of position on the methylation status of individual loci. We searched for MSAP markers that were differentially methylated between experimental blocks by comparing the fragment peak heights to survey for position effects on the plant methylation profile . Before differential methylation analysis, model-based normalization factors were calculated for the peak height libraries using the weighted trimmed mean method of Robinson and Oshlack (2010). For each variety and sampling method, peak heights were extracted and analyzed individually using the modeling approach of Mccarthy et al. (2012). To ensure the peak heights could be compared between positions, the individual models contained a term to account for variation between blocks as well as a term to capture the differences between the control and salt stress treatments. A likelihood ratio test was then performed to determine whether estimated coefficients for the positions were equal (Mccarthy et al., 2012). The p-values from these tests were then adjusted for multiple comparisons using the false discovery rate method of Benjamini and Hochberg (1995). Analyses were conducted using the R package edgeR , in the R statistical computing environment (R Core Team, 2019).
The extent of epigenetic divergence between salt treatments at the two developmental stages (4th leaf and anthesis) was assessed, first by performing a multiple correspondence analysis (MCA) on MSAP marker data. A linear discriminant analysis (LDA) was then performed on the MCA results. These analyses, referred to as MC-LDA thereafter, were done using the R packages FactoMineR and MASS (Lêet al., 2008;R Core Team, 2019). To visualize the results of comparisons involving more than two groups, the first two linear discriminant factors (LD1 and LD2) were plotted. Otherwise, a density plot of LD1 was performed.

Assessment of Correlations Between Epigenetic Profiles and Plant Phenotype
Epigenetic and phenotypic variability were estimated using averaged data per position for all nine barley varieties (Bishop et al., 2015). The software GraphPad Prism 6 v008 (Graph-Pad Software, San Diego, California, USA) was used to perform statistical analyses. Values of above-ground plant biomass were normalized by computing the ratio of plant biomass over the mean biomass for each individual experiencing the same treatment across all positions. The same formula was applied to grain yield. This normalization was intended to address quantitative variability between treatments and among barley genotypes. Then, biomass and yield distance matrices were generated using the difference between normalized values of any two individual plants.
We performed a Mantel Test (Mantel, 1967) to estimate the significance of the correlations between epigenetic distance and plant biomass, and position in the greenhouse. For this, we used matrices generated from epigenetic distance, physical distance and phenotypic (biomass or yield) differences estimated as described above. In all cases, the level of significance of the observed correlations was tested using 9,999 random permutations. Since both enzymes (HpaII, MspI) are methylation sensitive (Walder et al., 1983;Reyna-Loṕez et al., 1997), these enzymes can independently show epigenetic marks across the genome. Therefore, our inferences about plant epigenetic profile thereafter relate to results obtained using either enzyme or a combination of both.

Microclimatic Variability in the Greenhouse
Data quality control of climatic data provided 47,144 and 54,983 time-points of data recording for the periods of day and night, respectively. These correspond to time-points when recording was obtained simultaneously in all sensor nodes. There was clear evidence of both spatial and temporal variation for temperature, photosynthetically active radiation (PAR) and relative humidity (RH) within the experimental area (Figures 2 and 3).
The average dynamics of climatic data in the greenhouse showed a higher PAR between 8 AM and 10 AM at the East side than the rest of the greenhouse (node D, Figure 1). The PAR was also variable during the day between node positions, with sensor node B (Centre-West, Figure 1) recording the lowest PAR values around 12 PM (Figure 2A). The average temperatures evolved broadly in the same way at all node positions, with only around 1.5°C difference between the most divergent nodes at the warmest time of day ( Figure 2B). The RH was the highest at node A (West side of the greenhouse, Figure 1) during both day and night, and was significantly different from the rest of the positions during the day (Figures 2C and 3). The node D (East end of the greenhouse) presented the lowest RH during the day; it was not significantly different from nodes B and C ( Figure 3A).
Although there was no clear evidence of gradient between sensor nodes for any of the climatic factors (i.e. RH, temperature, VPD and PAR, the pairwise comparison of data from sensor nodes using Wilcoxon paired signed-rank test showed significant differences between positions for each variable (Figures 3A−G). Such differences were present during both day and night periods in the greenhouse. The RH appeared particularly variable at night between all positions of sensor nodes ( Figure 3B).

Correlation Between DNA Methylation Profile and Plant Position in the Greenhouse
As expected, the variation between MSAP profiles of the nine diverse varieties used in this study led to significant confounding clustering according to genotype. In subsequent analyses, we therefore elected either to consider perturbations to MSAP across all varieties collectively or else make comparisons on a variety by variety basis. The former included all confounding variation associated with genotype but sought to provide an indication of conserved effects across the panel. The latter analyses were intended to reveal the extent to which variability in the epigenetic response is influenced by genotype.
Plant DNA methylation profiles derived from MSAP data generated 269 alleles with sizes between 100 and 550 base pairs across samples from all nine barley varieties. PCoA of MSAP profiles for barley variety at anthesis showed grouping of samples more by plant position than salt treatment, regardless of the enzyme combination used (Figures 4A, B). The first coordinate Eigen space matched with the position of the plants in the greenhouse in the West-East direction (Figure 4). The Mantel test using all treatment samples together showed weak correlations between plant epigenetic profiles and plant positions in the greenhouse at 4th leaf stage, and more significant corrections at anthesis ( Table 3). For instance, for the variety Schooner, the Mantel test between pairwise epigenetic distance and plant position at the 4th leaf stage of barley development resulted in weak correlations for both HpaII (R 2 = 0.11, P-value = 0.025, Figure  5A) and MspI (R 2 = 0.12, P-value < 0.022, Figure 5C). Apart from two varieties (Buloke and Schooner), none of the remaining varieties showed a significant correlation between position and epigenetic profile at the 4th leaf stage ( Table 3, Figures S1).
Conversely, these correlations were stronger at anthesis for the same variety Schooner (R 2 = 0.48 and R 2 = 0.45, for HpaII and MspI, respectively, Figures 5B, D), with greater significance of the P-values (0.001). Additionally, all the remaining varieties showed significant correlation (P-value at least < 0.05) between DNA methylation profile at anthesis and the plant position in the greenhouse ( Table 3; Figure S1). The correlations at anthesis were high (R 2 > 0.3) for all varieties, except Buloke and Maritime ( Table 3).
The comparison of peak heights of MSAP markers generated from plants growing in different positions revealed significant differences between positions for some alleles ( Figure 6). In general, significant differences in peak height were observed between plants in position P1 and the other positions ( Figure 6). Overall, peak heights showed logarithmic trends (both positive and negative), significantly associated with the West-East distribution of the samples. A few markers were significantly different in peak heights over all positions (Table 4).
However, positional effect did not thwart the ability to differentiate between salt-stressed and control plants. The MC-LDA on MSAP marker data was able to separate salt stressed plants from those given control conditions ( Figures 7A, B). Furthermore, epigenetic divergence between treatment groups increased with time, with control and stress plants consistently more similar at the 4th leaf stage than at anthesis across all varieties ( Figures 7A, B and S2). MC-LDA of salt treatments could nevertheless discriminate treatments at both stages even though epigenetic divergence was strongly influenced by developmental stage (Figures 7C and S2).

Correlations Between Barley Phenotype, Epigenome, and Position
There was a clear trend in the final biomass of all nine barley varieties according to position, with a progressive increase from position P1 (west side of the greenhouse) to position P5 (East side) ( Figure 8A). This relationship was a logarithmic trend,  Figure 8B). However, when varieties were examined separately, both logarithmic and polynomial trends were observed ( Figure S3). Assessment of the relationship between pairwise differences in epigenetic distance and in grain yield showed significant correlations (P-values < 0.05) in control plants of six of nine varieties (Buloke, Commander, Fathom, Maritime, Schooner, Yarra), with R 2 varying between 0.247 and 0.907 (Table 5; Figure  S4). Likewise, stress plants showed significant correlations (P-values at least < 0.05) between grain yield and methylation profile in six varieties (Barque 73,Buloke,Commander,Flagship,Maritime,Schooner), with R 2 between 0.164 and 0.921 (Table 5; Figure  S4). An example of significant correlations between grain yield and epigenetic distance is presented in Figures 9A-D, for the variety Schooner.

Stochastic DNA Methylation Is Explained by Microclimatic Differences
The randomized block design aims to minimize unexplained variation between treatments, and has emerged as a preferred method in plant field trials and in controlled environment experiments (Edmondson, 1989;Guertal and Elkins, 1996;Brien et al., 2013). However, while block homogeneity is difficult to achieve, variability between blocks in the same experimental setting is often either ignored, attributed to randomness (Raj and Van ; temperature (C and D); vapour pressure deficit (VPD) (E, F); and photosynthetic active radiation (PAR) (G, PAR was deemed as null at night). Differences between positions were compared using Wilcoxon signed-rank test. *, ** and *** indicate the significance of the measured differences between positions for P value < 0.05, 0.01, and 0.001, respectively; ns, difference not significant. Oudenaarden, 2008; Karan et al., 2012;Tricker et al., 2012) or in the context of epigenetic research, explained by spontaneous occurrence of the methylation (Becker et al., 2011;Baulcombe and Dean, 2014;Van Der Graaf et al., 2015). In this study, we took care to control potentially confounding sources of variation between MSAP profiles by the selection of genetically diverse varieties with similar rates of growth and development, sourcing DNA from the same section of the same leaf from of all plants, and at two very distinct developmental stages . We nevertheless found evidence suggesting that microclimatic variation within a greenhouse was sufficient to trigger variability in the plant DNA methylation profile in a manner that was both independent of the experimental salt stress treatment and greater in magnitude. The clarity of the climatic variables measured across the experimental blocks, and the associated cline in methylation patterning is suggestive that each plant experienced a unique combination of climatic factors during the experimental period, and that this induces, at least partly, changes in methylation patterning. Similar observations were also reported for other greenhouse studies (Brien et al., 2013;Both et al., 2015;Cabrera-Bosquet et al., 2016). This finding is inconsistent with spontaneous DNA methylation being entirely responsible for the plant-plant variability in such experiments (Becker et al., 2011;Van Der Graaf et al., 2015), and throws into question how best to discriminate epigenetic responses to micro-environment fluctuations from those attributable to stochastic noise. Moreover, the effect of position can easily be overlooked in snap-shot exposure experiments, since the timeframe from stress exposure to induction of position-dependent methylation A B FIGURE 4 | Principal coordinates analysis (PCoA) of MSAP (methylation-sensitive amplified polymorphism) markers in barley variety Commander. MSAP markers were generated using five replicates of control (0 mM NaCl) and stress (75 mM NaCl) plant samples, for HpaII (A) and MspI (B). Positions 1 to 5 indicate experimental block numbers; symbols filled in black and hollow symbols represent salt stress (−S) and control (−C) samples, respectively. The PCoAs show sample distribution in the first two principal coordinates. Numbers in brackets represent the proportion of variation explained by the coordinate. Konate et al. Understanding Stochastic Epigenetic Variation Frontiers in Plant Science | www.frontiersin.org September 2020 | Volume 11 | Article 553907 markers is critical but also likely to vary between loci. Support for this reasoning can be taken from our findings that it was possible to separate salt and control samples by discriminate analysis at the 4th leaf stage and at anthesis but with higher divergence at the later stage. At the same time, correlation between epigenetic differences and physical distance among plants at anthesis (87 DAS) was stronger than at the 4th leaf stage (40 DAS), indicating that exposure to the stressor and positional microclimates both have a cumulative effect on the plant epigenome. These observations are congruent with the concept that plant adaptive improvisation, through DNA methylation, is proportional to the severity and duration of the environmental cue to which the plant was exposed (Soen et al., 2015). In this sense, the scale of the effect induced by intervention stress (salt) needs to be weighed against those imposed by coincidental stresses (microenvironment effects) but also by those associated with development or ageing, as was reported in humans (Gentilini et al., 2015). Any truly stochastic DNA methylation would represent residual variation. Previous studies have observed the influence of mesoclimatic conditions (Herrera and Bazaga, 2010) and factors such as temperature (Hashida et al., 2006), humidity (Tricker et al., 2012) or light (Barneche et al., 2014;Meyer, 2015) on methylome variability. However, the current study suggests, for the first time, that even slight variations in climatic factors (temperature, humidity or light) are sufficient to induce modifications in the plant DNA methylation profile, and that this can be sufficient to mask effects of mild stresses, Nine barley varieties were used, comprising ten individuals per variety, five replicates for control and stress plants. Samples were collected from the 4 th leaf (at 4 th leaf stage) and flag leaf (at anthesis). Epigenetic distances correspond to the Phi statistics of the MSAP markers between plant individuals. The coefficient of determination (R 2 ) was calculated according to Mantel (1967) using GenAlex 6.5. Asterisks (*), (**), (***) and (****) indicate significant correlation between treatments for P-value < 0.05, 0.01 and 0.001, respectively, estimated based on 9,999 permutations.

A B D C
FIGURE 5 | Correlation between pairwise epigenetic distance (Epi GD) and plant position in the greenhouse. The epigenetic distance was estimated at 4th leaf stage (A, C; 40 days after sowing) and anthesis (C, D; 87 days after sowing) of barley variety Schooner, using HpaII (A, B) and MspI (C, D) for the MSAP (methylation sensitive amplified polymorphism) analysis. Five replicates of control (0 mM NaCl) and stress (75 mM NaCl) were analyzed together and dots represent pairwise comparisons between individual plants. Equations represent the formula of the regression line, R 2 represents the coefficient of determination, calculated according to Mantel, 1967 using GenAlex 6.5. Asterisks (*) and (***) indicate significant correlation between treatments for P value < 0.05 and 0.001, respectively, estimated based on 9,999 permutations. as was observed here for salt stress. We certainly do not contend that all nascent methylation arises in response to environmental or biotic effectors but we do argue that far more care is needed before discounting unaccounted epigenetic variation as stochastic noise.

Positional Effect Affects Salt Stress-Induced DNA Methylation Changes in Barley
Positional effects in greenhouse experiments are well established and if not properly accounted for can generate uncharacterized background noise that can mask the effect of the experimental treatment (Edmondson, 1989;Guertal and Elkins, 1996;Brien et al., 2013). Spatial variability in coincident environmental factors has potential to introduce variability between replicate plants' development and response to experimental treatments (Edmondson, 1989;Guertal and Elkins, 1996). Such spatial variability is liable to introduce flaws in measurements and observations between replicates that, in fact, were not experiencing exactly the same constraints (Addelman, 1970). This can compromise the search for relationships between experimentally controlled stressors (in our study, soil salt stress) and perturbations in epigenetic profiles. Indeed, in the present work the observed FIGURE 6 | Exemplars of MSAP (methylation sensitive amplified polymorphism) alleles that show significant differences in peak height between positions in the greenhouse. Markers were detected in control (0 mM NaCl, red symbols) and stress (75 mM NaCl, blue symbols) plants; Vertical axis shows logarithm 2 (log 2) of peak height intensity and the horizontal axis represents positions in the greenhouse, in the west to east direction. The gray number in each plot represents −log10 of P values. The title of each plot shows the enzyme used (either HpaII (HPA) or MspI (MSP), the variety, and the allele identity number.  negative correlation between RH and differences in epigenetic differentiation between control and salt stressed pairs of plants growing in the different positions suggests that variations in environmental factors has interfered with reaction of the plant to mild salt stress. One possible mechanistic explanation is that the observed West to East decrease in RH changed the plant's requirement for water (Barnabaś et al., 2008;Verslues and Juenger, 2011), and this in turn may have affected the level of salt stress experienced by each plant. In this way, plants were grown under the same salt treatment but because they experienced different RH, are likely to exhibit different responses to the salt stress; hence the inconsistent salt-induced DNA methylation profiles.

Phenotypic Differences Associated to Greenhouse Microclimates Correlate With Epigenetic Differences
Plants have been long known to exhibit phenotypic symptoms of stress in organs that are not directly exposed to the stressor (Riley et al., 2002). Indeed, it is well-established that deficiency or toxicity of plant nutrients in the soil often becomes manifest as physical symptoms in the leaves in a wide range of plants, including barley (Uchida, 2000). Similar responses have been reported in the methylation profiles of DNA extracted from organs that are equally unconnected with the source of stress. For example, Konate et al. (2018) reported that exposure of the roots to mild soil salt stress impacts on the methylation profile of barley leaves. However, it is open to question is whether phenotypic symptoms of stress co-correlate with the Epigenetic distance between plants was calculated based on MSAP data generated using HpaII and MspI. Coefficients of determination (R 2 ) were computed according to Mantel (1967) using five replicates for each treatment per variety. Asterisks (*) and (**) indicate significant correlation between treatments for P value < 0.05, and 0.01, respectively, estimated based on 9,999 permutations.
A B D C FIGURE 9 | Correlation between pairwise epigenetic distance (EpiGD) and pairwise difference in grain yield between plants of the variety Schooner. The correlation was tested according to Mantel, 1967 using GenAlex 6.5. Epigenetic distance between plants was calculated based on MSAP (methylation sensitive amplified polymorphism) data generated using HpaII (A, B) and MspI (C, D). Pairwise differences in grain yield between plants were calculated separately for control (A, C) and stress (B, D) plants. Values of grain yield were normalized by computing the ratio of each individual plant grain yield over the mean grain yield for the same treatment across all positions. The dots represent pairwise comparisons between individual plants; equations represent the formulae of the regression line; R 2 represents the coefficient of determination of the Mantel test; asterisk (*) and (**) indicate significant correlation between treatments for P value < 0.05, and 0.01, respectively, estimated based on 9,999 permutations. epigenomic symptoms. The finding here of a plastic response by barley plants in terms of biomass and grain yield to subtle differences associated with greenhouse position corroborates earlier work by Lacaze et al. (2008) who suggested that barley is responsive to fluctuations in ambient conditions. We postulate that the irregularity of phenotypic variability patterns across barley varieties and treatments may have emerged from two complementary factors; 1) the genetic variability among barley varieties leading to differential responsiveness to positional effect, as reported elsewhere (Lacaze et al., 2008;Kren et al., 2015), and 2) the randomness of spatial microclimatic conditions, which did not have a linear spatial gradient. The influence of a genotype-by-environment effect on plant phenotype was expected (Gianoli and Palacio-López, 2009;Aspinwall et al., 2015), but the scale of phenotypic variation induced by small-scale environmental variation was not. Our findings highlight the possibility for plants to show substantial phenotypic responses to even slight variations in ambient conditions, and that homogeneity in temperature control does not have over-riding importance. Furthermore, our discovery of a significant correlation between barley MSAP profiles and grain yield suggests that DNA methylation could at least reflect and possibly contribute toward the plastic variation in plant phenotypes. These results are in accordance with a mounting body of evidence that plant plasticity is at least partly epigenetically governed (Boyko and Kovalchuk, 2008;Rois et al., 2013;Baulcombe and Dean, 2014;Aspinwall et al., 2015). Considered together, our results demonstrate a tight interplay between plant epigenome, environment and phenotype.

CONCLUSIONS
Homogeneity of environmental conditions is practically difficult to obtain in a greenhouse (Edmondson, 1989;Guertal and Elkins, 1996;Brien et al., 2013). Awareness of plant sensitivity to microclimate is therefore important, especially in epigenetic studies, where plant epigenomes seem to be extremely responsive to small fluctuations in environmental factors. This study reveals that at least some of the DNA methylation previously considered stochastic is likely to have been, at least partially, induced by 1) positional effects on growth conditions, 2) differences in the length of plant exposure to relatively trivial variations in environment and 3) synergistic effects of stress treatment (mild salt stress in this case) and microclimatic conditions. The correlation between phenotypic DNA methylation differentiations between plants grown in different microclimates suggests that position-induced DNA methylation, previously ignored or considered as stochastic, may be a substantial source of phenotypic variability. Accordingly, we advocate that future epigenetic analyses should take into account the effect of micro-variations in environmental factors by careful experimental design and by considering position-induced DNA methylation markers as strong candidates for finely-tuned response to small environmental changes. We also propose that further research is needed to untangle microclimate-induced epigenetic variations from epigenome instability due to experimental treatment and developmental stage. We also feel the possibility of a transgenerational transmission of these effects warrants urgent attention.

DATA AVAILABILITY STATEMENT
All datasets presented in this study are included in the article/ Supplementary Material.

AUTHOR CONTRIBUTIONS
MK performed the experiments, analyzed the data and wrote the manuscript. JT performed the statistical analysis of MSAP peak heights. MW, ES, BB, and CR conceived the experiments and supervised the work. All authors contributed to the article and approved the submitted version.