Interspecific variation across angiosperms in global DNA methylation: phylogeny, ecology and plant features in tropical and Mediterranean communities
Summary
- The interspecific range of epigenetic variation and the degree to which differences between angiosperm species are related to geography, evolutionary history, ecological settings or species-specific traits, remain essentially unexplored. Genome-wide global DNA cytosine methylation is a tractable ‘epiphenotypic’ feature suitable for exploring these relationships.
- Global cytosine methylation was estimated in 279 species from two distant, ecologically disparate geographical regions: Mediterranean Spain and tropical México. At each region, four distinct plant communities were analyzed.
- Global methylation spanned a 10-fold range among species (4.8–42.2%). Interspecific differences were related to evolutionary trajectories, as denoted by a strong phylogenetic signal. Genomes of tropical species were on average less methylated than those of Mediterranean ones. Woody plants have genomes with lower methylation than perennial herbs, and genomes of widespread species were less methylated than those of species with restricted geographical distribution.
- The eight communities studied exhibited broad and overlapping interspecific variances in global cytosine methylation and only two of them differed in average methylation. Altogether, our broad taxonomic survey supported global methylation as a plant ‘epiphenotypic’ trait largely associated with species evolutionary history, genome size, range size and woodiness. Additional studies are required for better understanding the environmental components underlying local and geographical variation.
Introduction
Angiosperms are the largest group of extant plants, dominating vegetation in most terrestrial environments and displaying an immense variety of phenotypic features. Despite the multivariate nature of plant phenotypes, recent across-species analyses suggest that extant plants have converged towards a relatively small set of successful trait combinations that can be defined straightforwardly (Díaz et al., 2016). Constraints operating in the multivariate phenotypic space illustrate the functional basis underlying some long-standing comparisons, such as those between woody and nonwoody plants, and between species that differ in leaf size, construction costs and growth rate (Díaz et al., 2016). Such functional phenotypes can be further related to both species-specific, intrinsic traits (e.g. genomic) and extrinsic factors (e.g. environmental), and the simultaneous consideration of these two classes of correlates offers a fruitful avenue for understanding ecological and evolutionary processes occurring at organismal, community and ecosystem scales (Soudzilovskaia et al., 2013; Kunstler et al., 2016; Funk et al., 2017).
Angiosperm nuclear genomes vary widely in size, chromosome number and ploidy level (Fedoroff, 2012; Wendel, 2015; Springer et al., 2016). Genome size is relatively easy to measure (Dolezel et al., 2007) and its broad variation has been investigated by plant biologists for decades (Greilhuber & Leitch, 2013; Leitch & Leitch, 2013). Current data indicate that genome size is related both to plant functional traits (e.g. growth rate, stomata size, phenology, perenniality) and environmental features (Vidic et al., 2009; Hodgson et al., 2010; Leitch & Leitch, 2013). Correlations also have been found with altitude, temperature and precipitation, although these ecological correlates should be confirmed by analyses controlling for phylogenetic relationships (Greilhuber & Leitch, 2013). Differences among species in genome size and complexity seem to be the outcome of cycles of upsurge and silencing of transposable elements, and whole-genome duplication events that have contributed to plant diversification and for which epigenetic regulation is crucial (Fedoroff, 2012; Wendel, 2015; Springer et al., 2016). Thus, epigenetic features involving histone modifications, noncoding RNAs and the methylation of DNA, although conceptually and functionally related to the well-studied topic of genome size, are much less investigated (Fedoroff, 2012).
Epigenetic features differ broadly within the angiosperms (Springer et al., 2016). However, there remains a need for a consensus on an easily measured variable that allows exploring the extent to which epigenetic differences between species are related to biogeography, evolutionary history, ecological settings and functionally relevant specific traits. DNA methylation is a key part of plant epigenomes (Fedoroff, 2012). Recent research has pointed out that global DNA cytosine methylation (‘global methylation’ hereafter) varies considerably among species and is evolutionarily related to variation in genome size (Alonso et al., 2015; Niederhuth et al., 2016; Vidalis et al., 2016). The global methylation level is a genomic feature that can be interpreted as an ‘epiphenotypic’ trait useful for undertaking comparative studies aimed at exploring the factors triggering its broad interspecific variation. It can be analyzed using a method based on high-performance liquid chromatography (HPLC) that does not require a full genome sequence (Alonso et al., 2016b) and provides cost-effective estimates of the proportion of genomic cytosines that are methylated, analogous to flow cytometry for quickly measuring genome size (Dolezel et al., 2007). Within species, global methylation is variable among populations, individuals within populations and different modules of a single plant, and such variation is related to key functional traits such as seed size and seed production (Alonso et al., 2014, 2018; Herrera et al., 2019), suggesting that its broad interspecific variance could have some ecological relevance. Although data on global methylation (epiphenotype) will be inherently unable to convey information on specific cytosine methylation patterns (epigenotype), it will still provide complementary, interpretable data bearing on the extent of global cytosine methylation that may be functionally linked to other epigenetic features (Vidalis et al., 2016). Comparisons of global methylation levels also can be useful for understanding the epigenomic changes that occur at both intra- and interspecific scales (Alonso et al., 2015, 2016a; Alonso et al., 2016b; Vidalis et al., 2016).
Understanding the ecological causes and consequences of epigenetic variation in natural populations has been singled out as a fundamental ecological question (Sutherland et al., 2013; Richards et al., 2017). In plants, analyses of extent and patterns of epigenetic variation at the intraspecific level have documented broad variation at different geographical scales that correlates to some extent with environmental gradients (Medrano et al., 2014; Schulz et al., 2014; Kawakatsu et al., 2016). To the best of the present authors’ knowledge, there is no previous study analyzing the ecological factors correlated with interspecific differences in epigenetic features. In this study, a plant community approach was adopted for broadening the taxonomic and environmental conditions that the study species encompass. Altogether, global methylation was analyzed in 279 angiosperm species from distinct communities distributed along an environmental gradient in two distant geographical regions: Mediterranean southern Spain and tropical southeastern México. The primary aim of this study was to assess the degree to which interspecific differences in global cytosine methylation were related to broad-scale factors such as phylogeny and geography, more proximate environmental conditions (i.e. local community), species-specific functional traits (e.g. growth form, leaf phenology) and size of geographical range. In particular, applying such phylogenetically informed methods was expected to establish: a strong divergence in global methylation levels among species of the Mediterranean and tropical study regions, according to the contrasting glacial histories and environmental features that plants have experienced in them; divergence in global methylation across communities in the two regions, according to their distinct environments and provided that coexisting species within a certain community would be more similar than across communities; and that genomes of species with divergent functional phenotypic traits also would differ in global methylation levels. Results showed that taxonomic affiliation, geography, distribution range and functional plant traits were much more important than small-scale environmental variation across plant communities as predictors of the global methylation level of individual species.
Materials and Methods
Study communities and species sampling
Two geographically distant, ecologically disparate regions with contrasting glacial histories were chosen for this study: Mediterranean Andalucía (southern Spain) and the tropical Yucatán Peninsula (southeastern México). Sampling species in these two regions aimed to broaden the taxonomic and environmental ranges of the analysis. In each region, four different vegetation types were selected along a rainfall gradient (Fig. 1; Table 1). One vegetation type, coastal sand dune, could be found in the two regions and, thus, was sampled in both Yucatán and Andalucía. The other three vegetation types were broadly defined by climate and orographic gradients specific to each region (Table 1). In Yucatán, the gradient was largely associated with rainfall seasonality and distance to the Atlantic coast (Flores & Espejel, 1994; García, 2004), whereas in Andalucía it encompassed a wider longitudinal and elevation range, including locations from the driest Atlantic coast to the wettest mountainous areas in Sierra de Cazorla (Ramos-Calzado et al., 2008; Gómez Mercado, 2011).
| Site | Code | Coordinates | Dist (km) | Altitude (m asl) | Precip (mm) | Climate | n | Vegetation type | Other features | References |
|---|---|---|---|---|---|---|---|---|---|---|
| Spain | ||||||||||
| Natural Area Arenales de Punta Umbría | M1 | 37°11′22″N/6°59′37″W | 0 | < 10 | 490 | Thermo-mediterranean | 29 | Coastal vegetation | Sandy soil, halophytic, windy | García-Mora et al. (1999) |
| Pinares Aznalcázar | M2 | 37°13′33″N/6°12′7″W | 40 | < 50 | 630 | Thermo-mediterranean | 45 | Pine woodland with cork oaks and xerophytic scrubs | Sandy soil, forest fragment | Aparicio (2008) |
| Sierra de Cazorla – Guadahornillos/Borosa valley | M3 | 38°0′46″N/2°51′35″W | 335 | 760–1100 | 985 | Meso-mediterranean | 57 | Arborescent matorral and scrubland | Mountainous pristine area, limestone | Gómez Mercado (2011) |
| Sierra de Cazorla -Nava Correhuelas/Cabeza del Tejo | M4 | 37°56′N/2°52′W | 340 | > 1500 | 1128 | Supra-mediterranean | 35 | Dwarf cushion-heaths and grasslands | Mountainous pristine area, snowfalls | Gómez Mercado (2011) |
| Mexico | ||||||||||
| Dunes Reserva Estatal El Palmar | T1 | 21°11′00″N/90°00′00″W | 0 | < 10 | 600 | Tropical warm, semi-arid | 35 | Coastal vegetation | Sandy soil, halophytic, windy | Flores & Espejel (1994) |
| Sierra Papacal | T2 | 21°08′13″N/89°46′54′′W | 20 | 8 | 650–800 | Tropical warm, subhumid | 20 | Lowland dry thorny forest | Sandy and clays soil | Flores & Espejel (1994) |
| Xmatkuil | T3 | 20°52′07″N/89°36′52″W | 45 | 10 | 728–1005 | Tropical warm, subhumid | 20 | Lowland dry forest | Predominant clay soil | Mizrahi et al. (1997) |
| Reserva Biocultural Kaxil Kiuic | T4 | 20°06′33″N/89°32′55″W | 130 | 60–160 | 1200 | Tropical warm, subhumid | 38 | Mid-elevation semi-deciduous forest with 13–18 m trees | Limestone hills alternate with flat areas with diverse soil types | Flores & Espejel (1994) |
At each site, 20–57 native plant species were chosen among the locally most abundant ones, encompassing a broad variety of life-forms and phylogenetic diversity as inferred by the number of genera represented. Species occurring in several communities were considered only in the one where it was more common. We collected 3–10 full-grown, current season leaves per species, each from a distinct, healthy adult individual. Samples were collected separately in numbered individual paper envelopes, dried at ambient temperature in silica gel and preserved in this way by periodic replacement of the desiccant until processing.
Species information search
Taxonomy followed Castroviejo (1986), Valdés et al. (1987), and/or Blanca et al. (2011) for the Spanish species and Standley (1930) for the Mexican ones. Order- and family-level classification followed the IV Angiosperm Phylogeny Group classification (Byng et al., 2016) according to the APG website (http://www.mobot.org/MOBOT/research/APWeb/, Stevens, 2001 onwards). Online repositories and plant databases Tropicos® (http://www.tropicos.org), the International Plant Name Index (http://www.ipni.org), the Encyclopedia of Life (http://www.eol.org), the Anthos project of Spanish Plant Information System (www.anthos.es, Anthos, 2016), and the eFlora project (http://www.efloras.org) were the main information sources used for verifying names and compiling data on growth form, leaf phenology and geographical distribution for each species (last accessed in 15 November 2016).
All species were categorized with regard to their growth form (woody vs nonwoody) with woody species defined as those maintaining an aboveground stem with the ability to produce hard lignified tissues. Woody species were further categorized by leaf phenology (deciduous vs evergreen). The area of the geographical distribution for each species was used as a proxy of its range size. Range size is one fundamental ecological and evolutionary characteristic of a species (Thompson et al., 1999; Gaston, 2003; Gaston & Fuller, 2009) that could be considered, at least partially, as a spatial representation of its degree of specialization (Devictor et al., 2010) or environmental tolerance (Sheth & Angert, 2014). As a fine quantification of range size was beyond the scope of the present study, information obtained about the area of geographical distribution for each species was converted to a categorical scale with three levels: (1) local, restricted to a small and well-defined region (e.g. íbero-magrebí, Yucatán and Belice), which may be indicative of species with particular habitat requirements or narrow environmental tolerance; (2) common, present across several European countries, or across several areas in North and Central America, species that may be able to grow in a greater set of environments; and (3) widespread, distributed over a wide geographical area (e.g. the whole Mediterranean Basin, Eurasia or America), species that may have the highest environmental tolerance.
Information on the genome size of the species included in the present sample was obtained from the Kew Royal Botanic Gardens Angiosperm C-value Database (Bennett & Leitch, 2012), which is the standard source for contemporary work on genome-size evolution in plants (Bennett & Leitch, 2005, 2012; Leitch & Leitch, 2013). C-value data, corresponding to the amount in picograms of DNA contained within a haploid nucleus, were available for 84 species in our sample (30.1% of total; Supporting Information Table S1).
Lab methods
Two individuals were analyzed per species. Dry leaves were homogenized to a fine powder using a Retsch MM 200 mill (Retsch, Haan, Germany). Total genomic DNA was extracted using Qiagen DNeasy Plant Mini Kit. A 100-ng aliquot was digested with three units of DNA Degradase PlusTM (Zymo Research, Irvine, CA, USA), a nuclease mix that degrades DNA to its individual nucleoside components. Digestion was carried out in a 40-μl volume at 37°C for 3 h, and terminated by heat inactivation at 70°C for 20 min. Digested samples were stored at −20°C until analyzed. Two independent replicates of digested DNA per sample were processed to estimate global cytosine methylation. Selective chemical derivatization of cytosine moieties with 2-bromoacetophenone under anhydrous conditions was carried out to induce fluorescence and subsequent reverse phase high-performance liquid chromatography (HPLC) with spectrofluorimetric detection allowed quantification of each separated product. The percentage of total cytosine methylation on each replicated sample was estimated as 100 × 5 mdC/(5 mdC + dC), where 5 mdC and dC are the integrated areas under the peaks for the methylated (5-methyl-2′-deoxycytidine) and nonmethylated (2′-deoxycytidine) forms, respectively (see Alonso et al., 2016b, for details).
Statistical analyses
Multi-species datasets challenge the key assumption of independent (uncorrelated) errors required for classical linear regression and ANOVA tests, because closely related species will tend to be more similar than distantly related species in many ways and functional features. Several methods have been developed to control for this (Lajeunesse & Fox, 2015; Uyeda et al., 2018). In the present study, three complementary analytical approaches were adopted to examine variation in global methylation – the focal variable of interest – namely: (1) analyzing its relationship with phylogeny by evaluating the phylogenetic signal in our multispecies sample; (2) applying Phylogenetically Independent Contrasts (PICs; Felsenstein, 1985; Pagel, 1999) for analyzing the covariance of the response with another continuous variable (i.e. genome size) along the inferred phylogeny; and (3) fitting linear mixed models to test for the significance of relationships between the response variable and other categorical traits modeling variance according to taxonomic affiliation (Maddison & FitzJohn, 2015). Each of these approaches is described in detail below. Before conducting the analyses, global methylation data were visually inspected to corroborate the absence of obvious outliers and normality of distributions. Unless otherwise stated analyses were conducted using the R environment (R Core Team, 2017).
Phylogenetic signal of global methylation
Evolutionary relationships among the 279 Angiosperm species included in the sample were inferred by constructing a phylogenetic tree using the phylomatic (Webb & Donoghue, 2005) tool bundled in Phylocom 4.2 software (Webb et al., 2008). All tree branch lengths were set to unity and polytomies (i.e. relationships that were not resolved) were solved randomly.
The relationship between phylogeny and interspecific variation in global methylation in the present sample was examined by testing for the presence of a phylogenetic signal in the data, defined as ‘a tendency for related species to resemble each other more than they resemble species drawn at random from the tree’ (Blomberg & Garland, 2002). The strength, or ‘effect size’, of phylogenetic structuring was assessed using Blomberg's K and Pagel's λ, two indices that are robust to missing branch length and assume a Brownian Motion (BM) model of trait evolution (see Münkemüller et al., 2012 for details). For both indices, values close to zero denote phylogenetic independence and a value of unity indicates a trait distribution as expected under BM (i.e. the trait evolves independently in each lineage and its variance between two species is proportional to their divergence time). Computations were performed with functions in the R/phytools (Revell, 2012). Statistical significance was tested by randomization with 9999 repetitions. Both λ and K indices assume constant rates of evolution across time and clades, which could be fairly unrealistic on the phylogenetic scale studied here. To test whether phylogenetic signal varies between clades and detect local hotspots of autocorrelation an adaptation of the local Moran's I (Local Indicators of Spatial Association; Anselin, 1995; Keck et al., 2016) was applied for phylogenetic data (Local Indicators of Phylogenetic Association (LIPA); Keck et al., 2016). Moran's I varies from −1 to +1, with positive values indicating that the trait of interest is more similar than expected by chance and negative values indicating an absence of correlation between the trait and interspecies phylogenetic distance. Local Moran's I (Ij) values were calculated based on patristic distances using the ‘lipaMoran’ function of R/phylosignal and mapped onto the phylogeny. Significance in Moran's I was assessed through 9999 permutations. R/ggtree was used for visualizing and annotating the phylogenetic tree with associated data (Yu et al., 2017).
Correlated evolution between global methylation and genome size
Correlated evolution between global methylation and genome size was tested in the subset of 84 species with complete data. Genome size data were log10-transformed for analyses. The PICs method was used, and contrasts were obtained with the ‘pic’ function in R/ape (Popescu et al., 2012). Departures from the BM model underlying the PICs method (continuous traits evolve randomly in any direction and amounts of change are normally distributed) may contribute to inflating Type I error (Uyeda et al., 2018). Normality of contrasts, a condition expected under BM (Paradis, 2012), was tested with the Shapiro–Wilk W test. The significance of the relationship between the contrasts for global methylation and the contrasts for genome size (log10-transformed) was tested by fitting a linear regression through the origin and testing its significance with a permutation procedure with 9999 repetitions, which does not require to make any assumptions about the sampling distribution (‘lmorigin’ function in R/ape; Paradis, 2012; Popescu et al., 2012).
Interspecific variation in global methylation: a linear mixed-model approach
In order to obtain a hierarchical partition of variance of the response variable, a fully random, nested linear model was fitted to global methylation data in which taxonomic affiliation levels (class/order/family/genus/species) and individuals (two replicates, nested within species) were treated as random effects.
Further, linear mixed models were used to analyze changes in global methylation levels as a function of the fixed predictor effects of interest (see the ‘Geographical and ecological correlates’ section below) while including order (n = 37) as the single random effect to account for taxonomic nonindependence without exhausting the degrees of freedom available in our sample (Bolker, 2015). The ‘lme’ function from R/nlme (Pinheiro et al., 2017) was used to fit mixed-effect models with Restricted Maximum Likelihood (REML) estimation. Four different models were fit to test for the effect of each of the following fixed categorical effects (1) geographical region (Mediterranean vs tropical), (2) growth form (woody vs nonwoody), (3) leaf phenology of woody plants (deciduous vs evergreen; n = 102) and (4) size of geographical range (local, common, widespread). Including order as random effect improved goodness-of-fit compared to the corresponding linear model without random effects and was included in all analyses of the full dataset, residuals were inspected with quantile–quantile plots and normal distribution of residuals was not rejected (P < 0.01) in any of the fitted models presented.
Results
Global methylation varied nearly 10-fold across species, ranging between 4.8% (Bauhinia divaricata, Fabaceae) and 42.2% (Pancratium maritimum, Amaryllidaceae) (interquartile range = 15.7–25.6%; mean ± SE = 20.6 ± 0.4, n = 279; Table S1). The phylogenetic, taxonomic, geographical and ecological correlates of this interspecific variation are dissected in the following sections.
Phylogenetic and taxonomic correlates
Indices used to test for a phylogenetic signal in the global methylation dataset yielded statistically significant results (Table 2), thus strongly supporting the view that interspecific variation in global methylation level was phylogenetically structured; that is, species that share a recent common ancestor tend, on average, to have genomes with global methylation level more similar to one another than to more distant relatives (Fig. 2). Pagel's λ estimate for the entire dataset (0.835) was close to unity, the expected value for trait evolution under a pure BM model. These findings were partially corroborated when species from each region were analyzed separately (Table 2). In particular, the phylogenetic signal was stronger in the Mediterranean species sample, which included 33 orders, 130 genera and 166 species, than in the taxonomically less diverse tropical sample, with only 22 orders, 103 genera and 113 species, and the full sample. LIPA analyses revealed significant local positive autocorrelation for global methylation (Fig. 2 right panel, red dots) in some basal lineages (Magnoliids and Monocots) and at several nodes within Eudicots, particularly within the Rosiids clade. Significant local associations within Rosiids were found in species groups with global methylation levels above (Fig. 2 middle panel, positive values) and below (Fig. 2 middle panel, negative values) the average methylation level.
| Dataset | Bloomberg's K | P-value | Pagel's λ | P-value |
|---|---|---|---|---|
| All species | 0.385 | 0.0001 | 0.835 | 0.0001 |
| Mediterranean | 0.779 | 0.0001 | 0.935 | 0.0001 |
| Tropical | 0.275 | 0.23 | 0.549 | 0.048 |
- Statistical significance was evaluated with 9999 randomization tests. Values of Bloomberg's K and Pagel's λ close to 0 indicate no relationship between variance in the trait and phylogenetic distance among species, whereas values close to 1 indicate that resemblance between relatives approaches that expected under Brownian motion evolution.
The PICs for global methylation and genome size were positively correlated in the subset of 84 species with genome size data available, indicating correlated evolution between both plant features. The regression through the origin of PIC global methylation on PIC log(Cvalue) was highly significant when tested by permutation (Fig. 3; b = 3.943 ± 0.67; F1,82 = 34.88, P < 0.0001; adjusted R2 = 0.29).The PIC for global methylation did not depart significantly from normality (W = 0.97, P = 0.077; Shapiro–Wilk test), whereas the PIC log(Cvalue) did (W = 0.90, P < 0.0001).
The hierarchical partition of variance in global methylation indicated that genera (39.1%) and order (19.9%) were the taxonomic levels accounting for most variance in the present species sample. Variance among species within the same genus was substantially lower (5.3%) and the error component attributable to individual replicates was negligible (1.96%).
Geographical and ecological correlates
The two geographical regions largely overlapped in the range of global methylation (Fig. 4). They differed, however, in the average value (F1,241 = 14.57, P < 0.0001), with genomes of tropical species being on average less methylated than those of Mediterranean ones (Fig. 4). This difference still held within woody (F1,139 = 7.13, P = 0.009) and nonwoody species (F1,83 = 6.57, P = 0.012) when these two growth forms were analyzed separately. A significantly lower average methylation in the genomes of tropical species also was found when the analysis was restricted to plants from the ecologically comparable, coastal vegetation type (19.1 ± 1.2 vs 22.6 ± 1.4 for tropical and Mediterranean coastal communities, respectively; F1,43 = 5.77, P = 0.02), which likewise supports the existence of a geographical effect on genome methylation levels.
Genomes of woody plants tended to be less methylated than those of nonwoody ones (F1,241 = 3.62, P = 0.058; Fig. 5a) both in tropical and Mediterranean species, and deciduous woody plants have genomes with lower methylation than evergreen ones (F1,102 = 3.23, P = 0.07; Fig 5b), although the difference was significant only in the Mediterranean plants (F1,60 = 7.35, P = 0.0087). Furthermore, geographically widespread species had genomes significantly less methylated than ORhad genomes with significantly lower methylation than species with comparatively restricted distributions (F2,240 = 6.35, P = 0.002; Fig. 5c).
Global methylation varied widely among species in all plant communities studied (Fig. 6). No evidence was found for significant variation in global methylation across Mediterranean communities (F3,130 = 0.93, P = 0.43). Besides, overall heterogeneity among the tropical communities studied (F3,88 = 2.88, P = 0.04) reflected an unforeseen difference between the intermediate positions in the environmental gradient analyzed, the lowland dry thorny and dry forests, that exhibited, respectively, the maximum (21.52 ± 1.48%) and minimum (16.05 ± 1.5%) average global methylation within this region (Fig. 6).
Discussion
The study of the role of epigenetics in plant adaptation and evolution should ideally be addressed using a diverse range of tools, approaches and scales, running from the dissection and experimental manipulation of molecular mechanisms to comparative analyses involving taxa growing in contrasting environments. The molecular understanding of epigenetic regulation, significance of DNA methylation and connections between different epigenetic traits in model plants, has considerably improved in the last decades (Diez et al., 2014; Matzke et al., 2015). Furthermore, recent comparative analyses of model and commercial plant species whose methylomes have been fully sequenced suggest that genomic distribution of methylated cytosines, regulation of chromatin modifications, and even the activity of different DNA methylation pathways, vary among distant lineages (Niederhuth et al., 2016; Takuno et al., 2016; Vidalis et al., 2016) and even closely related species (Seymour et al., 2014; Song & Cao, 2017). By contrast, next to nothing is known about the possible correlations between interspecific patterns of DNA methylation in wild plants, on the one hand, and plant features and ecological or geographical gradients, on the other. The present study constitutes a first attempt at identifying possible links among epigenome, phenotype and environment (Richards et al., 2017) in a broad sample of wild-growing species from diverse ecological scenarios. Results have revealed that the broad interspecific variation in global cytosine methylation was simultaneously related to phylogenetic patterns, environmental features and plant phenotypic traits, as discussed below.
Interspecific variation of global methylation, phylogenetic signal and correlated evolution with genome size
The present results confirmed that global methylation of plant genomes can vary 10-fold across angiosperm species, and that such variation is significantly associated with phylogenetic relatedness in a broader taxonomic and ecological sampling than addressed previously (Alonso et al., 2015; Vidalis et al., 2016). Tropical species in the present sample, however, exhibited a looser phylogenetic relationship suggesting lower strength of species differentiation strictly due to their divergence time (Brownian motion, BM) in this particular sample (Münkemüller et al., 2012). In addition, the analysis of local phylogenetic association revealed that the relationship between phylogenetic position and genome methylation tended to be stronger at some specific nodes located at different hierarchical levels within the tree, each of which might perhaps denote an evolutionary burst of rapid changes in genome-wide methylation.
The complex local phylogenetic association of global DNA methylation is consistent with the results of variance partitioning showing that most variance was explained by genus and order but not the intermediate family. Correlated evolution between global DNA methylation and genome size supported by Phylogenetically Independent Contrasts (PIC) analysis (Alonso et al., 2015; Vidalis et al., 2016) could in part explain some of these local patterns. For instance, Magnoliids have small genomes (Leitch & Leitch, 2013) and the three Magnoliid taxa studied here (Malmea depressa, Aristolochia baetica and A. paucinervis) had characteristically low genomic methylation levels, falling within the 5% percentile of the whole sample distribution. At the opposite extreme, most study taxa in monocots belonging to Liliales and Asparagales have big genomes (Leitch & Leitch, 2013), fell beyond the 95% percentile of global DNA methylation recorded (e.g. Pancratium maritimum, Polygonatum odoratum and Narcissus bulbocodium all have > 40% cytosines methylated), and exhibited highly significant local phylogenetic association. In eudicots, local phylogenetic association for either above or below average global methylation were found in the Rosiids clade (e.g. Fagales, Rosales, Myrtales) but not in the Asteriids. Some parallelism can be found with the evolutionary scenario of the plant-specific DNA methyltransferases recently proposed by Bewick et al. (2017) that supported two groups of chromomethylases in monocots/basal angiosperms, with several duplication events and specific losses within this clade. Further duplication, diversification and losses of genes encoding this methyl transferase family give rise to clade-specific, family-specific and species-specific chromomethylase combinations within eudicots (Bewick et al., 2017). Finally, there is a need to acknowledge the potential contribution of whole-genome duplication events, as intrinsic genomic features likely associated with variation in global methylation (Diez et al., 2014; Alonso et al., 2016a; Springer et al., 2016), and whose specific effect on generating interspecific variability in global methylation could not be evaluated here due to scarcity of documented polyploidy cases in the present dataset.
Environmental and phenotypic correlates of global methylation
The quantitative assessment of global DNA methylation with methods that do not require a full genome sequence provide an affordable method to explore world-wide biogeographical patterns of this epigenomic feature, as well as testing explicit ecological hypotheses, for example, whether plants growing in stressful environments, either at the global scale (e.g. arctic or desert areas) or within a certain region (e.g. edaphic islands or flooded lands), may exhibit contrasted levels of DNA methylation as has been found among populations at the intraspecific level (Lira-Medeiros et al., 2010; Medrano et al., 2014; Schulz et al., 2014). These possibilities are clearly exemplified by the results of the present investigation. Different spatial scales provided contrasting results regarding the association between environmental components and interspecific variation in global methylation. On the one hand, local variation in global methylation was extensive and overlapping in the eight study communities, suggesting that potential intercommunity differences are not related to the rainfall gradient studied and/or that species with contrasting epigenomic features can coexist, which is likely to be associated with substantial microenvironmental heterogeneity in communities of these seasonal environments. On the other hand, at the broadest geographical scale, genomes of tropical plants were, on average, less methylated than genomes of Mediterranean plants, the difference remaining significant when comparisons were restricted to the coastal vegetation type, or to woody and nonwoody growth forms. The present authors’ are not aware of previous reports of geographical divergence in epigenetic features of organisms with the exception of the negative relationship between global DNA methylation and environmental temperature shown for some ectothermic vertebrates (Varriale, 2014, and references therein). The analysis of plant lineages growing on different biomes and species with disjunct distributions (Villaverde et al., 2017) could help to clarify whether contrasting glacial histories or subsequent divergence of environmental conditions did significantly contribute to the observed geographical divergence in global methylation, on the implicit assumption that ecological history has some measurable effects on genomic methylation.
The magnitude of global DNA methylation was associated with woodiness, one of the main axes of variation of functional diversity of extant plants (Díaz et al., 2016). The genomes of woody species tended to be less methylated than those of nonwoody ones. This trend seemed independent of phylogenetic effects, because methylation differences were highly significant in the few genera with data available for both woody and nonwoody perennials (e.g. Viola, Hypericum). The frequently ancestral condition of woodiness in many angiosperm clades (Dodd et al., 1999), along with their reduced genome size (Ohri, 2005), support the tentative hypothesis that methylation-based epigenetic mechanisms might have contributed to key evolutionary transitions within the angiosperms (Fedoroff, 2012). It also was found that genomes of deciduous woody plants tended to be less methylated than those of evergreen species, although intrageneric comparisons revealed that this pattern was supported only in Quercus, not in Rhamnus, Viburnum and Pistacia, suggesting that this relationship is more labile and/or taxon-specific. Further analyses are needed before the possible association between genomic methylation and leaf phenology can be established.
One intriguing result of this study was the inverse relationship between genomic methylation and size of geographical range, with widespread species tending to have less-methylated genomes than those with narrower distributions. With due caution, because species’ range sizes were categorized using discrete classes rather than by direct quantification, this pattern suggests that reduced genomic methylation could perhaps favour the exploitation of a wider range of environmental conditions via enhanced epigenetic variability among populations or individuals (Richards et al., 2012), and/or that higher methylation could constrain the potential for colonization of new environments. These aspects deserve further study because of their potential implications for the success of biological invasions (Vogt, 2017; Hawes et al., 2018).
Conclusion
This broad taxonomic survey, including species from 37 orders and two distant geographical regions, supported the notion that global DNA cytosine methylation is a ‘epiphenotypic’ plant trait deserving consideration in ecological epigenetics studies of nonmodel plants. The considerable interspecific variation exhibited by this trait was found to be associated with evolutionary history, genome size, functional traits and range size. Results of analyses of local phylogenetic association further highlighted the importance of accounting for phylogeny when conducting epigenetic analyses of nonmodel plants in an interspecific context. Concurrent analyses of intra- and interspecific patterns of global methylation patterns in selected plant clades associated with contrasting environments or geographical gradients should foster improved understanding of the epigenetic contribution to evolutionary process.
Acknowledgements
We thank Pilar Bazaga, Laura Cabral, César Canché-Collí, Esmeralda López, Denis Marrufo and Rosalina Rodríguez for invaluable help during field and laboratory work; Filogonio May, Paulino Simá and José Luis Tapa for taxonomic identification of Mexican species; Rocío Astasio for assistance in obtaining permits, authorities and guide (Santos Uc) from Helen Moyers Biocultural Reserve (Kaxil-Kiuic) and Reserva Estatal El Palmar in Yucatán; Consejería de Medio Ambiente, Junta de Andalucía in Spain, for authorizing this work and providing invaluable facilities; and organisers and attendees of the 40th New Phytologist Symposium on ‘Plant epigenetics: from mechanisms to ecological relevance’ and three anonymous referees for discussion. Financial support was provided by grants SEV-2012-0262 (Estación Biológica de Doñana through the Severo Ochoa Program for Centres of Excellence), and CGL2013-43352-P and CGL2016-76605-P (Ministerio de Economía y Competitividad, Spanish Government).
Author contributions
CMH designed the research; MM, AC, VP-T, CMH and CA conducted field sampling; MM, CMH and CA contributed to data analysis; MM made the species information search and prepared all the figures; RP conducted all HPLC analyses and contributed to their interpretation; CA led the writing; all authors contributed to refine the manuscript.
