Methylation divergence of invasive Ciona ascidians: Significant population structure and local environmental influence

Abstract The geographical expansion of invasive species usually leads to temporary and/or permanent changes at multiple levels (genetics, epigenetics, gene expression, etc.) to acclimatize to abiotic and/or biotic stresses in novel environments. Epigenetic variation such as DNA methylation is often involved in response to diverse local environments, thus representing one crucial mechanism to promote invasion success. However, evidence is scant on the potential role of DNA methylation variation in rapid environmental response and invasion success during biological invasions. In particular, DNA methylation patterns and possible contributions of varied environmental factors to methylation differentiation have been largely unknown in many invaders, especially for invasive species in marine systems where extremely complex interactions exist between species and surrounding environments. Using the methylation‐sensitive amplification polymorphism (MSAP) technique, here we investigated population methylation structure at the genome level in two highly invasive model ascidians, Ciona robusta and C. intestinalis, collected from habitats with varied environmental factors such as temperature and salinity. We found high intrapopulation methylation diversity and significant population methylation differentiation in both species. Multiple analyses, such as variation partitioning analysis, showed that both genetic variation and environmental factors contributed to the observed DNA methylation variation. Further analyses found that 24 and 20 subepiloci were associated with temperature and/or salinity in C. robusta and C. intestinalis, respectively. All these results clearly showed significant methylation divergence among populations of both invasive ascidians, and varied local environmental factors, as well as genetic variation, were responsible for the observed DNA methylation patterns. The consistent findings in both species here suggest that DNA methylation, coupled with genetic variation, may facilitate local environmental adaptation during biological invasions, and DNA methylation variation molded by local environments may contribute to invasion success.

However, evidence is scant on the potential role of DNA methylation variation in rapid environmental response and invasion success during biological invasions. In particular, DNA methylation patterns and possible contributions of varied environmental factors to methylation differentiation have been largely unknown in many invaders, especially for invasive species in marine systems where extremely complex interactions exist between species and surrounding environments. Using the methylation-sensitive amplification polymorphism (MSAP) technique, here we investigated population methylation structure at the genome level in two highly invasive model ascidians, Ciona robusta and C. intestinalis, collected from habitats with varied environmental factors such as temperature and salinity. We found high intrapopulation methylation diversity and significant population methylation differentiation in both species. Multiple analyses, such as variation partitioning analysis, showed that both genetic variation and environmental factors contributed to the observed DNA methylation variation. Further analyses found that 24 and 20 subepiloci were associated with temperature and/or salinity in C. robusta and C. intestinalis, respectively. All these results clearly showed significant methylation divergence among populations of both invasive ascidians, and varied local environmental factors, as well as genetic variation, were responsible for the observed DNA methylation patterns. The consistent findings in both species here suggest that DNA methylation, coupled with genetic variation, may facilitate local environmental adaptation during biological invasions, and DNA methylation variation molded by local environments may contribute to invasion success.

