Genetic and epigenetic variations associated with adaptation to heterogeneous habitat conditions in a deciduous shrub

Abstract Environmentally induced phenotypic plasticity is thought to play an important role in the adaption of plant populations to heterogeneous habitat conditions, and yet the importance of epigenetic variation as a mechanism of adaptive plasticity in natural plant populations still merits further research. In this study, we investigated populations of Vitex negundo var. heterophylla (Chinese chastetree) from adjacent habitat types at seven sampling sites. Using several functional traits, we detected a significant differentiation between habitat types. With amplified fragment length polymorphisms (AFLP) and methylation‐sensitive AFLP (MSAP), we found relatively high levels of genetic and epigenetic diversity but very low genetic and epigenetic differences between habitats within sites. Bayesian clustering showed a remarkable habitat‐related differentiation and more genetic loci associated with the habitat type than epigenetic, suggesting that the adaptation to the habitat is genetically based. However, we did not find any significant correlation between genetic or epigenetic variation and habitat using simple and partial Mantel tests. Moreover, we found no correlation between genetic and ecologically relevant phenotypic variation and a significant correlation between epigenetic and phenotypic variation. Although we did not find any direct relationship between epigenetic variation and habitat environment, our findings suggest that epigenetic variation may complement genetic variation as a source of functional phenotypic diversity associated with adaptation to the heterogeneous habitat in natural plant populations.

Although there are several epigenetic mechanisms, including chemical modification of DNA and histones, position effects and interference by small noncoding RNAs (Richards, 2011), DNA methylation of cytosine is most the extensively studied epigenetic mechanism with important effects on ecologically relevant traits (Herrera & Bazaga, 2008, 2013Schrey et al., 2013;Wilschut et al., 2016;Xie et al., 2015). While DNA is predominantly methylated at CG sites, cytosine methylations in plants occurs throughout the genome in all sequence contexts (CG, CHG and CHH where H = A, C or T) (Law & Jacobsen, 2010) and could affect whether transposons are silenced and genes are expressed. Biologists have illuminated that the amount and pattern of DNA methylation in model plants is sensitive to various environmental stressors under laboratory conditions, and ecologists have focused on the variation in DNA methylation in wild populations to understand the role of DNA methylation in plant adaptation to real environmental stress in nature (Bossdorf, Richards, & Pigliucci, 2008;Kilvitis et al., 2014). The rapidly increasing number of publications has illustrated that variation in DNA methylation is correlated with herbivory in violets (Herrera & Bazaga, 2010), salinity in marsh perennials (Foust et al., 2016), artificial disturbance in Lavandula latifolia , metals in red maple (Kim, Im, & Nkongolo, 2016), and climate in Quercus lobata (Gugger, Fitz-Gibbon, PellEgrini, & Sork, 2016;Platt, Gugger, Pellegrini, & Sork, 2015).
There is a complex relationship between genetic and epigenetic variation in the wild. Epigenetic variants can be under genetic control (Slotkin & Martienssen, 2007), but environmental factors can also directly alter the epigenetic variation that may be inherited through meiosis over several generations (Jablonka & Raz, 2009). The pattern of epigenetic variation in wild populations contributing to phenotypic traits which cannot be explained by genetic variation may be a consequence of natural selection on pre-existing epigenetic variation, environmental induction of variable epigenetic variation, or both (Bossdorf et al., 2008;Verhoeven, vonHoldt, & Sork, 2016).
Vitex negundo var. heterophylla (Chinese chastetree, or five-leaved chaste tree), a native deciduous officinal shrub with a heterophylly leaf and entomophilous flower (Figure 1), is widely distributed in the hilly areas of North China (Hu et al., 2015). Vitex negundo var. heterophylla is a typical species that can survive in a broad range of habitats, including bush and understory, varying in an array of other biotic and abiotic parameters (Li, Yang, & Wu, 2008). Chinese chastetree in roadside and bush habitats is frequently cut by farmers as firewood or to prevent it from occupying farmlands. It can adapt to cardinal environmental factors through modification of a series of morphological and physiological characteristics (Du, Guo, Zhang, & Wang, 2010;Du et al., 2012). As a long long-lived perennial, V. negundo var. heterophylla may face variable environments through epigenetic processes (Bräutigam et al., 2013). Therefore, it provides an ideal study system to compare genetic and epigenetic differences in response to heterogeneous habitat conditions.
In this study, we measured several plant functional traits for individuals sampled from each plot to determine habitat differentiation in phenotype. We investigated the genetic variation using amplified fragment length polymorphism (AFLP) and epigenetic variation using methylation-sensitive AFLP (MSAP). These dominant markers provided a considerable number of anonymous loci, which generate powerful data for detecting the genetic and epigenetic structure across several heterogeneous habitats in a nonmodel species without a reference genome. We compared the patterns of genetic and epigenetic variation shaped by the habitat conditions to test the hypothesis that epigenetic variation plays a potential role independent of genetic variation in the adaptation of V. negundo var. heterophylla to the environment.

