Main {#s1}
====

To better understand how genotype and environment interact to affect DNA methylation and transcription, we grew 150 *Arabidopsis thaliana* accessions from Sweden ([@bib29]) in two different environments, 10°C and 16°C, chosen because they lead to very different flowering behavior ([@bib4]). Relying on existing genome sequence information ([@bib29]), methylome- and transcriptome-sequencing data were generated (see 'Materials and methods').

In plants, DNA methylation occurs on cytosines in the CG, CHG, and CHH contexts (where H is any nucleotide except for C), each of which is catalyzed by independent pathways ([@bib12]; [@bib39]). Consistent with previous results ([@bib42]; [@bib11]; [@bib32]; [@bib25]; [@bib35]; [@bib14]) we found considerable variation between accessions regardless of context, even at the level of genome-wide averages ([Figure 1A](#fig1){ref-type="fig"}). Temperature, on the other hand, did not appear to affect genome-wide CG or CHG methylation, but had a significant effect on CHH methylation, levels of which were 14% higher at 16°C than at 10°C, on average ([Figure 1A](#fig1){ref-type="fig"}). To investigate the genetic basis of DNA methylation, we performed genome-wide association studies (GWAS) using different facets of average methylation as the phenotype. For global CG and CHG methylation, no associations reached genome-wide significance, while for CHH methylation a clear peak of association was observed on chromosome 4 ([Figure 1---figure supplement 1](#fig1s1){ref-type="fig"}). The association was even more significant when restricting attention to average CHH methylation of large transposons ([Figure 1B](#fig1){ref-type="fig"}), in agreement with the notion that this type of methylation mostly occurs in transposons in *Arabidopsis* ([@bib40]).10.7554/eLife.05255.003Figure 1.The effect of CMT2 on genome-wide CHH methylation levels.(**A**) Genome-wide average methylation level reaction norms for each accession (156 samples at 10°C and 125 samples at 16°C). Only CHH levels differ significantly between temperatures (Wilcoxon rank sum test; p = 1.7e-16). (**B**) Manhattan plot of genome-wide association studies (GWAS) results using average levels of CHH methylation for 151 accessions at 10°C on large transposons as the phenotype (the peak is also seen at 16°C \[not shown\]). The threshold line indicates a Bonferroni-corrected p-value of 0.05. (**C**) CHH methylation on large (over 2 kb) transposons at 10°C by CMT2 two-locus genotype (population sizes are 36, 82, and 24 for CMT2a^nr/nr^/CMT2b^r/r^, CMT2a^r/r^/CMT2b^r/r^, CMT2a^r/r^/CMT2b^nr/nr^, respectively). The values plotted are the Best Linear Unbiased Predictor (BLUP) estimates after correcting for population structure. Since accessions are homozygous, only four genotypes are possible, of which only three exist due to complete linkage disequilibrium between *CMT2a* and *CMT2b*. [Figure 1---figure supplement 1](#fig1s1){ref-type="fig"} shows Manhattan plots of GWAS results for global methylation averages. [Figure 1---figure supplement 2](#fig1s2){ref-type="fig"} shows Stepwise GWAS using average CHH methylation of TE\'s.**DOI:** [http://dx.doi.org/10.7554/eLife.05255.003](10.7554/eLife.05255.003)10.7554/eLife.05255.004Figure 1---figure supplement 1.Manhattan plots of GWAS results for global methylation averages.(**A**) CG methylation at 10°C. (**B**) CHG methylation at 10°C. (**C**) CHH methylation at 10°C. Results for methylation at 16°C were similar (data not shown).**DOI:** [http://dx.doi.org/10.7554/eLife.05255.004](10.7554/eLife.05255.004)10.7554/eLife.05255.005Figure 1---figure supplement 2.Stepwise GWAS using average CHH methylation of TE\'s as a phenotype.(**A**) Without a cofactor. (**B**) Including SNP on chr 4 at position 10,459,127 (*CMT2a*) as a cofactor. (**C**) Including snps on chr 4 at 10,459,127 (*CMT2a*) and 10,454,628 (*CMT2b*) as cofactors. The threshold line indicates a Bonferroni-corrected p-value of 0.05.**DOI:** [http://dx.doi.org/10.7554/eLife.05255.005](10.7554/eLife.05255.005)

The association centered around a SNP at 10,459,127 on chromosome 4, 38 kb downstream from the locus AT4G19020, which encodes a homolog of the CHG methyltransferase chromo-methylase-3 ([@bib26]) that has recently been shown to catalyze both CHH and CHG methylation on transposons, and is thus an excellent candidate ([@bib45]; [@bib39]). A multi-locus mixed model ([@bib34]) that included the identified SNP (*CMT2a*) as a fixed effect revealed another SNP downstream of *CMT2*, at position 10,454,628 (*CMT2b*), 4.5 kb closer to *CMT2* than *CMT2a*, and in complete linkage disequilibrium with it (i.e., the non-reference alleles at *CMT2a* and *CMT2b* are never seen together). Repeating the GWAS with both *CMT2a* and *CMT2b* as cofactors identified no further loci ([Figure 1---figure supplement 2](#fig1s2){ref-type="fig"}). Both non-reference alleles are common in southern Sweden, but are also found in the north (22.6% vs 9.5% and 30.6% vs 7.9% for *CMT2a* and *CMT2b*, respectively). Accessions with the non-reference *CMT2a* allele have on average more CHH methylation on transposons than those with the reference haplotype (p = 1.1e-03), while those with the non-reference *CMT2b* allele have lower levels of CHH methylation than the reference haplotype (p = 8.1e-03; [Figure 1C](#fig1){ref-type="fig"}). The associations were readily confirmed using an F2 population generated by crossing accessions with the *CMT2a* and *CMT2b* non-reference alleles ([Figure 2](#fig2){ref-type="fig"}). No significant differences in CMT2 mRNA levels were observed between the alleles in our data and limited Sanger sequencing of cDNA showed no evidence of splicing variants (although, as will be discussed below, we did detect a putative rare null allele). Several non-synonymous polymorphisms in the methyltransferase and BAH domains of CMT2 were detected ([Supplementary files 1](#SD1-data){ref-type="supplementary-material"} and [2](#SD2-data){ref-type="supplementary-material"}) but they do not explain the phenotype as well as the *CMT2a* and *CMT2b* SNPs.10.7554/eLife.05255.006Figure 2.CHH methylation levels in an F2 population map to CMT2.(**A**) CHH methylation on large transposons by CMT2 genotype in an F2 population of 113 individuals (population sizes are 19, 52, and 38 for CMT2a^nr/nr^/CMT2b^r/r^, CMT2a^r/r^/CMT2b^r/r^, CMT2a^r/r^/CMT2b^nr/nr^, respectively; 4 individuals whose genotype at CMT2 could not be accurately inferred were omitted). (**B**) Mapping of CHH methylation of long TEs in the same population. The dotted line indicates a LOD threshold with a genome-wide p-value of 0.05 obtained using 1000 permutations, and the vertical blue line shows the marker interval that contains CMT2.**DOI:** [http://dx.doi.org/10.7554/eLife.05255.006](10.7554/eLife.05255.006)

The effect of genetic variation on local CHH methylation was examined by calculating the methylation level in 200 bp sliding windows across the genome (100 bp overlap between windows) and running GWAS for the 200,000 differentially methylated regions (DMRs; see 'Materials and methods') that varied most between individuals. 36023 DMRs had at least one genome-wide significant association (Bonferroni-corrected p-value of 0.05; 7273 remain after correcting for multiple GWAS using an FDR of 0.05). 45% (15,031) of the DMRs had a significant *cis-*association within 100 kb, while the rest showed evidence of *trans*-regulation, including the dramatic effect of *CMT2* on chromosome 4 which accounted for approximately 21% (7392) of all significant associations ([Figure 3A](#fig3){ref-type="fig"}).10.7554/eLife.05255.007Figure 3.Genetic basis CHH methylation variation.(**A**) GWAS for CHH differentially methylated regions (DMRs) at 10°C in 151 accessions, defined using 200 bp sliding windows across the genome and selecting the 200,000 most variable ones. For each DMR, SNPs significantly associated at the Bonferroni-corrected 0.05-level are plotted. (**B**) Variance-components analysis of the CHH DMRs. For each DMR, a mixed model with *cis*, *CMT2*, and genome-wide *trans* effects, plus environment and genetic interactions with environment was fitted (see 'Materials and methods'). DMRs were binned by the total variance explained by the model. The density of DMRs in each bin is shown at the top, and the bottom shows the average variance-decomposition for each bin.**DOI:** [http://dx.doi.org/10.7554/eLife.05255.007](10.7554/eLife.05255.007)

To quantify the regulation of DMRs, we partitioned the variance of CHH methylation into environmental (E), *CMT2*, *CMT2* X E, *cis*, *cis* X E, *trans*, and *trans* X E using a mixed model ([Figure 3B](#fig3){ref-type="fig"}). The analysis confirmed substantial *cis* and *trans* associations, with the environment modulating the genetic effects rather than having a major direct effect. At least for the *cis* associations, a possible explanation is that SNPs tag polymorphic TE insertions, with the insertion allele being silenced in a temperature-sensitive manner.

The effect of temperature on CHH methylation could also be seen at the local level. We defined 'temperature DMRs' by looking for windows that differed significantly between temperatures. Comparing 16°C--10°C, each accession on average gained CHH methylation at ∼400 temperature DMRs and lost it at ∼200 temperature DMRs (false discovery rate = 0.05). CHH methylation is associated with transposable elements (TEs; [@bib12]), and in agreement with this, 79% of the temperature DMRs where methylation was gained at 16°C were located within 500 bp of an annotated TE (with 60% directly overlapping one). These temperature DMRs were enriched in a small subset of TEs (835, or 2.7% of total, permutation based p-value = 0.05) that were more highly methylated than other transposons, but with lower methylation levels immediately adjacent ([Figure 4A](#fig4){ref-type="fig"}). Compared to TEs without temperature DMRs, these 'variable' TEs also tended to be euchromatic ([Figure 4B](#fig4){ref-type="fig"}), highly expressed ([Figure 4C](#fig4){ref-type="fig"}), and recently inserted ('evolutionarily young' TE insertions for which orthologs are not present in *Arabidopsis lyrata* \[[@bib46]\] comprised 75% of the variable TEs vs 68% of non-variable TEs). At the super-family level, members of the *SINE*, *SINE*-like, *Helitrons* and *Mutator*-like DNA TE superfamilies were over-represented among the variable transposons, and at the family level, 36 families were over-represented, including the *AtREP*, *Vandal* and *HAT* DNA transposons, as well as *COPIA78/ONSEN* and *META1* retroelements ([Table 1](#tbl1){ref-type="table"}). Interestingly, *COPIA78* has been shown to become active in response to heat stress ([@bib30]; [@bib19]) apparently due to heat-shock promoter elements in its LTR regions ([@bib7]).10.7554/eLife.05255.008Figure 4.CHH methylation varies with temperature.(**A**) Average methylation levels over variable transposons at 10°C (orange) vs 16°C (red), and over non-variable transposons at 10°C (purple) vs 16°C (dark blue). Methylation for variable TEs is significantly higher (permutation p-value for CHH methylation = 0.05). (**B**) The density of variable (red) and non-variable TEs along chromosomes in 500 kb windows. Density is defined as the percentage of the total number in either category in each window; pericentromeric regions are shaded grey. (**C**) The expression of TEs at both temperatures. Variable TEs are more highly expressed than non-variable TEs, but the difference is only statistically significant at 16°C (Wilcoxon: 10°C, p = 0.15; 16°C, p = 0.023).**DOI:** [http://dx.doi.org/10.7554/eLife.05255.008](10.7554/eLife.05255.008)10.7554/eLife.05255.009Table 1.Super-families (italics) and families that are over-represented among 'variable' TEs**DOI:** [http://dx.doi.org/10.7554/eLife.05255.009](10.7554/eLife.05255.009)TE (*super*-)familyExpectedObservedEnrichment95th Quantile*RathE1_cons5264.5610RathE3_cons293.236RathE2_cons152.524SINE372.007RC/Helitron3464441.28368DNA/MuDR1441841.27162*ATREP245312.078RP1_AT22711.595ATTIRX1C11211.493ATREP132289.876VANDAL221118.563SIMPLEHAT11117.344VANDAL2N11107.323ATREP82136.475VANDAL2176.083ATREP101105.934AT9NMU1175.813ATN9_11105.754SIMPLEHAT21115.634META13205.417ATDNAI27T9A3154.836ATREP2A3154.836ATCOPIA78034.672VANDAL18NA034.672RathE1_cons5264.5610VANDAL14034.482SIMPLEGUY13134.196ATDNATA1034.002TNAT2A143.933ATREP74163.648RathE3_cons293.236ATREP14143.183ATREP16143.183LIMPET1392.906ATREP64142.898ARNOLDY27222.8513ATSINE4272.675ATDNAI27T9C272.426ATREP338922.3949ARNOLDY16142.2111ATREP113231.7319HELITRONY337511.3648

In order to gain further insight into the mechanisms responsible for variation in CHH methylation, we bisulfite-sequenced knockout lines of CMT2 (SAIL_906_G03) and DCL3 (*dcl3-5* \[[@bib9]\], a component of the RdDM pathway), and identified 10,138 DCL3-dependant DMRs and 33,422 CMT2-dependent DMRs as described in section 'DMR calling on DNA methylation mutants' of the 'Materials and methods'. As expected under the assumption that CMT2 is responsible for the massive GWAS peak on chromosome 4, the GWAS peak at this locus remains if we consider only the CMT2-dependent DMRs, but not for DCL3-dependent DMRs ([Figure 5---figure supplement 1](#fig5s1){ref-type="fig"}). Furthermore, while CHH methylation varied with temperature at both DCL3- and CMT2-dependent DMRs ([Figure 5](#fig5){ref-type="fig"}), DCL3-dependent DMRs were much more strongly associated with previously identified temperature DMRs (4703 out of 10,138 DCL3-dependent DMRs, or 46%, overlapped temperature-sensitive DMRs, whereas the corresponding numbers for CMT2-dependent DMRs were 2299 out of 33,422, or 7%; Fisher\'s exact p-value \< 2.2e-16), suggesting that much of the temperature variation in CHH methylation is due to components of the RdDM pathway. This result is consistent with previous findings showing that RNA silencing is less active at lower temperatures ([@bib31]).10.7554/eLife.05255.010Figure 5.Temperature dependent CHH methylation variation at RdDM and CMT2 controlled DMRs.CHH methylation at CMT2- and DCL3-dependent DMRs in natural accessions grown at 10°C and 16°C (cf. [Figure 1A](#fig1){ref-type="fig"}, each population has 110 individuals). The difference between temperatures was highly significant for both types of DMR (Wilcoxon p-value = 9.1e-11 and p-value = 5.9e-12 respectively). Black points/grey lines indicate accessions with the *CMT2* reference allele; green, the *CMT2a* non-reference allele; and orange, the *CMT2b* non-reference allele. Red is the TAA-03 accession, which has a putative null allele of CMT2. Average methylation levels for each of the genotypes are shown in bars to the side [Figure 5---figure supplement 1](#fig5s1){ref-type="fig"} shows GWAS on CMT2 and DCL3 dependant DMRs. [Figure 5---figure supplement 2](#fig5s2){ref-type="fig"} shows a putative null allele of CMT2.**DOI:** [http://dx.doi.org/10.7554/eLife.05255.010](10.7554/eLife.05255.010)10.7554/eLife.05255.011Figure 5---figure supplement 1.GWAS on CMT2 and DCL3 dependent DMRs.(**A**) GWAS for CMT2-dependent DMRs at 10°C. (**B**) GWAS on DCL3-dependent DMRs at 10°C. Results from 16°C were similar in both cases. The threshold line indicates a Bonferoni-corrected p-value of 0.05.**DOI:** [http://dx.doi.org/10.7554/eLife.05255.011](10.7554/eLife.05255.011)10.7554/eLife.05255.012Figure 5---figure supplement 2.Putative null allele of CMT2.A screenshot from a genome browser indicating the lack of read coverage for CMT2 stretching from intron 7 to exon 16 in the accession TAA-03.**DOI:** [http://dx.doi.org/10.7554/eLife.05255.012](10.7554/eLife.05255.012)

Interestingly, we observed one accession from northern Sweden, TAA-03, with almost undetectable levels of CHH methylation at CMT2-dependant DMRs ([Figure 5](#fig5){ref-type="fig"}). Further investigation suggested that it has a deletion or rearrangement in CMT2, as we were unable to map reads between positions 2813 and 4944 (intron 7 to exon 16, [Figure 5---figure supplement 2](#fig5s2){ref-type="fig"}). Sanger-sequencing indicates the insertion of a stretch of TC dinucleotide repeats of at least 330 bp. The same deletion appears to be present in three more accessions from northern Sweden (TAA-14, TAA-18, and Gro-3) a situation reminiscent of the homologous CMT1 gene, which seems to be non-functional in most *Arabidopsis* accessions ([@bib17]). Although CMT2 null alleles have no obvious phenotype, the gene is highly conserved in plants (with the exception of maize; [@bib45]; [@bib43]).

It has recently been suggested that natural variation in CMT2 is associated with adaptation to climate ([@bib36]), although the alleles identified in that study do not overlap with the ones identified here. Given the sensitivity of CHH methylation to growth temperature observed here, we next investigated the correlation between DNA methylation and the climate of origin ([@bib15]). While CHH methylation was moderately correlated with photosynthetically active radiation (PAR) in spring (Pearson\'s r = 0.38), and CHG showed correlation with aridity (r = 0.35) and the number of frost-free days (Pearson\'s r = 0.30), by far the strongest signal was a strong positive correlation between CG methylation and latitude (Pearson\'s r = 0.70), as well as with a number of environmental variables that co-vary with latitude in our sample, such as minimum temperature and daytime length ([Table 2](#tbl2){ref-type="table"}, [Figure 6A](#fig6){ref-type="fig"}). As a result of the strong latitudinal correlation, accessions originating from northern Sweden (minimum temperature below −10°C) had on average 11% higher global CG methylation compared to those from the south ([Figure 6A](#fig6){ref-type="fig"}). The correlation appears to be driven by gene body methylation (GBM): as the correlation for CG methylation on transposons was much weaker ([Figure 6A](#fig6){ref-type="fig"}, [Figure 6---figure supplement 1](#fig6s1){ref-type="fig"}). Because the methylation variation observed for genes with average CG methylation below 5% appeared mostly to be noise ([Figure 6---figure supplement 2](#fig6s2){ref-type="fig"}, see also the 'Materials and methods' section), we classified genes into 'unmethylated' and 'having GBM' using this as a cutoff (5%). We also eliminated genes showing a transposon-like pattern of methylation in which not only CG, but also the CHH and CHG contexts are highly methylated ([@bib45]). In what follows, we use GBM to refer only to gene body CG methylation for this filtered set. In order to better understand the observed variation in GBM, we examined CG methylation at the single nucleotide resolution within GBM containing genes. Although methylation was detectable (using a cut-off of 1%) at a similar number of sites in the north and south (1085292 vs 1079443 CG sites), those in the north showed a distinct skew towards higher methylation levels ([Figure 6---figure supplement 3](#fig6s3){ref-type="fig"}). Likewise when the difference between average methylation levels in the north and south was calculated individually for each of the cytosines, the majority of cytosines showed a small increase in the north compared to the south ([Figure 6---figure supplement 4](#fig6s4){ref-type="fig"}). From this we concluded that there is a general small increase in methylation of most CG dinucleotides in GBM genes, rather than large changes in a specific subset. GBM primarily occurs on long, evolutionary conserved genes that tend to be moderately-to-highly expressed, and is positively correlated with gene expression ([@bib48]; [@bib41]). Genes with higher GBM tended to be more highly expressed in our data as well, and---more interestingly---accessions with higher average GBM showed slightly higher average expression of methylated genes (although the correlation was weak, [Figure 6---figure supplement 5](#fig6s5){ref-type="fig"}). Given that northern accessions had higher GBM, this meant that genes with GBM were on average more highly expressed in northern than in southern accessions, while unmethylated genes showed little difference ([Figure 6B](#fig6){ref-type="fig"}). GBM has previously been shown to be anti-correlated with temperature-dependent gene expression ([@bib23]). While no large-scale north-south expression differences were observed between 10°C and 16°C in our data, northern accessions showed considerably less variation in expression between the two temperatures for genes with GBM (Wilcoxon p-value = 1.2e-05), while no such difference was observed for genes without it ([Figure 6---figure supplement 6](#fig6s6){ref-type="fig"}).10.7554/eLife.05255.013Table 2.Correlation between methylation levels and environment-of-origin variables ([@bib15])**DOI:** [http://dx.doi.org/10.7554/eLife.05255.013](10.7554/eLife.05255.013)Environmental variableGrowing temp.CGCHGCHHrrhop-valuerrhop-valuerrhop-valueLatitude100.690.527.8E-11−0.24−0.192.7E-020.100.141.1E-01160.620.473.2E-07−0.21−0.204.2E-020.04−0.112.5E-01Longitude100.590.541.2E-11−0.14−0.093.1E-010.230.287.5E-04160.550.534.4E-09−0.12−0.037.4E-010.140.151.2E-01Temperature seasonality100.680.491.6E-09−0.27−0.244.8E-030.090.092.8E-01160.620.421.1E-05−0.23−0.266.6E-030.04−0.122.1E-01Max. temp. (warmest month)10−0.140.064.6E-01−0.07−0.131.3E-010.140.202.0E-0216−0.030.102.9E-01−0.10−0.203.8E-020.050.037.3E-01Min. temp. (coldest month)10−0.70−0.569.1E-130.270.211.2E-02−0.07−0.064.7E-0116−0.63−0.482.7E-070.240.241.4E-020.000.195.6E-02Precipitation (wettest month)100.450.521.2E-10−0.25−0.271.2E-03−0.20−0.121.7E-01160.290.434.0E-06−0.26−0.241.2E-02−0.22−0.195.8E-02Precipitation (driest month)100.310.401.5E-06−0.33−0.296.5E-04−0.24−0.211.6E-02160.210.327.4E-04−0.26−0.241.4E-02−0.15−0.186.0E-02Precipitation seasonality100.420.447.1E-08−0.07−0.165.4E-020.050.019.0E-01160.360.371.2E-04−0.13−0.161.1E-010.01−0.019.1E-01PAR (spring)100.040.228.9E-030.200.183.7E-020.240.237.3E-03160.030.186.6E-020.270.213.5E-020.380.352.8E-04Length of growing season10−0.59−0.575.5E-130.240.237.3E-03−0.16−0.183.3E-0216−0.58−0.544.0E-090.230.213.0E-02−0.040.018.9E-01No. consecutive cold days100.600.534.0E-11−0.19−0.131.2E-010.170.281.1E-03160.570.534.2E-09−0.17−0.093.7E-010.100.084.1E-01No. consecutive frost-free days10−0.59−0.491.2E-090.290.271.5E-030.020.037.1E-0116−0.51−0.394.9E-050.300.301.6E-030.070.131.9E-01Relative humidity (spring)100.620.475.6E-09−0.23−0.183.9E-020.090.064.5E-01160.530.371.2E-04−0.20−0.267.6E-030.04−0.084.3E-01Daylength (spring)100.690.507.2E-10−0.27−0.211.4E-020.080.055.7E-01160.630.411.5E-05−0.23−0.292.7E-030.04−0.178.7E-02Aridity100.530.498.4E-10−0.35−0.311.9E-04−0.18−0.211.3E-02160.430.428.4E-06−0.28−0.241.3E-02−0.13−0.203.8E-02[^2][^3]10.7554/eLife.05255.014Figure 6.Latitudinal difference in gene body methylation (GBM) and gene expression.(**A**) Global CG methylation levels at 10°C for 151 accessions are strongly correlated with minimum temperature at the location of origin. Results for 16°C are similar. (**B**) Genes with GBM are more highly expressed at 10°C in northern (blue) than in southern (red) accessions (wilcoxon rank sum test p = 2.1e-03), whereas genes without GBM show little difference (p = 1.9e-02). At 16°C the difference for genes with GBM is more significant (p = 6.4e-05), whereas the difference for genes without GBM is insignificant (p = 0.49).**DOI:** [http://dx.doi.org/10.7554/eLife.05255.014](10.7554/eLife.05255.014)10.7554/eLife.05255.015Figure 6---figure supplement 1.Correlation between CG methylation levels and the minimum temperature at location of origin.Above, GBM at 10°C and 16°C. Below, TE CG methylation at 10°C and 16°C.**DOI:** [http://dx.doi.org/10.7554/eLife.05255.015](10.7554/eLife.05255.015)10.7554/eLife.05255.016Figure 6---figure supplement 2.Filtering of GBM variation data.(**A**) Genes with low or no CHG methylation have variable levels of CG methylation, while genes with appreciable CHG methylation have very high CG (and CHH) methylation. (**B**) Among genes with only CG GBM, variance-component analysis reveals a bimodal distribution of the total variance explained: variation in methylation for genes with low levels of methylation typically does not appear to have a genetic basis.**DOI:** [http://dx.doi.org/10.7554/eLife.05255.016](10.7554/eLife.05255.016)10.7554/eLife.05255.017Figure 6---figure supplement 3.Distribution of methylation levels at individual CG dinucleotides within GBM genes.The histogram shows the average methylation level for each individual CG dinucleotide on GBM genes in all accessions in the north (blue) or in the south (red).**DOI:** [http://dx.doi.org/10.7554/eLife.05255.017](10.7554/eLife.05255.017)10.7554/eLife.05255.018Figure 6---figure supplement 4.Distribution of variation in methylation levels between the north and the south for individual CG dinucleotides within GBM genes.The histogram shows the average methylation level for each individual CG dinucleotide on GBM genes in all accessions in the north minus the average methylation level in the south for each dinucleotide.**DOI:** [http://dx.doi.org/10.7554/eLife.05255.018](10.7554/eLife.05255.018)10.7554/eLife.05255.019Figure 6---figure supplement 5.Accessions with higher average GBM tend to have higher average expression (of genes with GBM, normalized by genes without GBM; r = 0.131, p = 0.0386).**DOI:** [http://dx.doi.org/10.7554/eLife.05255.019](10.7554/eLife.05255.019)10.7554/eLife.05255.020Figure 6---figure supplement 6.Genes with GBM show less expression variation between temperatures.Mean per-gene variation in expression between 10°C and 16°C is reduced for GBM containing genes in northern (blue) accessions compared to southern (red) accessions (wilcoxon rank sum test p = 1.2e-05), whereas for genes without GBM the difference between north (light blue) and south (pink) is insignificant.**DOI:** [http://dx.doi.org/10.7554/eLife.05255.020](10.7554/eLife.05255.020)

As for CHH DMRs, the genetic basis of GBM was examined using a variance-component approach ([Figure 7A](#fig7){ref-type="fig"}). The results were dramatically different: relative to CHH methylation, the *trans* effects for GBM are massive, and the environment appears to have no effect (in agreement with the observation that only CHH methylation levels vary with temperature, see [Figure 1A](#fig1){ref-type="fig"}). To identify the genes responsible, we also performed GWAS for each gene with GBM ([Figure 7B](#fig7){ref-type="fig"}). A total of 3241 significant associations were found for 2315 genes. 43% of these genes had a significant *cis*-association (within 100 kb of the gene of interest), which could represent local variants affecting methylation directly, or indirectly by affecting gene expression ([@bib13]). No evidence for major *trans*-acting loci like *CMT2* was found, but 69% of all significant associations were in *trans*. A comparison of the direction of the effect of GBM-associated SNPs in *cis* and *trans* revealed a striking pattern ([Figure 7C](#fig7){ref-type="fig"}). While the non-reference alleles of *cis*-SNPs were 1.18 times more likely to be associated with decreased rather than increased GBM (p = 2.01e-04), the non-reference alleles of *trans*-SNPs were 3.48 times more likely to be associated with increased GBM (p = 2.2e-16), and the non-reference alleles at the 15 *trans*-SNPs that were associated with GBM at five or more genes were always positively correlated ([Figure 7C](#fig7){ref-type="fig"}). Furthermore, while *cis*-SNPs showed a wide distribution of allele frequencies similar to random SNPs, *trans*-SNPs showed a much more limited distribution of frequencies ([Figure 8A](#fig8){ref-type="fig"}) and were also much more strongly correlated with latitude ([Figure 8B,C](#fig8){ref-type="fig"}). The correlation between GBM and latitude thus appears mostly to be due to *trans*-acting SNPs.10.7554/eLife.05255.021Figure 7.The genetic basis of GBM.(**A**) Variance component analysis of GBM. (**B**) Significant associations (Bonferroni-corrected 0.05-level) from a GWAS of GBM for individual genes. (**C**) Correlation between non-reference allele at associated SNPs and GBM.**DOI:** [http://dx.doi.org/10.7554/eLife.05255.021](10.7554/eLife.05255.021)10.7554/eLife.05255.022Figure 8.Frequency and distribution of GBM associated SNPs.(**A**) Correlation between non-reference allele at associated SNPs and latitude. (**B**) Non-reference allele frequency distribution for *cis* and *trans* SNPs compared to random SNPs. (**C**) Accessions carrying the non-reference alleles are limited to northern Sweden (accessions with the non-reference allele at 8 or more of the 15 loci blue, remaining accessions are red).**DOI:** [http://dx.doi.org/10.7554/eLife.05255.022](10.7554/eLife.05255.022)10.7554/eLife.05255.023Figure 8---figure supplement 1.Linkage disequilibrium between the 15 highly associated trans-SNPs.**DOI:** [http://dx.doi.org/10.7554/eLife.05255.023](10.7554/eLife.05255.023)

The 15 highly associated *trans*-SNPs were largely limited to northern Sweden, and were in strong linkage disequilibrium with each other ([Figure 8---figure supplement 1](#fig8s1){ref-type="fig"}). *A. thaliana* from northern Sweden show signs of multiple strong selective sweeps ([@bib29]) and harbors many polymorphisms that appear to be involved in local adaptation (specifically to minimum temperature; [@bib15]). The 15 SNPs were more than ninefold over-represented in the previously identified sweep regions (empirical p-value = 5.1e-03) and over fivefold over-represented within 2 kb of SNPs in the 1% tail of those associated with minimum temperature ([@bib15]) (empirical p-value = 3.1e-04), ([Table 3](#tbl3){ref-type="table"}). The ancestral state could be determined for 10 of the 15 SNPs, and in 8 of these cases, the non-reference allele was derived, further supporting sweeps in northern Sweden.10.7554/eLife.05255.024Table 3.15 SNPs associated with gene body methylation (GBM) at 5 or more genes**DOI:** [http://dx.doi.org/10.7554/eLife.05255.024](10.7554/eLife.05255.024)ChrPositionAssociated with GBM at how many genes?Non-reference allele countSNP-latitude correlationOverlap with sweep ([@bib29])Overlap with min. temp. Assoc. SNPs ([@bib15])19122918420.73none1_914088_0.21144051035660.64nonenone176141015480.66nonenone1197559675880.75none1_19757140_0.24269986316550.872_6931030none276550166810.612_7613651none276604699550.782_76136512_7662427_0.30276660595690.722_76136512_7665507_0.25276808825820.632_7613651none279157126510.83none2_7913782_0.23293824955730.71none2_9383856_0.34296538789480.80nonenone34193098660.68nonenone45199828660.70nonenone4132900345740.74nonenone

That the difference in GBM between north and south is likely to reflect local adaptation is also clear from its relative magnitude. The north vs south divide explains a much higher fraction of the additive genetic variance for GBM (Qst = 0.772; see 'Materials and methods') than of the SNP variance (Fst = 0.187). This strongly suggest that the phenotypic differentiation is driven by selection rather than genetic drift ([@bib24]).

Identifying the causal variants is challenging, a gene-ontology analysis of genes within 100 kb (the average size of the sweep regions, [@bib29]), of the 15 *trans*-SNPs found enrichment of loci associated with mRNA transcription (GO0009299, p-value = 2.62e-03). Genes involved in epigenetic processes are not captured well by standard gene-ontology, but we found that genes from the plant chromatin database ([www.chromdb.org/](http://www.chromdb.org/)) were significantly overrepresented in these regions as well (permutation p-value = 0.012; [Table 4](#tbl4){ref-type="table"}).10.7554/eLife.05255.025Table 4.Genes in the plant chromatin database that are within 100 kb of one of the 15 SNPs associated with GBM at 5 or more genes**DOI:** [http://dx.doi.org/10.7554/eLife.05255.025](10.7554/eLife.05255.025)ChromDBLocusARID3AT2G17410ARP3AT1G13180CHB4AT1G21700CHR9AT1G03750CHR35AT2G16390CONS3AT3G02380DNG12AT1G21710FLCP39AT3G02310FLCP16AT2G22630FLCP9AT2G22540GTI1AT2G22720HMGB4AT2G17560JMJ27AT4G00990NFA1AT4G26110SDG23AT2G22740SDG37AT2G17900YDG2AT2G18000HON3AT2G18050

We also looked for genes whose expression variation is consistent with a causal role. We identified 68 genes within 100 kB of one of the 15-*trans* SNPs whose expression is highly correlated (Wilcoxon test p \< 0.001) with the adjacent SNP after correction for population structure ([Table 5](#tbl5){ref-type="table"}). No significant enrichment of Gene Ontology terms was observed among these genes, and manual inspection identified no proteins directly involved in DNA methylation. Instead, a number of proteins involved in the regulation of gene expression and/or chromatin accessibility were present ([Table 5](#tbl5){ref-type="table"}). This may suggest that the increased gene-body methylation observed in the north is not directly due to increased DNA methylation, but may be caused by increases in gene expression driven either by differences in transcription factors networks or chromatin compaction. Identification of the causal variants behind this phenomenon should provide insight into how plants adapt to their local environment.10.7554/eLife.05255.026Table 5.Genes within 100 kb of the 15 SNPs associated with GBM at 5 or more genes whose expression is also correlated with the SNP**DOI:** [http://dx.doi.org/10.7554/eLife.05255.026](10.7554/eLife.05255.026)SNPLocusDesciptionp-value1_19755967AT1G53030Encodes a copper chaperone4.72E-071_19755967AT1G52880**NO APICAL MERISTEM (NAM) Transcription factor with a NAC domain**5.47E-071_19755967AT1G52990Thioredoxin family protein2.36E-051_19755967AT1G52780Protein of unknown function (DUF2921)1.46E-041_4405103AT1G12750RHOMBOID-like protein 6 (RBL6); FUNCTIONS IN: serine-type endopeptidase activity3.74E-081_4405103AT1G12790RuvA domain 2-like2.76E-051_4405103AT1G12730GPI transamidase subunit2.81E-051_4405103AT1G13080CYTOCHROME P450 FAMILY 71 SUBFAMILY B POLYPEPTIDE 2 (CYP71B2)1.65E-041_7614101AT1G21790TRAM LAG1 and CLN8 (TLC) lipid-sensing domain containing protein1.10E-051_7614101AT1G21900Encodes an ER-localized p24 protein8.81E-051_7614101AT1G21760**F-BOX PROTEIN 7 (FBP7) putative translation regulator in temperature stress response**8.54E-041_912291AT1G03660Ankyrin-repeat containing protein1.26E-101_912291AT1G03770**RING1B protein with similarity to polycomb repressive core complex1 (PRC1)**5.76E-071_912291AT1G03940HXXXD-type acyl-transferase family protein1.18E-061_912291AT1G03610Protein of unknown function (DUF789)6.91E-061_912291AT1G03580Pseudogene with weak similarity to ubiquitin-specific protease 121.29E-051_912291AT1G03830Guanylate-binding family protein3.50E-052_6998631AT2G16340Unknown protein1.35E-082_6998631AT2G16210**Transcriptional factor B3 family protein**1.69E-042_7666059AT2G17630Pyridoxal phosphate (PLP)-dependent transferases superfamily protein2.47E-182_7660469AT2G17620Cyclin B2;1 (CYCB2;1)9.68E-072_7655016AT2G17740Cysteine/Histidine-rich C1 domain family protein1.22E-042_7655016AT2G17420NADPH-DEPENDENT THIOREDOXIN REDUCTASE 2 (NTR2)9.96E-042_7666059AT2G17430MILDEW RESISTANCE LOCUS O 7 (MLO7)7.56E-042_7915712AT2G18100Protein of unknown function (DUF726)1.73E-062_7915712AT2G17980ATSLY member of SLY1 Gene Family1.33E-052_7915712AT2G18400Ribosomal protein L6 family protein1.26E-042_7915712AT2G18150Haem peroxidase8.05E-042_7915712AT2G18050**HISTONE H1-3 (HIS1-3)**9.47E-042_9382495AT2G22260HOMOLOG OF *E. COLI* ALKB (ALKBH2) enzyme involved in DNA methylation damage repair1.21E-082_9382495AT2G21850Cysteine/Histidine-rich C1 domain family protein5.38E-062_9382495AT2G22240MYO-INOSITOL-1-PHOSPHATE SYNTHASE 1 (MIPS1)8.71E-052_9382495AT2G21940SHIKIMATE KINASE 1 (ATSK1) localized to the chloroplast1.80E-042_9653878AT2G22660Protein of unknown function (duplicated DUF1399)2.22E-142_9653878AT2G22900Galactosyl transferase GMA12/MNN10 family protein5.08E-092_9653878AT2G22830Squalene epoxidase 2 (SQE2)3.91E-062_9653878AT2G22640BRICK1 (BRK1)6.17E-052_9653878AT2G22540**SHORT VEGETATIVE PHASE (SVP) Floral repressor involved in thermosensory pathway**2.46E-042_9653878AT2G22570NICOTINAMIDASE 1 (NIC1)2.67E-042_9653878AT2G22770**NAI1 Transcription factor**7.71E-043_419309AT3G02220Protein of unknown function (DUF2039)2.06E-163_419309AT3G02230REVERSIBLY GLYCOSYLATED POLYPEPTIDE 1 (RGP1)4.58E-143_419309AT3G02300Regulator of chromosome condensation (RCC1) family protein1.25E-103_419309AT3G02120Hydroxyproline-rich glycoprotein family protein1.81E-093_419309AT3G01980Short-chain dehydrogenase/reductase (SDR)3.91E-093_419309AT3G02370Unknown protein4.53E-083_419309AT3G02020ASPARTATE KINASE 3 (AK3)4.18E-073_419309AT3G02160**Bromodomain transcription factor**2.60E-063_419309AT3G02390Unknown chloroplast protein5.60E-063_419309AT3G02050K+ UPTAKE TRANSPORTER 3 (KUP3)1.28E-053_419309AT3G02125Unknown chloroplast protein2.12E-053_419309AT3G02200Proteasome component (PCI) domain protein1.16E-043_419309AT3G02180SPIRAL1-LIKE3 Regulates cortical microtubule organization4.56E-043_419309AT3G02250O-fucosyltransferase family protein5.31E-043_419309AT3G02110Serine carboxypeptidase-like 25 (scpl25)6.18E-044_13290034AT4G26255**Non-coding RNA of unknown function**1.67E-134_13290034AT4G26450WPP DOMAIN INTERACTING PROTEIN 1 (WIP1)1.13E-044_13290034AT4G26230Ribosomal protein L31e family protein1.74E-044_13290034AT4G26160ATYPICAL CYS HIS RICH THIOREDOXIN 1 (ACHT1)5.72E-044_519982AT4G01090Protein of unknown function (DUF3133)1.23E-064_519982AT4G01230Reticulon family protein2.33E-054_519982AT4G01410Late embryogenesis abundant (LEA) hydroxyproline-rich glycoprotein family5.44E-054_519982AT4G01330Serine/threonine-protein kinase2.22E-044_519982AT4G01200Calcium-dependent lipid-binding (CaLB domain) family protein3.93E-044_519982AT4G01390TRAF-like family protein3.99E-044_519982AT4G01040Glycosyl hydrolase superfamily protein5.66E-044_519982AT4G01000Ubiquitin-like superfamily protein8.55E-04

In conclusion, genes with GBM are generally up-regulated and more heavily methylated in northern accessions, and the change appears to be due to *trans*-acting polymorphisms that have been subject to directional selection. The candidate regions show an overrepresentation of genes involved in transcriptional processes. We also found that CHH methylation of TEs is temperature sensitive, and identified a major *trans*-acting controller, *CMT2*. Taken together, these observations suggest that local adaptation in *A. thaliana* involves genome-wide changes in fundamental mechanisms of gene regulation, perhaps as a form of temperature compensation.

Materials and methods {#s2}
=====================

Raw data generation {#s2-1}
-------------------

### Plant growth {#s2-1-1}

A diverse set of 150 Swedish accessions were sown on soil and stratified for 3 days at 4°C in the dark. They were then transferred to environmentally controlled growth chambers set at 10°C or 16°C under long day conditions (04:00--20:00) and individual seedlings were transplanted to single pots after 1 week. When plants attained the 9-true-leaf stage of development, whole rosettes were collected between 15:00 and 16:00 hr and flash frozen in liquid nitrogen.

In addition, a cross between the T550 and Brösarp-11-135 accessions was created and F2 seed obtained. 113 individual F2 lines were grown in the same manner as the accessions.

### RNA-seq library preparation {#s2-1-2}

For each accession, three plants were pooled and total RNA was extracted by TRIzol (Invitrogen, Carlsbad, California, 15596-018), DNase treated and mRNA purified with oligo dT Dynabeads (Life Technologies, Carlsbad, California). RNA was then fragmented using Ambion Fragmentation buffer (Life Technologies) and first and second strand cDNA synthesis was carried out using Invitrogen kit 18064-071. The ends of sheared fragments were repaired using Epicentre (Madison, Wisconsion) kit ER81050. After A-tailing using exo-Klenow fragment (New England Biolabs, Ipswich, Massachusetts, NEB M0212L), barcoded adaptors were ligated with Epicentre Fast-Link DNA Ligation Kit (Epicentre LK6201H). Adaptor-ligated DNA was resolved on 1.5% low melt agarose gels for 1 hr at 100 V. DNA in the range of 200--250 bp was excised from the gel and purified with the Zymoclean Gel DNA recovery kit (Zymo Research). The libraries were amplified by PCR for 15 cycles with Illumina PCR primers 1.1 and 1.2 with Phusion polymerase (NEB F-530L).

Single-end 32 bp sequencing was performed at the University of Southern California Epigenome Center on an Illumina (San Diego, California) GAIIx instrument using fourfold multiplexing.

### MethylC-seq library preparation {#s2-1-3}

For each accession two individual plants were pooled and total DNA was extracted using CTAB and phenol-chloroform. Approximately two micrograms of genomic DNA was used for MethylC-seq library construction and sequencing (92 bp paired-end) by BGI.

Sequence analysis {#s2-2}
-----------------

### Genome sequences {#s2-2-1}

Illumina sequencing data from 180 published Swedish genomes ([@bib29]) were combined with sequencing data from another 79 ([1001genomes.org](http://1001genomes.org)), which had been processed using the same pipeline to yield polymorphism data for a total of 259 accessions (including those used for MethylC-seq and RNA-seq here). Linkage disequilibrium calculated using the R package LDHeatmap (version 0.9.1; [@bib37]).

### RNA-seq data processing {#s2-2-2}

#### Read mapping {#s2-2-2-1}

After demultiplexing, 36 bp RNA-Seq reads were trimmed from barcodes (4 nt) and mapped to the TAIR10 reference genome including known variation with the PALMapper aligner ([@bib20]) using a variant-aware alignment strategy. Two different sources of variants were used: (1) single nucleotide variants (SNV) and structural variants (SV) from genome sequencing (2.1) and (2) SNVs and SVs called in an initial alignment round of the RNA-Seq reads to the TAIR10 reference genome with PALMapper (relevant parameters: -M 4 -G 4 -E 6 -I 25000 -NI 1 -S). For both sources of variants we applied stringent filter criteria to reduce false calls: (1) genome variants had to appear in at least 40 strains with a minor allele count of at least 5 strains, (2) RNA-Seq variants had to be confirmed by at least 2 alignments within the same strain and had to have less than factor 2 many non-confirming alignments within the same strain. Variants from both sources were integrated into one file that was used for a second, variant-aware alignment round with PALMapper (relevant parameters: -M 2 -G 0 -E 2 -I 5000 -NI 0 -S). In variant-aware alignment mode, PALMapper builds an implicit representation of the reference genome that reflects all possible variant combinations that exist for a genomic region. The output is automatically projected to the TAIR10 coordinate system. To account for reads ambiguously mapping to multiple locations in the genome, we used a custom python script ([Supplementary file 3](#SD3-data){ref-type="supplementary-material"}) to remove all reads that showed at least one mapping additional to the best hit with the same edit distance. Additional hits were only counted as ambiguous, if they differed at least 3 nt in start and stop coordinates to the best hit. On average, 5.7 M reads were mapped per sample after removal of ambiguous reads. Low complexity libraries with less than 30% of mappable reads or samples with less than 800,000 mappable reads (6 in total) where excluded from further analysis.

#### Sample validation {#s2-2-2-2}

To correct for possible sample or data mix-ups, SNP were called from the RNA-seq alignments using a custom python script and compared to independently collected SNPs from the *Arabidopsis* 250K SNP array (Supplementary materials; [@bib22]). Samples not matching the expected genotype were reassigned to the correct genotype where possible or otherwise excluded from further analysis.

#### Filter for gene expression quantification {#s2-2-2-3}

We quantified gene expression by counting the number of reads that were longer than 24 bp and that mapped to genes on all non-chloroplast and non-mitochondrial chromosomes. To obtain a stable quantification, we only used those reads that were uniquely mapped into the exonic regions of genes. Furthermore, we required that the reads did not map completely into regions where two genes overlap in order to avoid mixing quantifications of different genes. In the later text we will refer to this estimate as the raw expression estimate.

We also quantified the gene expression when additionally accounting for SV, alternative splicing and repetitive sequences that can all bias gene expression quantification. This estimate will be referred to as sv-corrected expression. For this quantification we additionally filtered for reads that start in an insertion or deletions and their two neighboring bases, that mapped into regions that are not contained in all transcripts of a gene and reads which were mapped completely into regions which are repetitive based on a 50 bp window.

#### Quantification per ecotype and environment {#s2-2-2-4}

After filtering (see 'Filter for gene expression quantification'), there were 499 RNA-Seq libraries left for analysis. Next, we merged libraries per ecotype and environment, yielding 323 unique merged RNA-Seq libraries for a unique ecotype and environment (160 in 10°C, 163 in 16°C).

#### Estimation of library size and abundance estimates {#s2-2-2-5}

We followed the low level normalization proposed by [@bib2], jointly applied to the set of expression estimates across ecotypes and environmental backgrounds. First, we estimated effective library sizes as the median expression estimates across all genes. Based on this, we derived correction factors to adjust individual libraries for differences in size.

#### RKPM values {#s2-2-2-6}

Library-size adjusted raw counts were used to obtain standard read counts per million expression estimates for each gene.

### MethylC-seq data processing {#s2-2-3}

#### Read mapping {#s2-2-3-1}

Reads were aligned as previously described ([@bib10]) to a modified pseudo-reference chromosome in which SNPs were inserted into the TAIR10 reference genome using NextGenMap (version 0.4.3; [@bib33]) allowing up to 10% mismatch between the reads (-i 0.90) and the reference sequence and discarding reads that map equally well to more than one genomic location or have less than 45 nucleotides mapping without error to the reference sequence (-R 45). Average coverage was 12.6 X.

To correct for sample or data mix-ups, the raw data was also aligned to the first chromosome of the Columbia-0 TAIR reference genome as described above and SNP calling performed using the BISsnp package ([@bib28]). The polymorphism data were then compared to data from genome sequencing ([1001genomes.org](http://1001genomes.org)). Accessions that did not have the highest similarity to the expected genotype were excluded from further analysis.

#### DNA methylation analysis {#s2-2-3-2}

Methylation was estimated individually for each cytosine using a python script provided with the BSMAP software package ([@bib44]). Conversion efficiency was estimated from the fraction of methylated cytosines in chloroplasts using the R software package ([www.r-project.org](http://www.r-project.org/), version 2.15.2). After eliminating one outlier, the samples had conversion efficiencies ranging from 99.25%--99.80% (mean = 99.59%). Genome wide average methylation levels were calculated separately for the CG, CHG and CHH contexts. The average variance between 11 biological replicates was 2.2%, 3.2% and 7.3% for CG, CHG and CHH methylation respectively, while for identical genotypes grown at different temperatures (111 pairs) CG, CHG and CHH methylation variance was 2.7%, 4.6% and 15.9% respectively. The variance in genome wide methylation levels for the 152 accessions grown at 10°C was respectively 7.6%, 9.2% and 13.2% for CG, CHG and CHH methylation, while for the 121 accessions grown at 16°C genome wide CG, CHG and CHH methylation varied 8.5%, 9.5% and 14.3% respectively.

The Bioconductor package Repitools (version 0.6.0; [@bib38]) was used to average methylation over genomic features of interest (e.g., all genes, all transposons over 4 kB or a subset of transposons of interest). Pairwise DMRs were called individually for each accession using the R software package methylKit (version 0.5.6; [@bib1]) using a window size of 100 bp, an FDR rate of 0.05 and a minimum fold change of 0.3. Overlap of DMRs with (TAIR10) genomic features such as transposons and genes was calculated using the Bioconductor package ChIPpeakAnno (version 2.8.0; [@bib47]). For each accession, methylation data was smoothed independently for each context using the Bioconductor package BSmooth (version 0.4.5; [@bib16]) using the default settings. Average methylation was then calculated for (overlapping) 200 bp sliding windows centered every 100 bp across the genome. Further analysis was limited to the 200,000 windows showing the most variance among accessions.

Population genetic analysis {#s2-3}
---------------------------

### GWAS {#s2-3-1}

Linear mixed models that correct for confounding by the genetic background using a kinship matrix calculated from genetic data were used throughout ([@bib21]; [@bib34]). To examine the effect of genotype on local CHH methylation variation, DMRs were defined by filtering the 200 bp methylation windows to remove those containing missing data (no coverage) in one or more accessions, then selecting the 10^5^ remaining windows with the greatest variance in DNA methylation. For GBM, genes were filtered to remove those that had more than 0.05 average CHG methylation or less than 0.05 average CG methylation across the accessions ([Figure 6---figure supplement 2](#fig6s2){ref-type="fig"}).

### Variance component analysis {#s2-3-2}

To investigate the relative contributions of genetic and environmental effects to methylation differences we used LIMIX ([@bib27]), which efficiently estimates variance components using linear mixed models.

For each DMR, we considered a linear mixed model with a fixed effect for the environment and random effects for the contributions from *cis* and *trans* genetic variants and variants from the CMT2 locus. Indicating with N and E respectively the number of samples and environments (E = 2), the NxE multivariate phenotype **Y** can be written as$$\mathbf{Y} = \mathbf{1}_{N,1}\mu^{T} + \mathbf{U}^{\text{CMT}2} + \mathbf{U}^{\text{cis}} + \mathbf{U}^{\text{trans}} + \mathbf{\psi},$$where **μ** is a Ex1 vector of environment-specific mean values, and$$\mathbf{U}^{\text{CMT}2} \sim \text{MVN}\left( {\mathbf{0},\mathbf{C}^{\text{CMT}2},\mathbf{R}^{\text{CMT}2}} \right),\,\mathbf{U}^{\text{cis}} \sim \text{MVN}\left( {\mathbf{0},\mathbf{C}^{\text{cis}},\mathbf{R}^{\text{cis}}} \right),$$$$\mathbf{U}^{\text{trans}} \sim \text{MVN}\left( {\mathbf{0},\mathbf{C}^{\text{trans}},\mathbf{R}^{\text{trans}}} \right),\,\mathbf{\psi} \sim \text{MVN}\left( {\mathbf{0},\mathbf{\Sigma},\mathbf{I}_{N}} \right),$$where MVN(**0**,**C**,**R**) denotes a matrix normal distribution with mean **0**, column covariance matrix **C** and row covariance matrix **R**. **R**^cis^ and **R**^trans^ indicate the genetic relatedness matrices based on cis and trans variants respectively, where all variants within 50 kb from the DMR were defined as *cis*-acting and all others as *trans*--acting. Similarly, **R**^CMT2^ denotes the genetic relatedness matrix based on genotypes at the CMT2 locus. The row covariance of the noise component ***I~N~*** corresponds to an *N* x *N* identity matrix.

The covariance matrices **C**^CMT2^, **C**^cis^, **C**^trans^ and **Σ** describe phenotypic correlations across environments due to these contributions, and were estimated from the data using maximum likelihood. For each DMR, we considered up to 10 random restarts for the optimization and stopped as soon as convergence was achieved. DMRs for which no convergence was achieved were discarded from genome-wide summary statistics.

Once the model parameters have been estimated, the variance explained by environment can be calculated from **μ**, while environment-persistent and environment-specific effects from a given random effect can be estimated by decomposing the corresponding trait covariance into a shared and an independent component ([@bib27]).

### QTL mapping {#s2-3-3}

MethylC-seq data for the 113 F2 individuals was mapped as described in section 'Read mapping' to the Columbia-0 TAIR reference genome. SNP-calling was done directly from the methylC-seq data using the BIS-SNP package ([@bib28]). From these SNPs local haplotype was inferred for sequential 500 Mb windows which were then used to create a haplotype map using the R package R/qtl ([@bib6]). Mapping was done using Haley-Knot regression ([@bib3]) with a 4 centimorgan steps size. Genome wide significance was estimated by permutation testing (1000 permutations).

### DMR calling on DNA methylation mutants {#s2-3-4}

Pairwise DMRs were called for T-DNA mutants vs the wild-type control using the R software package methylKit (version 0.5.6; [@bib1]) using a window size of 100 bp, an FDR rate of 0.05 and a minimum fold change of 0.3. Overlap between these DMRs and 'temperature DMRs" calculated for the accessions was calculated and significance testing (Fisher\'s exact test) was calculated using R software.

### *Qst-Fst* test {#s2-3-5}

*Fst* was computed using the Hudson estimate as suggested in [@bib5]. We note that our estimate of 0.187 is consistent with the recent estimate of [@bib18] (although the samples only overlap in part). For *Qst*, we first estimated northern, southern, and overall additive variance using the Hasemann-Elston regression, and a SNP-based identity-by-state matrix ([@bib8]), then calculated *Qst* as $\sigma_{B}^{2}/\left( \sigma_{B}^{2} + 2\sigma_{w}^{2} \right)$, where $\sigma_{w}^{2}$ is the weighted average of variance within north and south populations, and $\sigma_{B}^{2}$ is the variance between populations, obtained by subtracting $\sigma_{w}^{2}$ from the overall additive variance.

Funding Information
===================

This paper was supported by the following grants:

-   http://dx.doi.org/10.13039/100000051National Human Genome Research Institute (NHGRI) P50HG002790 to Richard M Clark, Magnus Nordborg.

-   http://dx.doi.org/10.13039/501100000781European Research Council (ERC) 268962 to Magnus Nordborg.

-   http://dx.doi.org/10.13039/501100000780European Commission Marie Curie FP7 fellowship 253524 to Oliver Stegle.

-   http://dx.doi.org/10.13039/501100000780European Commission European Community Framework Programme 7 283496 to Magnus Nordborg.

-   http://dx.doi.org/10.13039/100000002National Institutes of Health (NIH) GM07464 to Edward J Osborne.

-   http://dx.doi.org/10.13039/501100002428Austrian Science Fund (FWF) M1369 to Manu J Dubin.

This work was supported multiple sources: the National Human Genome Research Institute of the US National Institutes of Health (P50HG002790 to MN and RMC; PIS Tavaré); the European Research Council (268962 MAXMAP to MN, Marie Curie FP7 fellowship 253524 to OS); the European Community Framework Programme 7 (283496 transPLANT to MN); the National Institutes of Health Genetics (Training Grant GM07464) to EJO; the Austrian Science Fund (FWF M1369) to MJD; as well as core funding at the GMI. Sequencing was carried out by the Epigenome Center at: the University of Southern California; BGI; and the Campus Support Facilities at the Vienna Biocenter. The authors wish to thank: Andrew Smith for preliminary analyses and effectively supervising the project during the transition of MN from USC to GMI; J Bergelson for providing seed; many members of MN\'s lab for discussion (in particular D Filiault, F Rabanal, P Novikova, E Kerdaffrec and Ü Seren for help with various analysis tasks); A Hancock and C Huber for discussions about selection; F Sedlazeck, and P Rescheneder for advice on alignment; and F Berger and O Mittelsten Scheid for discussions and comments on the paper.

Additional information {#s3}
======================

The authors declare that no competing interests exist.

MJD, Analysis and interpretation of data, Drafting or revising the article.

PZ, Acquisition of data, Analysis and interpretation of data.

DM, Analysis and interpretation of data.

GJ, Analysis and interpretation of data.

BV, Analysis and interpretation of data.

QS, Analysis and interpretation of data.

M-SR, Acquisition of data.

JJ, Acquisition of data.

SI, Acquisition of data.

VV, Acquisition of data.

EJO, Analysis and interpretation of data, Contributed unpublished essential data or reagents.

FPC, Analysis and interpretation of data, Contributed unpublished essential data or reagents.

PD, Analysis and interpretation of data, Contributed unpublished essential data or reagents.

AK, Analysis and interpretation of data, Contributed unpublished essential data or reagents.

GR, Analysis and interpretation of data, Contributed unpublished essential data or reagents.

OS, Analysis and interpretation of data, Contributed unpublished essential data or reagents.

QL, Provided pre-publication access to genome sequencing data.

RMC, Conception and design, Contributed unpublished essential data or reagents.

MN, Conception and design, Analysis and interpretation of data, Drafting or revising the article.

Additional files {#s4}
================

10.7554/eLife.05255.027

###### 

Multiple sequence alignment of CMT2 sequences from different accessions.

**DOI:** [http://dx.doi.org/10.7554/eLife.05255.027](10.7554/eLife.05255.027)

10.7554/eLife.05255.028

###### 

Multiple sequence alignment of CMT2 sequences from different accessions in FASTA format.

**DOI:** [http://dx.doi.org/10.7554/eLife.05255.028](10.7554/eLife.05255.028)

10.7554/eLife.05255.029

###### 

Source for scripts used to extract SNPs from RNA-seq data and comparing to existing 250k SNP data.

**DOI:** [http://dx.doi.org/10.7554/eLife.05255.029](10.7554/eLife.05255.029)

Major dataset {#s4-1}
-------------

The following datasets were generated:

Dubin MJ, Nordborg M, 2014,DNA methylation variation in Arabidopsis has a genetic basis and appears to be involved in local adaptation,<http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE54292>,Publicly available at NCBI Gene Expression Omnibus GSE54292.

Osborne EJ, Zhang P, Remigereau MS, Drewe P, Rätsch G, Nordborg M, Clark RM, 2014,DNA methylation variation in Arabidopsis has a genetic basis and appears to be involved in local adaptation,<http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE54680>,Publicly available at NCBI Gene Expression Omnibus GSE54680.

Dubin MJ, Nordborg M, 2014,Population epigenomics of Swedish Arabidopsis,<http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE65685>,Publicly available at NCBI Gene Expression Omnibus GSE65685.

Dubin MJ, Nordborg M, 2014,Effect of Chromomethylase-2 and Dicer-like-3 knockouts on DNA methylation,<http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE66017>,Publicly available at NCBI Gene Expression Omnibus GSE66017.

Long Q, Rabanal FA, Meng D, Huber CD, Farlow A, Platzer A, Zhang Q, Vilhjálmsson BJ, Korte A, Nizhynska V, Voronin V, Korte P, Sedman L, Mandáková T, Lysak MA, Seren Ü, Hellmann I, Nordborg M, 2013,Massive genomic variation and strong selection in Arabidopsis thaliana lines from Sweden,<http://plone.gmi.oeaw.ac.at/downloads/nordborg/data-release-for-massive-genomic-variation-and-strong-selection-in-arabidopsis-thaliana-lines-from-sweden>,Publicly available at Gregor Mendel Institute of Molecular Plant Biology.

10.7554/eLife.05255.030

Decision letter

Dermitzakis

Emmanouil T

Reviewing editor

University of Geneva Medical School

,

Switzerland

eLife posts the editorial decision letter and author response on a selection of the published articles (subject to the approval of the authors). An edited version of the letter sent to the authors after peer review is shown, indicating the substantive concerns or comments; minor concerns are not usually shown. Reviewers have the opportunity to discuss the decision before the letter is sent (see [review process](http://elifesciences.org/review-process)). Similarly, the author response typically shows only responses to the major concerns raised by the reviewers.

Thank you for sending your work entitled "DNA methylation variation in Arabidopsis has a genetic basis and shows evidence of local adaptation" for consideration at *eLife*. Your article has been favorably evaluated by Stylianos Antonarakis (Senior editor), a Reviewing editor, and three reviewers..

The Reviewing editor and the other reviewers discussed their comments before we reached this decision, and the Reviewing editor has assembled the following comments to help you prepare a revised submission.

The work was considered important and novel and no major technical concerns were raised.

Below are the key issues to be addressed in revision:

1\) What is the significance of the identified variation in CMT2-dependent CHH methylation (in contrast to others species like maize that don\'t have this enzyme) and the same for the GBM variation? Do these results reflect differences in chromatin compaction? What does an increase in ∼10% gene body methylation actually mean? Is there cell type specificity? Is the CG methylation higher at each CG site or are more CG sites methylated or is it a combination of both? What is the explanation for differential CHH between 10 and 16C?

2\) The present data does not show that the levels of DNA methylation influence local adaptation. Rather, it shows that variation in methylation has a strong genetic basis. Furthermore, there is no evidence showing that this genetic variation underlies local adaptation. The statement about adaptation should be toned down or removed.

10.7554/eLife.05255.031

Author response

*Below are the key issues to be addressed in revision*:

1\) What is the significance of the identified variation in CMT2-dependent CHH methylation (in contrast to others species like maize that don\'t have this enzyme)...

CMT2 appears to be present in most plants, with the notable exception of maize (at least for the reference strain), suggesting it plays an important role. It has been suggested (Shen, 2014) that differences in CMT2 dependent CHH methylation affect tolerance to thermal stress (which would certainly fit with maize being very tolerant to heat stress). It could also potentially affect TE expression/expansion, although we have been unable to find evidence for this in our data. We now mention these possibilities in the text.

... and the same for the GBM variation?

The function of GBM is generally poorly understood, but some involvement in transcriptional regulation seems highly likely (GutierrezArcelus, 2013). In other data, the variation in GBM between north and south is associated with differences in the level as well as variance of expression and is also strongly correlated with climate. It maps to loci that have previously been shown to be involved in local adaptation to climate. To strengthen this further, we have added an FstQst analysis that demonstrates that the trait is under divergent selection. A deeper understanding of the mechanisms involved will have to await the cloning and characterization of the transacting loci affecting GBM.

Do these results reflect differences in chromatin compaction?

Given that most of the variation in DNA methylation was fairly mild (max ∼10%) we consider it unlikely that they cause largescale changes in chromatin organization, especially when one considers that very few null mutants involved in DNA methylation (with a few notable exceptions such as mom1, ddm1 and met1) seem to affect chromatin compaction. That said, natural variation in chromatin compaction clearly does exist, at least in the Cvi accession (Tessodori, 2009) and we have mentioned the possibility. Future studies should investigate this.

What does an increase in ∼10% gene body methylation actually mean? Is there cell type specificity? Is the CG methylation higher at each CG site or are more CG sites methylated or is it a combination of both?

While DNA methylation variation between cell types has been described in the endosperm and pollen (for example Calarco, 2012) and also the germline (Baubec, 2014), in all of these cases it seems to be limited to transposons and mostly affect CG methylation. Furthermore, several investigations have shown that CG gene body methylation is very highly conserved between different organs (root, shoot, rosette, cauline leaf etc; ColemanDerr & Zilberman, 2012; Seymore DK, 2014). Given that our samples were all rosettes sampled at the same developmental stage, cell type should have very little if any affect on GBM. We have expanded our analysis to demonstrate that roughly the same number of CG sites are methylated in north and south, with the difference being a skew towards higher values in the north ([Figure 6--figure supplement 3](#fig6s3){ref-type="fig"} and [Figure 6--figure supplement 4](#fig6s4){ref-type="fig"}).

What is the explanation for differential CHH between 10 and 16C?

The difference appears to affect primarily RdDM target sites. Given that the temperature sensitivity of the smallRNA processing pathway is well documented (see, e.g., Romon, 2013), this would seem to be the most plausible explanation. This is noted in the text.

*2) The present data does not show that the levels of DNA methylation influence local adaptation. Rather, it shows that variation in methylation has a strong genetic basis. Furthermore, there is no evidence showing that this genetic variation underlies local adaptation. The statement about adaptation should be toned down or removed*.

Here we must respectfully disagree. It is certainly correct that we have no evidence that DNA methylation is somehow important for adaptation independent of genetic variation, and indeed we show that it has a strong genetic basis. However, it is wrong to say that there is no evidence for selection. The variants affecting GBM are strongly overrepresented in sweep regions, and show precisely the kinds of climate associations that would be expected if they were involved in local adaption. Furthermore, the northsouth divergence in GBM is far too great to be due to genetic drift (we have added an FstQst analysis to make this clearer). Of course the evidence is all indirect, and of course all these analyses come with significant caveats, but this is usually the case in population genetics. Generating direct evidence of fitness differences is extremely difficult and time consuming.

We are well aware that a number of strange claims concerning the role of DNA methylation are currently being made, and it precisely for this reason that we try to make it very clear that we are talking about genetic effects in this paper (DNA methylation is essentially treated as a phenotype).

[^1]: These authors contributed equally to this work.

[^2]: r = Pearson\'s correlation, rho = Spearman\'s rank correlation, p-value = significance of rho.

[^3]: PAR = photosynthetically active radiation.
