Epigenetic patterns associated with an ascidian invasion: a comparison of closely related clades in their native and introduced ranges

Environmentally induced epigenetic modifications have been proposed as one mechanism underlying rapid adaptive evolution of invasive species. Didemnum vexillum is an invasive colonial ascidian that has established in many coastal waters worldwide. Phylogenetic analyses have revealed that D. vexillum populations consist of two distinct clades; clade B appears to be restricted to the native range (Japan), whereas clade A is found in many regions throughout the world, including New Zealand. The spread of D. vexillum clade A suggests that it might be intrinsically more invasive than clade B, despite low levels of genetic diversity compared to populations from the native region. This study investigated whether D. vexillum clade A exhibits epigenetic signatures (specifically differences in DNA methylation) associated with invasiveness. Global DNA methylation patterns were significantly different between introduced clade A colonies, and both clades A and B in the native range. Introduced colonies also showed a significant reduction in DNA methylation levels, which could be a mechanism for increasing phenotypic plasticity. High levels of DNA methylation diversity were maintained in the introduced population, despite reduced levels of genetic diversity, which may allow invasive populations to respond quickly to changes in new environments. Epigenetic changes induced during the invasion process could provide a means for rapid adaptation despite low levels of genetic variation in introduced populations.

A major threat to biodiversity worldwide is the introduction, establishment and spread of invasive species beyond their native range 1,2 . Due to the considerable ecological and economic impacts of biological invasions 3 , research has focused on the prevention, control and eradication of invasive species 4 . Molecular-based approaches are often utilized to characterize the ecological and evolutionary factors involved in the invasion process [5][6][7] , and are effective for estimations of propagule pressure 8 , early detection of introduced populations 9 , identifying the source of introductions 10 , and identifying genetic characteristics associated with successful invaders 11 . A major goal of invasion biology is to understand why certain species are successful invaders in order to facilitate assessment and management of risks 7 .
Phenotypic plasticity, a wide tolerance to environmental stress, rapid growth and reproduction, and a high level of competitive ability have all been proposed to facilitate invasiveness 12 . Higher genetic diversity within populations is also thought to aid in the successful colonization of new environments 13,14 . Specifically, low levels of genetic variation are known to increase the risk of inbreeding depression and limit the adaptive and evolutionary potential of a species 12 . Yet, invasive populations with little or no genetic variation can be extremely successful in new environments 15 , sometimes outcompeting native species 16 . Reduced competition from conspecifics 17,18 , the maintenance of a high level of quantitative trait diversity 19 , and the release from natural enemies 20 have been proposed as key factors determining invasion success when genetic diversity is low. More recently, epigenetic mechanisms (e.g., DNA methylation), which can increase both phenotypic plasticity and heritable variation, have been suggested to play a key role in invasion success [21][22][23] . Epigenetic mechanisms can alter a given genotype's influence on an organism's phenotype, without changes in the underlying DNA nucleotide sequences.
www.nature.com/scientificreports www.nature.com/scientificreports/ By increasing, decreasing, or silencing the activity of genes, epigenetic changes can lead to novel phenotypic variation in the absence of additional genetic diversity 24 . Despite this, few studies have investigated epigenetic patterns in invasive populations, which is an important step in understanding the role of epigenetic modifications in successful biological invasions.
A valuable approach for identifying genotypes and phenotypes associated with successful invasive species is to compare populations in their native and non-native habitats 25,26 . Molecular analyses have revealed that invasive species are often comprised of distinct sibling species or clades across their expanded and native ranges [27][28][29] , and in some cases, only certain clades become invasive 18,25,30 . Didemnum vexillum is an invasive marine colonial ascidian that has become established in temperate coastal waters worldwide, including several regions in New Zealand 31,32 , North America 33 , and Europe [34][35][36] . D. vexillum is believed to have been translocated throughout the world as fouling on aquaculture species and equipment, such as cultured Pacific oysters (Crassostrea gigas), and on ship hulls [37][38][39] . Additional to this well-defined history of spread, morphological 37 and molecular comparisons 18,40,41 have concluded that the native range of D. vexillum is likely to be the Northwest Pacific Ocean, including Japan. Phylogenetic analyses have revealed that D. vexillum populations consist of two distinct clades 18 . One clade (referred to as clade B) appears to be geographically restricted to D. vexillum's native region (Northwest Pacific Ocean), whereas colonies belonging to clade A have become established in temperate coastal areas around the world, including New Zealand. Clades A and B are morphologically identical 37 , and the level of the mitochondrial cytochrome c oxidase subunit I (COI) gene sequence divergence between the two clades is far below the comparable divergence values for other ascidian species, such as the recently classified cryptic species Ciona robusta and Ciona intestinalis (formally C. intestinalis type A and type B, respectively) 26,[42][43][44] . Therefore, evidence to date suggests that D. vexillum clades A and B represent an ongoing speciation event 18,26 . Diversification of clades A and B within the Northwest Pacific following a human-mediated introduction from another region is not consistent with known evolutionary rates at COI, as divergence of the A and B COI clade lineages is estimated to have occurred 1.5-2.5 million years ago, which pre-dates anthropogenic vectors 18 . The persistence of clade B's restricted distribution may indicate that it is intrinsically less invasive than clade A 18 . Introduced populations of clade A in New Zealand have low levels of genetic diversity compared to populations within the native region of Japan, indicating a founder effect 18 . Despite this, introduced D. vexillum populations have been extremely successful in their new environments, forming large colonies that can extend over significant areas, including commercially important aquaculture farms 39 .
The most well-studied epigenetic mechanism related to invasion success is DNA methylation, which generally refers to the addition of a methyl group to cytosine nucleotides to form 5-methylcytosine. DNA methylation can respond to environmental cues [45][46][47] , and has been associated with adaptation to different environmental conditions in invaded habitats 48,49 . Additionally, some introduced populations have higher DNA methylation diversity than populations in their native range despite reductions in genetic diversity due to founder effects 50 , and lower levels of global or genome-wide DNA methylation 51 , both of which have been associated with phenotypic variation and phenotypic plasticity [51][52][53][54] . Given the tendency of D. vexillum clade A to successfully invade multiple regions throughout the world, we hypothesized that (1) introduced populations of D. vexillum clade A would have high levels of DNA methylation diversity, despite low genetic diversity, (2) introduced populations of D. vexillum clade A would have lower levels of global DNA methylation compared to populations within the native region, and (3) populations would be differentially methylated when collected from sites with different environmental conditions. To test these hypotheses, we compared whole genome DNA methylation in D. vexillum colonies collected from New Zealand (clade A) with those from Japan (clade A and clade B), and between sites with different environmental conditions (differences in temperature and salinity). Samples used in the present study were a subset of samples previously collected and analysed by Smith, et al. 18 using mitochondrial cytochrome c oxidase subunit I (COI) partial coding sequences.