| Sampling design
We selected plots for this study from seven sites in the south hills of Shandong Province, eastern China, separated by approximately 20-150 km ( Figure 2). For each site, we collected samples from two adjacent plots with contrasting habitat conditions at a distance of <2 km (Table 1) in May 2016. At each plot, 13-15 widely spaced individuals were randomly chosen and marked with permanent tags. To avoid developmental epigenetic variation confusing the differences among heterogeneous environments, all leaf samples for molecular analyses were picked at the same position of the plant before the flowering phase. Young expanding leaves were collected from each plant, placed in paper envelopes, and dried immediately with silica gel.
This material was used for genomic DNA extraction, and additional expanding leaves were randomly collected to measure the leaf traits.

| Measurement of plant functional traits
For each sampled individual, plant maximum height (H), basal diameter (D), and number of resprouts (NR) were measured, and more than five fully expanded compound leaves were collected F I G U R E 1 The study organism Vitex negundo var. heterophylla, a deciduous shrub with entomophilous flowers and digitate leaves containing five lanceolate leaflets, sometimes three at different positions of the stem. Both leaf fresh weight (LFW) and leaf dry weight (LDW) were measured and leaf water content [LWC = (LFW − LDW)/LDW] was obtained. Leaflet area (Area) was determined with image processing program ImageJ (Abràmoff, Magalhães, & Ram, 2004), and the specific leaf area (SLA = LA/ LDW) was calculated. The number of resprouts (NR) was divided into three levels to assess the environmental disturbance (Pérez-Harguindeguy et al., 2013): level 0 for one stem, level 1 for more than one but <5 stems, and level 2 for more than five stems. Finally, the coefficient of height-diameter allometry (HDA) was calculated for each population using function MSA(D~H, log = "xy") in R package smart (Warton, Duursma, Falster, & Taskinen, 2012).

| AFLP and MS-AFLP protocol
We investigated a total of 195 individuals for genetic and epigenetic variation with AFLP and MSAP using the same DNA sample for each individual. Total genomic DNA was extracted from dried leaf tissue according to the cetyl trimethylammonium bromide (CTAB) method (Doyle & Doyle, 1987) with some modifications. PVP 40,000 was used to improve DNA yield and quality. DNA was quantified with both 0.8% agarose gels and microscopic spectrophotometry.
The protocol for MSAP was adapted from a standard AFLP (Vos et al. 1995), replacing the MseI enzyme in two separate runs with the methylation-sensitive enzymes HpaII and MspI using appropriate  and HpaII/MspI + TA) analyses (Table 2), of which the EcoRI + 3 primers were 5′-end labeled with FAM dye.
The final selective PCR productions were separated and visualized on an ABI3730XL DNA capillary sequencer (Applied Biosystems, Foster City, USA) with Rox-500 internal size standard (Applied Biosystems) at the Shandong Academy of Agricultural Sciences. We used PEAKSCANNER v1.0 (Applied Biosystems) to analyze the AFLP and MSAP profiles. Binning of fragments was performed using a peak height threshold of 50 relative fluorescence units and a minimal size of 50 base pairs. Peak height data were scored with a binary code, zero for band absent, and one for band present. For every polymorphic locus, each allele must exist in more than two individuals (>1% of all samples).