K E Y W O R D S
biological invasion, DNA methylation, methylation divergence, methylation-sensitive amplified polymorphism, tunicate
Epigenetic modification represent one important rapid mechanism for invasion success, as epigenetic modifications can influence phenotypes by regulating gene expression without changing the underlying DNA sequences (Bossdorf, Richards, & Pigliucci, 2008;Hawes et al., 2018;Rapp & Wendel, 2005). Although epigenetic variation can be predicted from its genotypic contexts or specific genetic polymorphisms elsewhere in the genome, still part of epigenetic variation can be independent from genetic control and is involved in response to environmental challenges as an autonomous system (Richards, 2006(Richards, , 2008. Various types of epigenetic modifications, including modifications of DNA, histones and other chromosomal proteins, and the generation of extrachromosomal regulatory small RNAs and noncoding RNAs, have been found in many eukaryotes (Richards, 2008). Among these epigenetic modifications, DNA methylation, the addition of a methyl group to a specific base (mostly cytosine for eukaryotes), is one of the most intensively studied for its ubiquity in eukaryotes and essential biological functions (Jones, 2012;Suzuki & Bird, 2008). DNA methylation is involved in numerous biological processes, such as blocking transcription initiation, controlling transcription elongation, silencing transposon, inactivating X chromosome, and imprinting on genes (Duncan, Gluckman, & Dearden, 2014;Jones, 2012). In addition, studies of DNA methylation have been advanced by the development of new techniques, such as bisulfite conversion-based methods (e.g., whole genome bisulfite sequencing-WGBS; reduced representation bisulfite sequencing-RRBS), immunoprecipitation-based methods (e.g., methylated DNA immune-precipitation sequencing-MeDIP-seq; methylated DNA-binding domain sequencing-MBD-seq), and restriction enzyme-based methods (e.g., methylation-sensitive amplified polymorphism-MSAP; Methyl-seq) (Suzuki & Bird, 2008;Wang et al., 2015). Although bisulfite conversion-based sequencing techniques have been considered as promising tools for detecting methylation modifications, these methods are limited in studies at the population level (i.e., large-scale population epigenomics studies), mainly owing to high experimental and computational costs as well as high-quality reference genomes sequenced and assembled for genomewide analyses (Schulz, Eckstein, & Durka, 2013;Suzuki & Bird, 2008). Instead, the methylation-sensitive amplification polymorphism (MSAP) technique, which can detect different methylation status at random restriction sites over a genome, is favorable as it allows to evaluate genomic methylation patterns of wild populations and identify associations between environmental conditions, phenotypic traits, and methylation status with low cost (Alonso, Pérez, Bazaga, Medrano, & Herrera, 2016). Thus, MSAP is the most widely used method in ecological epigenetics studies at the population level in a wide range taxa (Schulz et al., 2013).
During the process of local accommodation and adaptation, the most important feature of DNA methylation is its environmental effects (Verhoeven, VonHoldt, & Sork, 2016 (Dowen et al., 2012;Huang et al., 2017;Morán, Marco-Rius, Megías, Covelo-Soto, & Pérez-Figueroa, 2013;Platt, Gugger, Pellegrini, & Sork, 2015;Verhoeven, Jansen, van Dijk, & Biere, 2010). In the wild, DNA methylation differentiation among populations has been frequently observed in different environments (Gugger, Fitz-Gibbon, Pellegrini, & Sork, 2016;Paun et al., 2010;Wenzel & Piertney, 2014). In addition, maintenance and removal of DNA methylation are enzymatically mediated, thus the state of DNA methylation can change rapidly throughout the whole genome (Law & Jacobsen, 2010). Furthermore, changes in DNA methylation states can often cause important phenotypic consequences, leading to rapid accommodation and even adaptation to different local environments (Gao, Geng, Li, Chen, & Yang, 2010;Lea, Altmann, Alberts, & Tung, 2016;Pu & Zhan, 2017). Therefore, DNA methylation can integrate environmental signals into genomes rapidly and thereby modify gene expression and/or phenotypic variation to adapt to environmental changes/stresses. DNA methylation variation is expected to be involved in the process of accommodation/adaptation to heterogeneous environments for invasive species (Hawes et al., 2018 and references therein;Pu & Zhan, 2017). Recent studies in invasive species showed that methylation diversity was much greater than genetic diversity in introduced populations of Japanese knotweed Fallopia species (Richards, Schrey, & Pigliucci, 2012) and Alligator weed Alternanthera philoxeroides (Gao et al., 2010). Liebl, Schrey, Richards, and Martin (2013) inferred that high methylation diversity could compensate for reduced genetic diversity during the invasion process of the house sparrow Passer domesticus. Interestingly, introduced populations of bluegrass Poa annua tended to have higher levels of methylation diversity than native populations (Chwedorzewska & Bednarek, 2012), and studies of the pygmy mussel (Xenostrobus securis) and tubeworm (Ficopomatus enigmaticus) found that recently introduced populations seemed to be less methylated in comparison with older introduced populations (Ardura, Zaiko, Morán, Planes, & Garcia-Vazquez, 2017). Importantly, studies of invasive plants have also detected the association between DNA methylation variation and different habitats (Gao et al., 2010;Richards et al., 2012), and the authors further identified habitat-related methylated loci (Richards et al., 2012). In order to further analyze possible contributions of local environments to methylation differentiation, studies assessed the degree of methylation differentiation among populations sampled from different environments. Interestingly, the level of methylation differentiation highly varied among invasive species, for example, strong methylation divergence was observed among populations of bluegrass (differentiation = 0.5; Chwedorzewska & Bednarek, 2012) and Japanese knotweed (differentiation = 0.5-0.8; Richards et al., 2012), while low methylation differentiation was found in introduced house sparrow populations in Kenya (differentiation = 0.004; Liebl et al., 2013). Collectively, it still remains largely unexplored on methylation patterns and possible contributions of varied environmental factors to methylation differentiation during range expansions in many invasive species.
Here, we used invasive model ascidians C. robusta (=C. intestinalis sp. A, Brunetti et al., 2015) and C. intestinalis (=C. intestinalis sp. B, Brunetti et al., 2015) to investigate population methylation structures and possible responsible environmental factors for observed patterns in diverse habitats. C. robusta and C. intestinalis, two morphologically similar but genetically distinct species of the C. intestinalis complex, have widely invaded temperate and warm-temperate coastal zones over the past century (Carver, Mallet, & Vercaemer, 2006;Zhan, Macisaac, & Cristescu, 2010;Zhan et al., 2015). Due to ambiguous taxonomy within the genus Ciona, the native/invaded ranges and invasion histories of both C. robusta and C. intestinalis remain uncertain (Zhan et al., 2015 and references therein). C. robusta, which is generally considered as an East Asian native, has invaded coasts of Mediterranean Sea and many regions of Atlantic and Pacific coasts, while C. intestinalis, likely a native of Europe, has widely col- peaks per year (Carver et al., 2006 and reference therein). The broad geographical distribution, as well as varied environments in their native and/or invaded habitats, reflects the wide range environmental tolerance of both invasive species, especially for temperature and salinity (Dybern, 1967;Therriault & Herborg, 2008;Zhan et al., 2015).
Both C. robusta and C. intestinalis have a short pelagic larval phase (1-5 days) and a sessile adult stage, suggesting that natural dispersal can only occur over relatively limited geographical ranges . Therefore, the wide range expansion of both invasive species, particularly at the regional and continental scales, mainly is owing to human-mediated pathways, which results in sudden shifts of habitats with varied environments. As such, these two ascidians provide good models to study mechanisms of rapid microevolution during biological invasions Pu & Zhan, 2017;Zhan et al., 2015). Our previous study revealed significant genetic differentiation among C. robusta populations collected globally by using genomewide gene-associated microsatellites, and genetic signatures and loci under selection were associated with varied local environmental conditions . Furthermore, when we analyzed the DNA methylation variation of salinity and temperature-related genes using bisulfite sequencing at the population level, we found significant variation of DNA methylation among populations (Pu & Zhan, 2017). Interestingly, frequencies of several CpG loci were significantly correlated with local environmental factors (Pu & Zhan, 2017). However, Pu and Zhan's (2017) only focused on five genes which were putatively responsible for changes of temperature and salinity, and it remains largely unknown how environmental changes associated with habitat invasions shape methylation divergence among Ciona populations in different environments at the genome level.
In this study, we investigated the DNA methylation variation in C. robusta and C. intestinalis populations using the MSAP technique.
We analyzed the DNA methylation patterns for populations of both species collected from different continents with varied local environments. In order to interpret the significant population methylation differentiation among populations, the possible contributions of genetic variation and two crucial environmental factors in marine ecosystems (i.e., temperature and salinity) were further tested.

| Sample collection
To cover the substantial environmental gradients in habitats of C. robusta and C. intestinalis, sampling sites were selected based on our former studies (Zhan et al., 2010(Zhan et al., , 2012) across four continents (Asia, Africa, Europe, and Oceania) for C. robusta and two continents (Europe and North America) for C. intestinalis (Table 1; All collected adult specimens were preserved in 100% ethanol at 4°C before analyses. All specimens were identified and confirmed to the species level by using one mitochondrial DNA fragment, cytochrome c oxidase subunit 3-NADH dehydrogenase subunit 1 (COX3-ND1; Zhan et al., 2010).

| Methylation-sensitive amplification polymorphism assay
Total genomic DNA was isolated from approximately 50 mg of siphon tissues following the proteinase K method (Waters, Dijkstra, & Wallis, 2000). The quality and quantity of DNA were measured using Nanodrop 2000 spectrophotometer (Thermo Scientific).
Total genomic DNA (300 ng) was digested at 37°C for 3 hr in two parallel reactions using 5 U of EcoRI-HF and 5 U of either MspI or HpaII (New England Biolabs) in a final volume of 10 μl. MspI and HpaII can recognize and cleave the same sequence 5′-CCGG-3′, but their susceptibility is different depending on the methylation state of cytosines at restriction sites (Schulz et al., 2013). After digestion, enzymes were heat-deactivated at 80°C for 10 min. EcoRI_ Based on the number of amplified fragments, five primer pairs ( were excluded from further analyses. Replicated samples (6% of the total) starting from DNA extraction were included in all steps to test the reproducibility of MSAP assay (Bonin et al., 2004). Epiloci with more than two mismatches across the replicated samples were discarded from our datasets. The error rate based on replicated samples was estimated as 5.71% and 5.66% for C. robusta and C. intestinalis, respectively.

| Data analyses
Individuals were scored based on the presence (as "1") or absence ; type II and type III scored as "1," others scored as "0," i.e., methylated epiloci). Thus, the multistate raw data matrices were recorded into two binary data sets: dataset U and dataset M following the "Mixed Scoring 1" approach (Figure 2; Schulz et al., 2013). This approach can temper the confusion of data interpretation between types II and III, which may be caused by the methylation variation at two closely spaced CCGG sites (Fulneček & Kovařík, 2014).
Intrapopulation methylation diversity was characterized by the percentage of private loci, percentage of polymorphic loci, and Shannon's information index for both u-and m-subepiloci. Using GENALEX version 6.5 (Peakall & Smouse, 2006), population methylation differentiation (Φ PT ) was determined based on both u-and m-subepiloci for all intraspecific population pairs in both species by pairwise analysis of molecular variance (AMOVA) with 999 randomizations (Excoffier, Smouse, & Quattro, 1992). Pairwise distances between intraspecific individuals were computed by Euclidean measure for both u-and m-subepiloci with GENALEX version 6.5.
To visualize the methylation differentiation among populations of C. robusta and C. intestinalis, principal coordinates analyses (PCoA) of u-and m-subepiloci were conducted separately with the individual-by-individual methylation distance matrices using GENALEX version 6.5.
As genetic analyses based on genomewide microsatellites have been conducted for C. robusta , we used C. robusta as an example to evaluate the relationship between genetic and methylation differentiation for populations using Genodive (Meirmans & Van Tienderen, 2004) with 99,999 permutations.
The genetic differentiation matrix was obtained based on pairwise

| Intrapopulation methylation diversity
We generated 332 and 456 epiloci for C. robusta and C.  (Table 4). Interestingly, Shannon index of m-subepiloci was much higher than that of u-subepiloci in both C. robusta and C. intestinalis (Table 4) (Table 4).

| Interpopulation methylation differentiation
For C. intestinalis, higher Φ PT values in U profile than those in M profile were observed in all population pairs (Table 5B). Putative native populations from Europe (SL and SC) were highly and significantly differentiated from invasive populations collected from North America (HF, MR and YM) in both U (Φ PT = 0.133-0.216, p < 0.05; Table 5B) and M profiles (Φ PT = 0.053-0.083, p < 0.05; Table 5B). The differentiation among three Canadian invasive populations (HF, MR & YM) was relatively low but significant in both U (Φ PT = 0.018-0.069, p < 0.05; Table 5B) and M profiles (Φ PT = 0.010-0.023, p < 0.05; Table 5B). The highest methylation differentiation was detected between two adjacent European native populations SC and SL in the U profile (Φ PT = 0.224, p < 0.05; Table 5B), while in M profile, the highest value was found between an European native population SL and a Canadian invasive population HF (Φ PT = 0.083, p < 0.05; Table 5B).
When interpopulation methylation differentiation was assessed by pairwise AMOVA, significant methylation differentiation in C. robusta was observed between all population pairs in both U and M profiles (p < 0.05). Interestingly, overall pairwise Φ PT values in U profiles were higher than their corresponding values in M profiles (Table 5A).
The highest methylation differentiation for both U (Φ PT = 0.274;  Mantel tests showed significant correlation between genetic differentiation (Table 6) with methylation differentiation (Table 5A)

| Environment-related subepiloci
SAM analyses detected that 24 (4.0%) and 20 (2.4%) subepiloci were correlated with at least one of the six environmental factors for C. robusta and C. intestinalis, respectively (Table 9). BAYESCAN identified six and 10 outliers for C. robusta and C. intestinalis, respectively ( Figure 5; Table 9). Three environment-related subepiloci in C. robusta and two in C. intestinalis were also identified as outliers in BAYESCAN ( Figure 6; Table 9). The frequencies of environmentrelated subepiloci varied along with the environmental factors and showed significant difference among populations ( Figure 6). For C. robusta, seven subepiloci were temperature-related, of which five were correlated with the minimum temperature. We detected 20 salinity-related subepiloci in C. robusta, and 85% were related to the mean minimum salinity (Table 9). Furthermore, seven subepiloci were related to more than one salinity parameter, and three subepiloci were associated with both temperature-and salinity-related parameters. For C. intestinalis, nine subepiloci were temperaturerelated and 16 were correlated with salinity. Interestingly, 15 of 16 salinity-related epiloci in C. intestinalis were correlated with the mean maximum salinity, mean minimum salinity, and average salinity at the same time, and five of nine temperature-related subepiloci were also associated with salinity (Table 9).

| D ISCUSS I ON
In this study, we explored the population methylation patterns of C. robusta and C. intestinalis at the genome level. Our results revealed highly consistent patterns in both invasive species. We detected a high degree of intrapopulation methylated polymorphism for both invasive species. In addition, both C. robusta and C. intestinalis showed significant population methylation differentiation and TA B L E 5 Estimates of population methylation differentiation in two highly invasive ascidians Ciona robusta (A) and C. intestinalis (B). Above diagonal: pairwise Φ PT based on m-subepiloci, below diagonal: pairwise Φ PT based on u-subepiloci both genetic variation and varied local environments could partly explain the DNA methylation variation of C. robusta. We identified subepiloci presumably under selection and some of them were significantly correlated with environmental factors such as temperature and/or salinity. All results in our study suggest the putative contribution of DNA methylation variation to environmental response during biological invasions.

| Methylation diversity within populations
For C. robusta and C. intestinalis, a great proportion of subepiloci was polymorphic, especially for m-subepiloci (79.03%-92.40% in C. robusta, 81.39%-94.39% in C. intestinalis). The ratios of polymorphic loci were higher than those in other species, such as the tubeworm F. enigmaticus where only 50% of the methylation-susceptible loci were polymorphic (Ardura et al., 2017). A study on a rare floodplain herb Viola elatior showed that the polymorphism among populations ranged from 22.7% to 61.8% for u-subepiloci and varied from 30.6% to 61.8% for methylated subepiloci (Schulz, Eckstein, & Durka, 2014). Spontaneous epi-mutations, which generally occur as the consequence of incorrect replication of the methylation pattern, may be one important reason for the high level of methylation polymorphism, as epi-mutation rate can be several orders of magnitude higher than that of genetic mutations (Becker et al., 2011;van der Graaf et al., 2015;Hagmann et al., 2015;Schmitz et al., 2011). In addition, previous studies showed that Ciona ascidians are among the most genetically diverse animal investigated to date: The allelic polymorphism across the entire genome was 1.2% in C. robusta and the genomic average synonymous diversity was estimated to 0.05 per site in C. intestinalis (Dehal et al., 2002;Tsagkogeorga, Cahais, & Galtier, 2012). The high genetic variation can cause the absence of restriction enzyme sites, resulting in the high degree of methylated polymorphism.
The high level of intrapopulation methylation diversity can also be seen from the high values of Shannon index. Recent studies showed that environmental stresses could increase random methylation changes, leading to the high methylation divergence among individuals in the same populations Morán et al., 2013;Richards et al., 2012). In a salt-enriched experiment of brown trout (Salmo trutta L.), genomewide methylation changes were triggered after feeding on salt-enriched diets for 2 days, and the largest methylation differentiation was observed at the fourth day after treatment (Morán et al., 2013). Interestingly, our recent study showed that Ciona ascidians could respond to rapid changing environments in the methylation profile, such as after 1 hr of high temperature exposure or after 3 hr of low salinity challenge, and environmental changes increased intrapopulation methylation variation . Many studies clearly showed that environment-induced DNA methylation variation was heritable, thus such a type of DNA methylation variation can be relatively steadily  .

| Remarkable methylation differentiation among populations
In our study, significant population methylation differentiation was found among populations for both C. robusta and C. intestinalis. The remarkable methylation differentiation can be influenced by various factors, including genetic background, spontaneous epi-mutations, environmental induction and natural selection, and drift (Richards, 2006;Richards et al., 2017). For C. robusta, our results based on pRDA indicated that part of DNA methylation variation was explained by genetic variation. In addition to separate mechanisms shaping and maintaining genetic and DNA methylation variation, genetic and stable DNA methylation variation can be shaped simultaneously by the same local environmental pressures or neutral processes, leading to the correlation between genetic and methylation patterns (Preite et al., 2015). In the results of pRDA, more than 85% of the observed DNA methylation variation was not explained by available environmental factors and genetic variation, suggesting that the generation and/or maintenance of DNA methylation variation are complex processes and more variables including environmental and molecular ones may participate in such complex processes. More local environmental factors may influence DNA methylation variation at the population level, leading to a high degree of DNA methylation differentiation among populations collected from different environments. In addition, several recent studies suggest substantial genetic impacts on DNA methylation variation (Dubin et al., 2015;Roadmap Epigenetics Consortium et al., 2015). However, the contribution of genetic variation to DNA methylation variation was not as high as ex-     random DNA methylation variation has been recorded in natural populations to cope with unpredictably changing environments (Leung et al., 2016). Studies on Arabidopsis thaliana showed that the stochastic changes in environment-related methylation could be maintained over several generations and the accumulation of DNA methylation variance finally led to differentiation among populations (van der Graaf et al., 2015;Hagmann et al., 2015;Schmitz et al., 2011). Stress-related methylation differentiation has been documented in experimental populations (Gao et al., 2010;Huang et al., 2017;Morán et al., 2013;Verhoeven et al., 2010). The stability of stress-related DNA methylation variation varied from several hours to several generations (Herman, Spencer, Donohue, & Sultan, 2014;Huang et al., 2017). If the same advantageous phenotype arose every generation repeatedly by transient methylated modifications, such changes are expected to contribute to the interpopulation methylation differentiation in different environments (Schulz et al., 2014). The environment-related DNA methylation variation may have an effect to maximize individual fitness to local ecological conditions. Studies of three acroporid corals showed that DNA methylation variation may influence their tolerance to the thermal stress and ocean acidification (Dimond & Roberts, 2016). In wild baboons Papio cynocephalus, a large number of differentially methylated regions were found near metabolism-related genes, suggesting the epigenetic regulation for catering to different resource availability (Lea et al., 2016). For invasive species, recent works showed the correlation between methylation differentiation and diverse phenotypes in different invaded habitats, suggesting the potential role of DNA methylation variation in shaping different phenotypes to respond to varied environmental stresses (Gao et al., 2010;Richards et al., 2012). Despite that only a small proportion of DNA methylation variation was explained solely by collected environmental factors, fitness effects of some loci may be sufficient to make Ciona populations adapt to local environments. Our previous study also showed that a limited number of CpG loci and/or genes were involved in environmental adaptation Pu & Zhan, 2017).

| Environment-related epiloci
We detected that subepiloci of C. robusta (4.0%) and C. intestinalis (2.4%) significantly correlated with changes of environmental factors by SAM. Similarly, 14 epiloci (4.5%) were identified to significantly correlate with ecoclimatic variables in the study of Orchids Dactylorhiza (Paun et al., 2010). However, the proportion of environmental-related subepiloci in our study was much less than that in red grouse Lagopus lagopus scotica associated with gastrointestinal parasite load (13.6%; Wenzel & Piertney, 2014). Likewise, 12% of polymorphic epiloci in Spartina alterniflora were found to be correlated with oil exposure (Robertson, Schrey, Shayter, Moss, & Richards, 2017). The different ratios among studies may result from the difference of stresses and degree of stresses. In concert, findings in related studies indicated that environmental factors modulated the establishment and maintenance of methylation modifications to specific stress-related genes (Paun et al., 2010;Pu & Zhan, 2017).
Although epiloci related with temperature and salinity in C. robusta and C. intestinalis were identified, it was still unknown whether these environment-related epiloci were identical in certain groups responding to temperature or salinity in this study. This comparison cannot be clarified by using fragments acquired from MSAP, because fragments from different species with the same sizes may represent different sequences (Caballero et al., 2008), even they are the products of the same selective primers.
In this study, subepiloci with significant methylation differentiation among populations were identified by BAYESCAN, which may be involved in local adaptation (Foll & Gaggiotti, 2008). The overlap at subepiloci of SAM with BAYESCAN suggests that these environment-related subepiloci may contribute to the divergence of methylation pattern among populations. This hypothesis was confirmed by another specific study in populations of C. robusta, where we detected the correlation between methylation patterns of key genes and environmental factors (Pu & Zhan, 2017). Significant correlation was observed between methylation levels and water temperature at several CpG sites in heat shock protein 90 (HSP90), while one CpG site in Na + -K + -2Cl − cotransporter (NKCC) showed significant association with salinity (Pu & Zhan, 2017). As the technique of MSAP can only produce anonymous fragments with unknown genetic contexts, we cannot identify whether the outliers detected in these two studies play a consistent role in the process of response to environmental changes. All these results demonstrated that varied local environments may have important influence on the methylation patterns in populations. Therefore, we infer that DNA methylation may enhance invaders' adaptive capacity to different habitats by extending the flexibility of a genotype to respond differentially under variable environmental conditions.

| Technical issues
MSAP is a method based on PCR amplification of selected restriction fragments, which can be affected by both methylation modifications and genetic mutations of restriction sites. Therefore, MSAP data matrices inevitably include part of genetic information.
The ambiguous definition of type IV fragments in MSAP makes it difficult to distinguish between genetic and methylation variation. Although the polymorphism level of epiloci and population methylation differentiation may be affected by the mixed genetic variation, the role of local environmental factors in shaping DNA methylation variation can still be inferred from results of pRDA and SAM analyses. In order to reduce the confusion in data interpretation between type II and type III, the "Mixing Scoring 1" approach was used in our study to combine both type II (HPA+/ MSP−) and type III (HPA−/MSP+) as methylated epiloci. The use of this strategy resulted in more conservative results but some degree of loss of methylation variation among individuals. As MSAP has its technical disadvantages (Schulz et al., 2013), whole epigenome sequencing-based methods are needed in future studies to obtain more accurate epigenetic information.