Materials and Methods
tissue sample collection. Didemnum vexillum tissue samples were collected between April 2008 and July 2009 from four sites in Japan (n = 39) and three sites in New Zealand (n = 35) by Smith, et al. 18 . Samples were collected in the summer months of the respective countries. Tissue samples were preserved in approximately 2.0 ml of 96% (v/v) ethanol and stored at -20 °C before being assigned to clade A or clade B by mitochondrial COI coding region analysis 18 . All colonies from New Zealand were clade A (n = 35), whereas samples collected from Japan were comprised of clade A (n = 16) and clade B (n = 23). In Japan, the distributions of clades A and B overlap at several sites (Fig. 1). The earliest worldwide record of D. vexillum was in Japan in 1926, but it was only identified as Didemnum sp. at the time 37 . In New Zealand, D. vexillum was first reported in 2001 in the North Island, in Whangamata Harbour 31 and was described initially as a native species to New Zealand. Two months later, D. vexillum was discovered in the South Island, in Queen Charlotte Sound 55 . It is thought to have been transported to the South Island as fouling the hull of a large barge that had earlier been anchored in Whangamata Harbour 55 . Due to the continued movement of vessels, D. vexillum subsequently spread throughout the Marlborough Sounds, including to Pelorus Sound 39 . environmental conditions. Site-specific environmental data (temperature and salinity) were collated from publicly available reports and websites (Table 1). Based on experiments identifying D. vexillum tolerance limits and the optimal temperature and salinity range for growth 46 , we allocated a temperature and salinity 'stress score' to each site, with 1 being the lowest stress score and 3 being the highest ( Table 2). Temperature and salinity stress scores related to tolerance ranges were as follows; Temperature stress: D. vexillum grows well in temperatures less than 23 °C, so a stress score of 1 was allocated to sites with maximum temperatures of <23 °C, D. vexillum growth is usually reduced at temperatures higher than 24 °C, so sites with a maximum temperature of 23-25 °C were given a stress score of 2, growth is further reduced above 25 °C, and 28 °C is close to the upper tolerance limit for www.nature.com/scientificreports www.nature.com/scientificreports/ D. vexillum 46,56 , so sites with a maximum summer temperature of 25-28 °C were given a stress score of 3; Salinity stress: D. vexillum can grow well in sanities above 26 PSU 46 , so sites with minimum salinity values ≥26 were given a stress score of 1, growth is generally reduced below 26 PSU, and ascidians are rarely found in salinities lower than 25 PSU 57 , so a stress score of 2 was given to sites with minimum salinity values of 23-<26 PSU, a stress score of 3 was allocated to sites that drop below 23 PSU, as 20 PSU is thought to be close to the low salinity limit for D. vexillum, with extended periods of exposure resulting in mortality 58,59 . We then combined the temperature and salinity stress score at each site to obtain a 'summative stress score' which was used for further analysis (see Table 2). This is a similar approach to the method described in Ardura et al. 60 .  18 . The date D. vexillum populations were first recorded in each location in New Zealand is indicated below each site and the number of colonies sampled is indicated in brackets. www.nature.com/scientificreports www.nature.com/scientificreports/ The sampling sites represent different environments within both New Zealand and Japan (Table 1). New Zealand seawater temperatures are similar between sites and within the optimal range for D. vexillum, with winter minima around 11-12 °C, and maxima between 19 °C and 23 °C. Therefore, we allocated all New Zealand sites a temperature stress score of 1 ( Table 2). While the mean and maximum salinity values are similar between sites in New Zealand (mean approx. 33-35 PSU), the minimum salinity is more variable, with salinity above 30 PSU in Queen Charlotte Sound but dropping below 25 PSU in Pelorus Sound and Whangamata Harbour. The lowest salinity values were reported in Whangamata (13 PSU). Minimum and maximum temperatures and salinities are most likely to influence global DNA methylation differences (see Hawes et al. 46 ). Therefore, each site in New Zealand was allocated a different salinity stress score based on minimum salinity values (1, 2 and 3 respectively; Table 2).