| Data analysis
To identify how many different genetic groups are represented in our collection regardless of the geographic sampling location, we performed Bayesian clustering of genetic data using STRUCTURE v2.3.4 (Falush, Stephens, & Pritchard, 2007;Pritchard, Stephens, & Donnelly, 2000). Structure estimates the number of groups (K) present among individuals and assigns individuals to each K using Bayesian modeling. We tested 20 populations (k = 1-20), which are more than the maximum anticipated based on sampling location, with 20 runs at each k. We used both the log probability of observing the data (Ln Pr(x|k)) method of Structure and Delta K (Evanno, Regnaut, & Goudet, 2005) with the online program STRUCTURE HARVESTER for visualizing STRUCTURE output and implementing (Earl & vonHoldt, 2012), which determines the number of populations that best fit the data.
We incorporated sampling locations in our model to assist in detecting sampling weak differentiation. We performed clustering with the admixture model, 30,000 burn-in steps, 100,000 postburn-in steps, and allowed correlated allele frequencies. We assigned individuals to groups based on the highest q-value.
All basic statistical analyses were carried out using the R environment. The MSAP profiles were analyzed with the R script MSAP_calc  To estimate the amount of epigenetic and genetic variation, Shannon's diversity index (H) and percentage of polymorphisms (PPL) for each population were calculated with GENALEX 6.5 (Peakall & Smouse, 2006). We also used GENALEX to estimate genetic and epigenetic differentiation using hierarchical analysis of molecular variation (AMOVA) to determine whether spatial location structured genetic or epigenetic differences by comparing the genetic and epigenetic variation among sites (ФRT), among populations (i.e., habitats) within sites (ФPR), and within populations (ФPT). We used 9,999 permutations to estimate statistical significance and an initial alpha of 0.05. Moreover, we used generalized linear models (GLM) for each loci to determine the specific AFLP and MSAP loci correlated with the habitat type using function glm in R.
We analyzed the correlation between genetic variation, epigenetic variation, and habitat by performing Mantel and partial Mantel tests using zt software (Bonnet & Van de Peer, 2002 and HDA, performed using function dist (method = "euclidean") in R software. Both strategies make the assumption that differences between habitats will be essentially the same magnitude regardless of individual population differences. In all cases, we used the Euclidean genetic and epigenetic distance matrices generated by GENALEX.
We also constructed Euclidian geographic distances and Nei unbiased genetic and epigenetic distances to test the roles of major role for population differentiation. As simple and partial Mantel tests have been questioned for a number of drawbacks (Bradburd, Ralph, & Coop, 2013;Guillot & Rousset, 2013;Legendre, Fortin, & Borcard, 2015), we applied multiple matrix regression with randomization (MMRR) (Wang, 2013) as an alternative approach to Mantel procedures. Computations were implemented in with the MMRR function using 9,999 permutations.

| Differences in plant functional traits
At the population level, we found significant variation for D [paired t- and SLA [paired t-test: t(6) = 2.96, p = .025], and similar patterns of disturbance between two populations were revealed among sites other than Dafotou and Yaoxiang (Figure 3a). At the individual level, significant differences between two habitats were detected for LWC and SLA (Figure 3b) at four sites: Dafotou, Yaoxiang, Liantai, and Taishan (Table 3)

| Comparison between genetic and epigenetic variation
Simple Mantel tests showed significant correlations between genetic and epigenetic variation across all populations (Table 6). Both simple Mantel tests (Table 6) and MMRR (Table 7) showed that phenotypic variation was correlated with epigenetic variation, but MMRR showed a possible association between genetic differentiation and habitat types (r = .081, p = .0485). We also detected the geographic structure of genetic and epigenetic variation. At the level of population, we found a lower geographic regression slope for epigenetic structure than genetic structure using the Mantel method (β AFLP = 0.0234, β MSAP-u = 0.0042, β MSAP-m = 0.0062; Figure 5). The above results obtained from u-MSAP and m-MSAP markers separately showed the same trends (Table 6). MMRR also revealed the analogous genetic and epigenetic spatial structure (Table 8).
Partial Mantel tests always displayed a significant correlation (r > .10, p < .05 for u-MSAP; r ≥ .05, p < .05 for m-MSAP) between genetic and epigenetic variation when we controlled for the additional factor ( Figure 6). However, we detected no relationship between habitat type and genetic variation when we controlled for epigenetic variation and no relationship between habitat type and genetic variation when we controlled for epigenetic variation.
Moreover, we did not find a significant correlation between phenotypic and genetic variation when we controlled for epigenetic variation, but we found a significant correlation (r = .041, p = .033 for u-MSAP; r = .059, p < .001 for m-MSAP) between phenotypic and epigenetic variation when we controlled for genetic variation.
Additionally, we detected the spatial structure of genetic (r = .058, p = .002 with controlled u-MSAP; r = .058, p = .005 with controlled T A B L E 4 Values of diversity (h) and percent polymorphisms (%P) of populations

| DISCUSSION
Consistent with the previous studies of V. negundo (Su et al., 2003;Zhang, Zheng, & Ge, 2007), we found high levels of genetic diversity and low population genetic differentiation both among and within sites. In particular, Bayesian clustering revealed a significant genetic differentiation shaped by habitat heterogeneity. We also found high levels of epigenetic diversity but little epigenetic structure among populations within site, and random levels of DNA methylation across all populations. This result is consistent with the previous findings in salt marsh perennials (Foust et al., 2016), suggesting that genetic variation may be more strongly structured by environment than epigenetic variation. Although the scoring approaches are largely consistent when comparing the levels of diversity among populations, T A B L E 7 Summary of multiple matrix regression analysis with randomization (MMRR) relating the phenotypic distance matrix with genetic and epigenetic distance matrices Gen, genetic variation using AFLP markers; Epi_u, epigenetic variation using u-MSAP; Epi_m, epigenetic variation using m-MSAP; Phe, phenotype. p values were calculated with 9,999 permutations.
T A B L E 8 Summary of multiple matrix regression analyses with randomization (MMRR) relating genetic, epigenetic, and phenotypic distance matrices with geographical and environmental distance matrices estimates of epigenetic diversity varied strongly (Schulz et al., 2013), When running simple and partial Mantel tests, we obtained the same pattern as for Borrichia frutescens (Foust et al., 2016), failing to detect any relationship between epigenetic and habitat environment using crude binary data to describe complex habitats, as we had a larger number and range of sampling plots than the previous studies (Foust et al., 2016;Kim et al., 2016;Robertson, Schrey, Shayter, Moss, & Richards, 2017;Schulz et al., 2014). Limited studies have been conducted of epigenetic differentiation in natural populations across heterogeneous habitat conditions, and there is an urgent need to develop or replace the binary method to comprehensively characterize the complex habitat conditions. Here, we tried to introduce the functional traits for representing the environment.
The selected traits, including leaf water content, specific leaf area, level of resprout, and height-diameter allometry, were plastic and could respond to the heterogeneous environment. For example, SLA, an indicator of ecophysiological characteristics with a strong phenotypic plasticity and substantial genetic effects (Scheepens, Frei, & Stöcklin, 2010), would increase in the understory habitat to optimize light harvesting (McIntyre & Strauss, 2014), as supported by our data. In addition, NR could be considered as a typical parameter F I G U R E 6 Genetic and epigenetic correlations to variation in habitat environment, functional phenotype, and geographic distance using partial Mantel tests with the Euclidean genetic and epigenetic distance matrices, habitat distance matrices, and trait distance and geographical distance matrices for all individuals across all populations. The correlations between genetic and epigenetic variation were calculated in separate partial Mantel tests. Gen = genetic variation, Epi-u = Epigenetic variation based on MSAP-u, Epi-m = Epigenetic variation based on MSAP-M, Env = environmental variation, Phe = phenotypic variation, Geo = geographical distance, NS = not significant, r = correlation coefficient when significant and p-value of branching architecture for the shrub, which differs on the basis of browsing history, fire history, and other types of disturbance (Pérez-Harguindeguy et al., 2013), especially access to light, water stress, and the human-caused cutting for V. negundo var. heterophylla in our research region. We found phenotypic plasticity in response to controlled drought and shade treatment in our previous glasshouse experiments using V. negundo var. heterophylla (Du et al., 2010(Du et al., , 2012, and this field experiment also displayed adaptive plasticity in response to heterogeneous wild habitats. Nevertheless, assessing plant functional trait data separately may misrepresent the adaptive response to the complex habitat. Thus, we combined several typical fitness-related phenotypes to quantify the divergence of adaptions. Amazingly, we found 15 AFLP loci statistically correlated with habitat type, and the Bayesian clustering suggested a parallel genetic divergence that may result in microevolution from independent origins. This was also supported by MMRR showing potential effect of habitat environment on genetic differentiation. Although strong gene flow has homogenized the genetic differentiation across the genome, adaptive loci will remain due to similar pressure (Trucchi, Frajman, Haverkamp, Schönswetter, & Paun, 2017). Only several plastic functional traits were investigated, and very limited loci were detected by AFLP. Therefore, to test this hypothesis, it is necessary to find the adaptive phenotypes and adopt next-generation sequencing methods. We did expect epigenetic differentiation between open and understory populations, but we only found one locus showing differentiation due to habitat type. Moreover, it remains a question whether this locus plays an independent role from genetic variation in adaptation to diverse habitats. It is possible either that the understory environment did not change any epigenetic variation between habitats, or that any existing epigenetic signature was too weak to be detected given the high epigenetic diversity between individuals among all populations.
The epigenetic mechanism may be restricted when natural plant populations endure some discrete human-caused disturbance, such as heavy metal pollution (Kim et al., 2016), experimental disturbance , and oil spills (Robertson et al., 2017).
According to our field investigations and planting experiments (Du et al., 2012) use statistical approaches to uncover patterns of epigenetic variation that are not predictable from patterns of genetic variation (Foust et al., 2016;Gugger et al., 2016;Herrera, Medrano, & Bazaga, 2016).
Simple and partial Mantel tests were proposed to combine the analyses of genetic and epigenetic variation and present the correlation between genetic or epigenetic variation and habitat while controlling for the correlation between genetic and epigenetic variation (Foust et al., 2016). A controlled planting experiment investigating epigenetic variation in response to warming failed to find significant correlations between epigenetic variation and phenotype or habitat using Mantel tests (Nicotra et al., 2015). In contrast, we found a weak but significant correlation between epigenetic variation and adaptive phenotype, with no significant correlation between genetic and adaptive phenotypes.
The results indicated that epigenetic mechanisms may play a more important role in adaptive plasticity than genetic variation in some scenarios.
Finding direct casual links between epigenetic and variation and ecological phenotype is a key challenge in the study of epigenetic adaptation with nonmodel species . A glasshouse experiment with various epigenetic recombinant inbred lines of Arabidopsis thaliana concluded that heritable epigenetic variation can cause substantial variation in ecologically important plant traits (Zhang et al., 2013). To investigate the relationship between epigenetic variation and functional plant diversity in wild populations, genetic and epigenetic marker-trait association analyses for 20 functional traits in a perennial herb were conducted, showing that more MSAP markers than AFLP involved in significant association (Medrano et al., 2014). Due to the limited loci and phenotypes in our study, we did not perform a genotype-phenotype association analysis. Bayesian clustering and GLM revealed the genetic differentiation between habitats, but Mantel tests and MMRR indicated a significant correlation between epigenetic and phenotypic variation. This contradictory result suggests that the key stable phenotypes associated with habitat type were not considered in our study, as we only measured several plastic functional traits varying not only between habitats but also among all sites. However, our results added the indirect evidence that epigenetic variation can serve as an important source of intraspecific functional diversity in nature.
According to the illustration of genetic and epigenetic isolationby-distance scenarios along with hypothesized causes , our data revealed the moderate recent epigenetic variation between generations in natural populations of V. negundo var. heterophylla. Consistent with the results from a perennial herb (Herrera, Medrano, & Bazaga, 2017), geographic distance explained genetic differentiation better than epigenetic differentiation, and epigenetic variation contributed to the divergence in functional traits in the perennial shrub. It is necessary to conduct transplantation or common garden experiments with offspring from populations in different habitats to distinguish the contribution of induced and inherited epigenetic variation to adaptive phenotypes. As the MSAP marker only provides limited anonymous loci that are difficult to link to functional genomic elements or phenotype, a reduced representation bisulfite sequencing approach based on next-generation sequencing methods is the next level of epigenetic study in natural populations (Gugger et al., 2016;Platt et al., 2015;Robertson & Richards, 2015;Trucchi et al., 2016).

| CONCLUSION
Our study used functional traits, AFLP, and MSAP to analyze phenotypic, genetic, and epigenetic variation in natural populations of V. negundo var. heterophylla from an extensive variety of habitat environments. Two scoring approaches for MSAP marker data were conducted, and the results were consistent in most cases. The analysis demonstrated significant habitat-related adaptation in phenotypic and genetic differentiation, suggesting an evident process of natural selection. However, we did not find a direct correlation between functional phenotypic and genetic variation, and Mantel tests and MMRR revealed a significant relationship between epigenetic and genetic or phenotypic variation. This result ultimately implied a potential intermediary role of epigenetic mechanisms in the adaption to heterogeneous habitats.

ACKNOWLEDGMENTS
We want to thank Cui Wenqiang, Liu Shuna, Chen Bin, and Ma Haoran for their assistance with fieldwork. We are also grateful to Dr. Cai Yunfei for his guidance in the molecular biology experiments.

CONFLICT OF INTEREST
None declared.

AUTHOR CONTRIBUTION
LL, WG, and ND conceived the idea and designed the study. LL, CP, and ND contributed to the field survey. CP and LL performed the molecular laboratory work. LL, XG, and WG analyzed and interpreted the data. LL drafted the manuscript and all authors participated in manuscript modifications and gave final approval for publication.