Site of collection
For Japan, minimum water temperatures were similar, or slightly cooler, than minimum water temperatures at sites in New Zealand. Minimum temperatures ranged from 9-12.5 °C at all sites, except for Shizugawa Bay which can be <7 °C in winter. The maximum temperature in Shizugawa Bay was slightly warmer than the maximum water temperatures at New Zealand sites (19-23.5 °C), while Sagami Bay, Ise Bay, and Izu Peninsular all reach maximum water temperatures of 26-28 °C (Table 1). Therefore, all sites in Japan were allocated stress scores of 2 or 3 ( Table 2). For salinity, Shizugawa Bay salinity values were most similar to Queen Charlotte Sound, while Ise Bay salinity values were lower, and more similar to Pelorus Sound and Whangamata Harbour. Major rivers empty into Pelorus Sound (Pelorus River), Whangamata Harbour (Wentworth River) and Ise Bay (Kiso Three Rivers) which can result in low minimum salinity values following precipitation events, so these sites were given stress scores of 2 or 3 while all other sites obtained a salinity stress score of 1 ( Table 2). DNA extraction and methylation-sensitive amplified polymorphism. Tissue samples were finely chopped up using a sterile scalpel blade and genomic DNA was obtained using G-spin Total DNA extraction kits (animal tissue protocol; Intron, Gyeonggi-do, South Korea). Methylation-sensitive amplified polymorphism (MSAP) was then used to determine changes in global methylation patterns as previously described in detail in Hawes et al. 46 .
Data analysis. The Peak Scanner software (Applied Biosystems, Foster City, CA, USA) was used to assign the MSAP fragments peak height and size. To establish the peak height with the lowest error rate (the lowest number of fragment presence/absence mismatches between triplicate samples), and the error rate for msap (v. 1.1.8) analysis 61 described below, we completed the MSAP procedure described above in triplicate for one individual sample. The following settings were associated with the lowest error rate between triplicates: error rate per primer, 0.07; analysis range, 50-500 base pairs (bp); minimum peak height, 50 relative fluorescent units (RFU). Peak presence/absence data corresponding to HpaII and MspI fragments was then converted to a binary matrix (presence = 1, absence = 0), so that the methylation state of each restriction site could be identified. The MSAP profiles were then assessed in R 62 using the msap package (v. 1.1.8) 61 .
The msap package determines whether individual fragments (loci) are methylation-susceptible loci (MSL) or non-methylated loci (NML) by analyzing the contents of the binary matrix. Subsequent analyses are performed independently, with MSL used to assess DNA methylation variation and NML used to assess genetic variation 61 . While NML has been used to assess genetic variation in several studies 51,60,63 , it is only a proxy for genetic variation and should be interpreted with some caution 64 . Methylation state is categorized as either: Type I = unmethylated state, Type II = methylation of internal cytosine, Type III = hemimethylation, and Type IV, which is considered uninformative as this type could be due to either hypermethylation or the absence of the target restriction site 61,65 . Within msap, epigenetic (MSL) and genetic (NML) differentiation between populations were assessed by principal coordinates analyses (PCoA) followed by analyses of molecular variance AMOVA 66 with 10000 permutations 61 . A Mantel test was used to establish the correlation between MSL and NML, an indication of how much epigenetic variation might be influenced by the genetic background 61 . The amount of genetic and epigenetic variation was estimated using the Shannon diversity index (S). Differences between NML and MSL diversity were tested using the Wilcoxon Rank Sum test with continuity correction.  Table 2. Site-specific stress scores related to Didemnum vexillum tolerance ranges based on method described in Ardura et al. 60 . Summative score = temperature score + salinity score.

Site of collection
www.nature.com/scientificreports www.nature.com/scientificreports/ Following analysis with msap, the overall MSL and NML variation estimated by the Shannon diversity index was compared between populations using the number of samples (N), mean diversity (I) and the standard deviation (obtained from msap) to perform and one-way analysis of variance (ANOVA) and post-hoc pairwise comparisons using Tukey's honest significance difference test (Tukey's HSD). The global methylation level was calculated as the proportion of methylated loci (the number of methylated loci (Types II and III) divided by the scorable loci (Types I, II and III)) following 51,67 . Using R (v. 3.3.0) 62 , global methylation levels were compared among populations using a generalised linear model (GLM) with a quasibinomial distribution to correct for overdispersion. Residuals were checked and assumptions validated. We then compared global methylation differentiation in colonies collected from sites with the same summative stress scores using msap (v. 1.1.8) 61 (PCoA followed by AMOVA) to determine whether site-specific MSL differences were associated with differences in environmental conditions (temperature and salinity). Only summative stress scores '3' and '4' could be analysed with msap because scores '2' and '6' were only associated with one site respectively. To assess the relationship between global methylation levels and environmental conditions, we used a GLM with a quasibinomial distribution to compare the proportion of methylated loci in colonies collected from sites with different summative stress scores, followed by post-hoc pairwise comparisons using Tukey's HSD tests.

Results
Environmental differences between sites. DNA methylation patterns in colonies from sites with the same summative stress scores were significantly differentiated (Supplementary Table S1), suggesting that global DNA methylation differences between sites in this study are not clearly associated with summative temperature and salinity scores. There were also no significant differences between the proportion of methylated loci (Types II and III/Types I, II and III) associated with different summative stress scores (Supplementary Table S2).

DNA methylation and genetic differentiation.
In New Zealand colonies (clade A), 1784 loci were identified, 1309 of these were MSL (73.4%, 94% polymorphic) and 475 were NML (26.6%, 97% polymorphic). In Japan clade A, 1780 loci were identified, 1496 of these were MSL (84%, 92% polymorphic), and 284 were NML (16.6%, 97% polymorphic). Finally, in Japan clade B, 1785 loci were identified, 1658 were MSL (92.9%, 93% polymorphic) and 127 were NML (7.1%, 96% polymorphic). Combined, 1799 loci were identified. 1617 were MSL (89%, 99% polymorphic) and 182 were NML (10%, 99% polymorphic). We found no significant differences in MSL or NML patterns between Japan clade A (JAA) and Japan clade B (JAB). Instead, significant differences were due to colonies collected in New Zealand being different to colonies collected in Japan, regardless of clade affiliation (Table 3 and Fig. 2). When clades collected in Japan (JAA and JAB) were compared, the methylation patterns almost fully overlapped (Fig. 2a). However, when colonies collected in New Zealand (NZA) and Japan (JAA and JAB) were compared, they were partially differentiated, regardless of clade (Fig. 2b,c). MSL differences between NZA and JAA were due to colonies collected from Pelorus (PelA) and Whangamata (WhaA) (New Zealand) being different to colonies collected from Japan (Shizugawa (ShiA) and Ise (IseA)). In contrast, there were no  There were no significant MSL or NML differences between colonies collected from different sites within Japan DnA methylation and genetic diversity. MSL diversity was similar in all three populations (JAA, JAB, and NZA), and DNA methylation diversity (MSL) was significantly higher than genetic diversity (NML) (Fig. 3a and Supplementary Table S7). In contrast, there were significant differences in NML diversity between populations (Supplementary Table S8). NML diversity was significantly higher in JAA compared to both JAB and NZA (Supplementary Table S9).

DnA methylation levels.
The proportion of MSL classified as unmethylated, internally methylated, hemimethylated or hypermethylated for JAA, JAB, and NZA is shown in Fig. 3b. In all three groups, a greater proportion of Type IV methylation (hypermethylation) was detected compared to unmethylated or methylated loci. However, as Type IV methylation can also be due to absence of the restriction target (mutation) it is considered uninformative and is not used to calculate the proportion of methylated loci 51,60 . For the proportion of methylated loci, recently introduced New Zealand colonies had significantly lower methylation levels than colonies collected from within the native region (Japan) (p < 0.05, Supplementary Table S10). For JAA, JAB and NZA, the methylated loci (Types II and III) were 49%, 53% and 41% of the scoreable loci respectively (Fig. 3c). There were no significant differences in the proportion of methylated loci between JAA and JAB.

Discussion
Comparing closely related clades, only one of which has become invasive, provides a unique opportunity to identify genotypes and phenotypes that might be associated with invasion success 26 . Significant global DNA methylation (MSL) and genetic (NML) differences were detected between Didemnum vexillum colonies introduced to New Zealand, and both clades in the native range (Japan). These differences were evident when looking at both differentiation and global DNA methylation levels, which will be discussed in detail below. In contrast, we found no significant differences in MSL or NML patterns between D. vexillum clades A and B collected from within the native range, despite the predicted divergence of D. vexillum COI clade lineages being between 12,000 and 2.6 million years ago 18 . A recent study of genome-wide DNA methylation variation within and between five house sparrow subspecies in their native range (Middle East, Iranian plateau) also found the majority of subspecies were not significantly differentiated, despite also inhabiting different environmental conditions (e.g., arid desert, humid forest, tropical coast). Most of the observed DNA methylation variation could be attributed to differences within subspecies rather than among subspecies, and appeared to follow a similar pattern to genetic differentiation among populations, which was also very low 68 . In contrast, a recent study of global DNA methylation in co-occurring subspecies of the facultatively clonal plant, P. australis, found significant differentiation of DNA www.nature.com/scientificreports www.nature.com/scientificreports/ methylation patterns between subspecies that inhabit the same environment 69 . However, these subspecies have very different life histories; one subspecies is native and the other introduced. Therefore, epigenetic differences between these subspecies might be associated with changes induced by invasion history 69 . In D. vexillum, the level of COI sequence variation between clades A and B (as estimated by percentage of polymorphic sites) is low compared to values reported for other ascidian species 18 , and colonies from the two clades are morphologically identical and occur in sympatry at several sites in Japan 37 . The lack of global DNA methylation differentiation between clades A and B in Japan might, therefore, represent similar genetic and life history backgrounds.
Clade A colonies collected from Japan had significantly higher levels of genetic (NML) diversity compared to the two other populations studied (NZA and JAB). While differences in genetic variation do not necessarily reflect differences in quantitative traits, it has been theorized that a high level of standing genetic variation will increase adaptive and evolutionary potential, which may facilitate the colonisation of new environments 70 . High levels of genetic variation within the native range (as found in JAA) may facilitate the spread of invasive species like D. vexillum. In the copepod, Eurytemora affinis, two clades have overlapping distributions in several estuaries and salt marshes, but one clade has invaded freshwater habitats while the other remains restricted to the native  (Type III), and hypermethylation or absence of restriction site (Type IV); and (c) global methylation levels measured as the proportion of methylated loci (Type II + III/scorable loci (Types I + II + III)) following Nicotra et al. 67 and Ardura et al. 51 . Significant differences between treatments are denoted by different letters (p < 0.05).
www.nature.com/scientificreports www.nature.com/scientificreports/ range 25 . The invasive clade was found to have higher DNA nucleotide diversity than the non-invasive clade within the native range 25 , and has evolved increased haemolymph osmolality following the invasion of freshwater habitats 71,72 . Comparisons of geographically restricted and widespread species also suggest that higher levels of genetic variation in ecologically important traits contribute to larger niche breadth 73,74 . While it is possible that future sampling might reveal that clade B has also invaded other regions, the rate of discovery suggests that if it is present, it is rare. Clade A has been identified in five broad sampling areas (New Zealand, Japan, West Coast North America, East Coast North America and Europe), whereas clade B haplotypes have been exclusively amplified from Japan 18 . It is also possible that clade A is only invasive 'by chance' . Clade A and clade B haplotypes are not equally abundant nor equally distributed among sites in Japan 41 , and the probability of a haplotype being transported during an invasion event is proportional to the frequency of that haplotype in the native population where the transport vector originated 75 . This effect would be amplified if introductions to new regions were secondary invasion events rather than multiple primary introductions originating from the native range 41 .
Genetic variation based on NML should be interpreted with caution 64 , but lower genetic diversity in NZA is consistent with previous studies that have concluded that introduced populations of D. vexillum have experienced a reduction in genetic diversity, indicating a founder effect 18 . If lower levels of genetic diversity are in fact limiting adaptive potential, alternative mechanisms for increasing phenotypic variation, such as epigenetic variation, might be especially important during the introduction and establishment phases. Despite differences in genetic diversity, there were no significant differences in DNA methylation diversity between NZA, JAA, and JAB. A similar result was found when DNA methylation diversity was compared across house sparrow populations with varying invasion histories 53 . Combined, these results support the idea that epigenetic variation could preserve the ability of invasive populations to adapt following reductions in genetic diversity that can be associated with colonisation 53,54 , although the functional significance of epigenetic diversity remains unknown.
DNA methylation patterns were significantly differentiated between sites in the introduced range of D. vexillum (New Zealand), but not in the native range (Japan). Recent studies of invasive populations have found DNA methylation differences between environments that may play a role in adaptation to different habitats 49,76,77 . In the present study, comparisons of DNA methylation patterns between sites with similar or different environmental stress (temperature and salinity) scores indicate that DNA methylation patterns are not driven by environmental variations alone. However, it remains possible that environmental stress is an important driver of epigenetic structure in the invasive New Zealand population. Random DNA methylation changes have been shown to be more likely following exposure to environmental stressors 78 . Invasive populations may accumulate novel epigenetic marks more quickly than populations in the native range if the invasion process can be considered a stress 51 , especially if combined with additional environmental challenges 60 . When unpredictable environmental conditions are encountered, epigenetically variable offspring will be favored, serving as a 'bet-hedging' strategy if some individuals are suited to the invaded habitat [79][80][81][82] . Additionally, it is possible that the founding genotypes that successfully established in New Zealand exhibit a higher propensity for methylation plasticity due to environmental filtering. Site-specific differentiation in the introduced populations could be due to the introduction and spread of epigenetic variants to different sites which, like genetic drift, is likely to have a large effect in small populations. However, these processes rely on a level of stability and heritability of methylation. It is not yet clear how much DNA methylation variation underlies selectable phenotypic variation, or how it responds to evolutionary forces 83 .
While we found DNA methylation differences between colonies from New Zealand and Japan, and between sites within New Zealand, these were associated with significant genetic (NML) differences. Invasive populations can include a different admixture of genes or genotypes compared to those in the native range, due to founder effects or multiple introductions from different source populations. Therefore, DNA methylation differences might reflect genetic differences between the native and introduced populations and between sites within the introduced range. Invasive populations can also be exposed to variable selection pressures across the newly colonised range, and have a faster evolutionary response to selection 79 . Models have shown that during early adaptation to new environments, epigenetic shifts could provide a "stepping stone", facilitating the persistence of populations and allowing for the accumulation of neutral genetic variation 84 . Some epigenetic variation might become assimilated, meaning plastic responses or phenotypes are epigenetically encoded first, followed by genetic evolution 85,86 . Stress-induced reductions in DNA methylation can also affect genetic variation by increasing transposition events 87 . The results from this study cannot establish whether DNA methylation differences between sites are pure, obligatory or facilitated, but all three types are likely to play a role in range expansions 23,69 . Selection on epigenetic variation, in addition to genetic variation, can allow populations to adapt faster than selection on genetic variation alone 84 .
The introduced New Zealand population had lower methylation levels compared to populations in the native range. Reduced levels of DNA methylation have been suggested to be a common 'signature of invasion' , increasing phenotypic plasticity during the expansive phase of invasion when it is most required 51,60 . However, the mechanistic link between levels of methylation, environmental stress, biological invasions and adaptation to new environments is not well understood. Few studies have reported methylation levels in aquatic invertebrates and, of those that have, it appears that increased genome-wide methylation is the most common response to environmental stress 47,[88][89][90][91] . Evidence suggests that high levels of gene body methylation can reduce spurious transcription and transcriptional noise (variability of gene expression levels) and may be a defense mechanism to combat the damaging effects of environmental challenges 88,92 . In contrast, genes with reduced levels of germline DNA methylation have been linked with increased phenotypic plasticity via associations with alternative transcription start sites, exon skipping, alternative splicing and transient methylation 52,92,93 . Hypomethylation has also been associated with increased activity of transposable elements in several species, which can introduce phenotypic variation and may facilitate adaptation to new environments 87,94 . While stable and accurate gene expression could be beneficial within the native range of a species, where genetic diversity is generally high and environmental conditions are predictable from one generation to the next, following invasion, noisy or unstable genomes might be especially www.nature.com/scientificreports www.nature.com/scientificreports/ beneficial. If environmental stress is associated with higher levels of DNA methylation, it should be considered that lower levels of methylation might be associated with a reduction in stress. However, there is little evidence for this at present. One study has compared invasive populations of pygmy mussels (Xenostrobus secures) from two sites with apparently varying levels of environmental stress (an international port and a protected lagoon) 51 . It was expected that the population from the less stressful environment (the lagoon) would be less methylated, however it was found that methylation levels between sites were similar. In this case, it was suggested that the 'invasive signature' might be overriding any differences associated with the environment. It should be noted that while Type IV methylation (hypermethylation) is considered uninformative because it could be due to target absence, a large amount could also be due to methylation polymorphisms, so could also be important when comparing methylation levels in invasive and native populations.
conclusion This is the first study investigating DNA methylation differences between closely related invertebrate clades, one of which has become invasive and one of which has not. We found further evidence for a 'signature of invasion' associated with a recent population expansion 51 , with introduced populations having distinct DNA methylation patterns and levels compared to populations in the native range. These results suggest that genetic diversity might facilitate the spread of invasive populations outside the native range, while DNA methylation might be especially important in introduced populations when genetic diversity is reduced. The MSAP method used in this study only allows for analysis of global genomic patterns, so there may be important functional methylation differences between the invasive and non-invasive clades that were not detected here. Our future research will focus on examining the functional properties and stability of DNA methylation changes in response to different environmental conditions, which will aid our understanding of ecological adaptation and evolutionary dynamics, including the epigenetic basis of invasiveness.

Data Availability
The datasets analysed during the current study are available from the corresponding author on reasonable request.