Solidago canadensis s.l. is a perennial Asteraceae native to North America and introduced to Europe in the seventeenth century (Aiton, 1813; Weber, 2000), and to several other continents (van Kleunen et al., 2019). This species reaches heights of over two meters and flowers in Europe from July to October (Schmeil et al., 1993). Multiple shoots re-sprout in the next spring from rhizomes after the aboveground parts have died back in winter (Weber, 2000). Flowers are self-sterile (Melville and Morton, 1982; Schmid and Dolt, 1994), although some degree of autofertility has been found (Razanajatovo and van Kleunen, 2016). This species heavily relies on anemochory with seeds equipped with a feathery pappus (Tackenberg et al., 2003). Solidago canadensis s.l. is represented by several cytotypes in its native range in North America, but so far, only diploid individuals have been found in the invasive range in Central Europe (van Kleunen and Schmid, 2003). In its non-native range, S. canadensis is frequently found in ruderal and disturbed sites (van Kleunen and Schmid, 2003) and, due to its ability to outcompete local flora, has been classified as an invasive species in Central Europe (Freisetzungsverordnung AS 2008 4377, 2008; Bundesamt für Naturschutz, 2013).

Leaf samples were collected from plants grown in a common-garden experiment in the Botanical Garden of the University of Konstanz. The experiment is described in detail in Eckert et al. (2021). In brief, S. canadensis seeds were collected between 2014 and 2015 from 25 wild populations along a latitudinal gradient in Central Europe, with at least five mother plants per population. Information on exact sampling locations, elevation (m a.s.l.), estimated numbers of ramets in the source populations, and the sample sizes of maternal lines are given in Supplementary Table 1. Two batches of seed subsamples were germinated per mother plant, one was treated with an aqueous solution of 25 μM of zebularine during germination and the other was mock-treated with water. Zebularine treatment results in non-specific and dose-dependent genome-wide cytosine demethylation (Zhou et al., 2002), whereby zebularine inhibits DNA methyltransferases leading to hypomethylation during mitosis (Baubec et al., 2009; Champion et al., 2010). Zebularine-derived hypomethylation (hereafter demethylation) has reportedly been shown to be heritable and to modify abiotic and biotic stress responses in different plant species without changing the underlying DNA sequence (Baubec et al., 2009; Verhoeven and van Gurp, 2012).

In total, 400 S. canadensis plants (206 control plants; 194 zebularine−treated plants; 224 maternal lines) were grown in an (incomplete) randomized-block design for 2 years (2015–2016). For each individual, the following phenotypic traits were obtained: initial seedling size, phenology (as days until flowering), plant height, the number of ramets counted in spring of the second year, and total and reproductive-to-total aboveground biomass harvested at the end of the first year. Eckert et al. (2021) found latitudinal clines in all traits (even pronounced in some traits measured on zebularine-treated individuals) except for the number of ramets and the initial seedling size. Latitude of source populations was significantly correlated with climate variables derived from the WorldClim v2.0 database (Fick and Hijmans, 2017). Therefore, we restricted the phenotypic variables in our study to the ones that showed latitudinal variation in Eckert et al. (2021). We refer to Eckert et al. (2021) for further details about the zebularine treatment, the corresponding common-garden experiment and the scoring of plant traits.

For this study, we collected intact leaves from 395 surviving and healthy S. canadensis plants (Eckert et al., 2021) at the onset of flowering in August 2015, and stored them in dry silica gel at room temperature. Genomic DNA was extracted from 10 to 20 mg of leaf tissue using an extraction kit (E.Z.N.A. Plant DNA Kit, omega BIO-TEK Inc., Norcross, Georgia, United States) and about 15% of the samples were randomly chosen to serve as duplicates for later error rate analysis. To assess genetic and epigenetic variation, we performed AFLP and MSAP analysis using protocols detailed in Schulz et al. (2014) and Gáspár et al. (2019) (see also Supplementary Method 1 for changes applied to these protocols). For AFLP analysis, both the restriction and ligation step were combined in one reaction using 500 ng of genomic DNA with the standard enzymes EcoRI and MseI in the AFLP analysis. In each analysis, we used four primer combinations in the selective amplification step (see Supplementary Table 2 for primer sequences) with the selective EcoRI primers labeled with the fluorescent dyes FAM, VIC, NED, and PET. The amplified products were measured on an ABI 3130 Genetic Analyzer (Applied Biosystems, Waltham, Massachusetts, United States) using the GeneScan LIZ 500 sizing standard (Thermo Fisher Scientific, Waltham, Massachusetts, United States).

For MSAP analyses, we followed the same procedures, except that MseI was replaced in two parallel reactions by the isoschizomers HpaII and MspI, which differ in their sensitivity to methylated cytosines (Schulz et al., 2013). The combination of these enzymes with different sensitivity allows the distinction between a maximum of four methylation states for each locus, which is non-methylated if a signal is present in both cases, fully or partially methylated if a signal is present in either of these cases and ambiguous in its status if no signal is present at the specific locus and individual. Selective primer combinations were selected based on pre-testing and quality checking.

After peak sizing of each fragment, each primer combination was analyzed separately and binned using GeneMapper Software v5.1 (Thermo Fisher Scientific, Waltham, Massachusetts, United States). Only peaks with a minimum relative fluorescence unit (RFU) height of 10 within the range of 50–500 base-pairs (bp) were binned as peaks. Binned peaks did not exceed a width of 1.5 bp and peaks were ignored if they overlapped with shoulder peaks or had an odd shape indicating contamination. After peak quality assessment, the final AFLP dataset received an overall error rate of 3.92% on average based on 44 replicates. The final MSAP dataset received an overall error rate of 7.16% on average based on 40 replicates. Details of the genotyping approach are provided in Supplementary Method 2.

To determine the methylation state per peak (hereafter locus), binary tables were analyzed using the MSAP_calc.R script by Schulz et al. (2013) and the scoring option Mix1, which for each locus distinguishes two epigenetic states, methylated (hereafter MSAP-m) and non-methylated loci (hereafter MSAP-n; Schulz et al., 2013). In total, 385 S. canadensis plants (199 control plants and 186 zebularine-treated plants, Supplementary Table 1) were successfully genotyped with an initial number of 856 peaks and, after quality scoring, comprised a final number of 361 polymorphic AFLP loci. Additionally, 336 of these plants (172 control plants and 164 zebularine-treated plants) were successfully epigenotyped with an initial number of 1,396 peaks each and a final number of 187 polymorphic methylated MSAP-m loci and 182 non-methylated MSAP-n loci. Both datasets are openly available in a public repository (Eckert et al., 2022).

To check for population structure along the sampled latitudinal gradient, we used the software STRUCTURE v2.3.4 (Pritchard et al., 2000) on AFLP data of control plants assuming K = 1–11 genetic clusters (Supplementary Figure 1). For each STRUCTURE run, we applied the admixture model and considered the sampling location (LOCPRIOR parameter) with a burn-in of 1,000,000 and 500,000 repetitions, respectively (Porras-Hurtado et al., 2013). Additionally, we tested for isolation-by-distance using the Mantel test with n = 9,999 permutations correlating the geographical log-distance between populations with their corresponding genetic distance as the standardized Wright’s F_ST (Rousset, 1997; Figure 1).

To infer the impact of uneven gene flow in the invasive range, we checked for spatial genetic and epigenetic patterns using Moran’s eigenvector maps (MEM; Dray et al., 2006; Griffith and Peres-Neto, 2006) and the mgQuick-function in the R package MEMGENE v1.0.1 (Galpern et al., 2014; Figure 1). To find the set of eigenvectors that are significantly associated with spatial genetic and epigenetic variation, respectively, we used the geographic coordinates of the source populations and Nei’s unbiased between-population distance (Nei, 1977) calculated with GENALEX v6.5 (Peakall and Smouse, 2012). The mgQuick-function calculates all Moran’s eigenvectors based on the population coordinates and reduces them to a significant set using both forward and backward selection of positive and negative eigenvectors against genetic and epigenetic distance, respectively, via permutation tests. Explained genetic and epigenetic variation of each eigenvector is then assessed using redundancy analysis (RDA). Finally, mgQuick calculates MEMGENE axes, i.e., spatial genetic neighborhoods from eigenvectors based on a principal components analysis (PCA) of the fitted values from RDA representing orthogonal vectors maximizing spatial autocorrelation. We checked significant MEMGENE axes separately for genetic (AFLP) and epigenetic (MSAP-m, MSAP-n) datasets using non-treated control samples. To test the extent to which MEMGENE axes were associated with coordinates of source populations, we calculated a correlation matrix using Pearson’s correlation coefficient. Details of MEMGENE calculations are provided in Supplementary Figures 2, 3.

To investigate whether the zebularine treatment has shifted the presence/absence patterns of AFLP, MSAP-m and MSAP-n loci per sample, we conducted a principal coordinates analysis (PCoA; also known as multi-dimensional scaling) coupled with a multivariate dispersion test (Anderson et al., 2006) based on 9,999 permutations using the betadisper-function from the VEGAN v2.5-6 package (Figure 2). We expected that the zebularine treatment would not affect genetic variation (AFLP and MSAP-n loci) but reduce the number of MSAP-m loci present. To assess between-individual dissimilarity, we used the Sørensen-Dice coefficient (Dice, 1945; Sørensen, 1948) which counts pairwise presence signals as matches. We used the Sørensen-Dice coefficient because we expected that the zebularine treatment would reduce the presence of methylated markers due to its demethylating characteristic. To infer among-population variation against within-population variation, we conducted an analysis of molecular variance (AMOVA; Table 1) according to Excoffier et al. (1992) on control plants based on 9,999 permutations using the poppr.amova-function from the R package POPPR v2.8.6. For both, MSAP-m and MSAP-n, we assessed — at the population level — the total number of loci present, the percentage of polymorphic loci, and the Shannon diversity H_S per population and treatment using the MSAP_CALC v1.3 R script (Schulz et al., 2013; see details in Supplementary Figures 4, 5). For epigenetic diagnostics, we used all corresponding plants per population but also confined this analysis using only sibling pairs (maternal lines) present in both control and zebularine-treatment groups. To check whether the frequencies of epigenetic loci (MSAP-m, MSAP-n) changed due to the zebularine-treatment, we compared control plants vs. zebularine-treated individuals for each locus on a population level. Therefore, we calculated the log₂-fold change (log₂FC) which is common to evaluate the direction of gene expressions (e.g., McCarthy and Smyth, 2009) and visualized log₂FCs per locus and population via heatmaps and box-plots (see details in Supplementary Figure 6). For the calculation of log₂FC, we used only maternal lines that were present in both control and zebularine-treated subsamples.

To check whether zebularine treatment affects the relationship between (epi)genetic and phenotypic population-level variation, we conducted (partial) Mantel tests with n = 9,999 permutations. All (partial) Mantel tests were conducted separately for control and zebularine-treated plants using only plants that were both successfully genotyped and epigenotyped (169 control plants; 159 zebularine-treated plants; Table 2). In partial Mantel tests comparing epigenetic (MSAP-m, MSAP-n) and phenotypic variation, we controlled for genetic variation (AFLP) and spatial genetic neighborhoods from MEMGENE analysis. We applied Nei’s fixation index (Nei, 1977) as (epi)genetic population-level distance and population-level Euclidean distance for both trait variation and spatial genetic autocorrelation data. Trait variation was assessed using the population-average empirical quantile dispersion coefficient (Bonett, 2006) based on plant height and flowering phenology in both years, first-year total biomass and first-year reproductive-to-total biomass ratio, respectively. In addition to the unadjusted p-values from (partial) Mantel tests, we also provide adjusted p-values based on the false discovery rate (Benjamini and Hochberg, 1995; Yekutieli and Benjamini, 1999), but only the unadjusted p-values were used for interpretation.

To screen for (epi)genetic loci of potential adaptive significance (hereafter referred to as outlier loci), we ran genome scans. As suggested by Meirmans (2015), we applied four complementary genome-scan methods that are commonly combined to minimize the risk of false positives: BAYESCAN v2.1 (Foll and Gaggiotti, 2008), BAYESCENV v1.1 (de Villemereuil and Gaggiotti, 2015), redundancy analysis (RDA; Capblancq et al., 2018) and latent factor mixed models (LFMM; Frichot et al., 2013). We ran all methods separately on genetic (AFLP; except RDA) and epigenetic (MSAP-m, MSAP-n) data using non-treated control samples. Detected putative outlier loci per dataset were pooled from each genome-scan approach and visualized as Venn diagrams (Venn, 1881; Figure 3). Conservatively, we only retained outlier loci for further statistical analyses that were detected by at least two approaches.

The applied approaches vary in their algorithm but also in whether they are capable of including environmental variables in their calculations. BAYESCAN uses logistic regression to estimate locus-specific and population-specific contributions to the multilocus F_ST, distinguishing between two selection models, diversifying and balancing selection (Supplementary Figure 7). In contrast to BAYESCAN, both BAYESCENV and LFMM rely on environmental variables as a reference to check for loci of potential adaptive significance (Supplementary Figures 8–13). BAYESCENV implements the same algorithm as BAYESCAN on a modified set of models, a neutral model, a locus-specific model, and the local adaptation model linked to the corresponding environmental variables of interest. LFMM corrects for confounding effects due to underlying population structure by applying admixture coefficient estimation as a first step. We used the same WorldClim 2.0 variables as in Eckert et al. (2021) but accounted for collinearity of these variables via PCA and retained only the first three principal components (PC1, PC2, PC3) for the above-mentioned genome-scan approaches according to the broken-stick criterion (Supplementary Figure 8). To infer outlier loci that might be associated with spatial genetic neighborhoods of populations, we applied RDA on the detected MEMGENE axes as environmental predictors on both epigenetic datasets (MSAP-m, MSAP-n; Supplementary Figure 14) following Forester et al. (2018). Details and settings of the different genome-scan approaches are provided in Supplementary Method 3.

To infer the variation of outlier loci along the latitudinal climatic gradient and across spatial genetic neighborhoods, and whether these were affected by the zebularine treatment, we applied logistic mixed-effects models (LogMM). All models were fitted separately using the glmer-function from the LME4 v1.1-23 package and incorporated the binary (presence/absence) information of each outlier locus as response variable. The sampled source population and maternal lines nested within populations were included as random-intercept factors to account for non-independence of plants from the same maternal lines within populations. For each model, the second-order Akaike information criterion, and the marginal and the conditional R² were calculated using the AICc-function and the r.squaredGLMM-function, respectively, from the MUMIN v1.43.17 package. We provide details of the model parameters in Supplementary Tables 3, 4.

Fixed factors of each fitted model depended on the environmental variables with which the outlier loci were detected (AFLP: loc58, loc286; MSAP-m: loc189, loc282; MSAP-n: loc135, loc176, loc222; see summary in Supplementary Table 5). Models fitted with loci detected with genome-scan approaches that were based solely on PCs from a principal components analysis (PCA) derived from WorldClim 2.0 variables (WorldClim-PCA) included zebularine treatment, corresponding principal components (PCs) and their interaction as fixed factors. We provide details on dimensionality reduction using PCA based on WorldClim 2.0 variables in Supplementary Figure 8. Models fitted with outlier loci detected with genome-scan approaches based on both PCs from WorldClim-PCA and spatial genetic neighborhoods included zebularine treatment, all three MEMGENE axes and their interaction as fixed factors. Because spatial autocorrelation could already point to variation in local climatic conditions (Herrera et al., 2016), PCs from WorldClim-PCA were not included in these latter models. To facilitate model convergence, fixed factors were z-transformed and shifted to positive space by adding a constant, i.e., the highest negative value plus 0.1. For each model, we checked the assumption that the log-odds of the response and each continuous fixed factor was linear applying the Box-Tidwell procedure (Box and Tidwell, 1962). In rare cases where linearity was violated, fixed factors were squared. Models showed no influential values, which was checked for each model separately using Cook’s distance (Cook, 1986), and only moderate levels of multicollinearity among covariates, which was checked with the variance inflation factor (Dormann et al., 2013).

Fixed effects were tested for significance using likelihood-ratio tests (LRT) by comparing the full model to a model without interaction and the non-interaction model with models where single terms have been iteratively removed (Lewis et al., 2011). In addition to the unadjusted p-values from LRTs, we also provide adjusted p-values using the false discovery rate (Benjamini and Hochberg, 1995; Yekutieli and Benjamini, 1999), but interpretation was based only on the unadjusted p-values. Predicted marginal-effect values were visualized with the GGPLOT2 v3.3.0 package and all results were obtained using R v4.0.2 (R Core Team, 2020).

We found weak but significant spatial genetic structuring in the analyzed S. canadensis populations along the sampled latitudinal gradient. Based on coordinates of S. canadensis source populations, three MEMGENE axes (hereafter MEMGENE1, MEMGENE2, MEMGENE3) significantly explained 13.2% of the genetic variation underlying Moran’s eigenvector maps (MEMs; n = 1,000 permutations; Figure 1C). On each MEMGENE axis, shared genetic neighborhoods among populations are given as similar population-level values and their corresponding signs (negative or positive; Supplementary Figure 2). MEMGENE1 explained 48% of variation underlying MEMs and divided the analyzed populations into a southern (populations 1–13) and a northern genetic subgroup (populations 14–25). This corresponds to results of admixture coefficients estimation via LEA within the genome-scan framework (Supplementary Figure 10 and Supplementary Method 3) as well as to the STRUCTURE analysis where K = 2 was detected as the most probable number of genetic clusters according to Evanno et al. (2005) (see population-averaged cluster probabilities from STRUCTURE analysis in Figure 1A). MEMGENE1 significantly correlated to latitude of the sampled source populations (R² = −0.94; p < 0.001), whereas MEMGENE2 and MEMGENE3 axes were independent of the latitudinal gradient (Supplementary Figure 3).

MEMGENE2 (29% of variation explained underlying MEMs) and MEMGENE3 (23% of variation explained underlying MEMs) split populations into multiple alternating subgroups where subgroups at the top of the sampled latitudinal gradient as well as in the southern part (MEMGENE2) or only in the southern part (MEMGENE3) split the remaining subgroup (Figure 1C). AMOVA showed that there was significant genetic differentiation among populations explaining around 8% of genetic variation (Table 1), but this was not reflected in IBD along the sampled latitudinal gradient (Figure 1B). Coordinates of source populations did not significantly explain any epigenetic variation in the MEMGENE analysis and no pronounced epigenetic structuring was found with LEA (Supplementary Figures 11, 12). In contrast, AMOVA showed significant epigenetic differentiation between populations explaining 2.7% (MSAP-m) and 4.3% (MSAP-n), respectively, of epigenetic variation (Table 1).

We found that more than one-third (37.9%; 66 loci) of all MSAP-m loci showed a fourfold decrease (log₂FC ≤ −2) in occurence across zebularine-treated individuals compared to control plants, but also a fourfold increase (log₂FC ≥ 2) was present in 22.9% loci (40 loci) pointing to a bias for decreased methylation from the zebularine treatment. In addition, half of the MSAP-n loci (49.4%; 77 loci) showed a fourfold decrease in occurrence, whereas only 18.5% (29 loci) showed a fourfold increase. On average, up to 20.2 epigenetic loci per population at least halved in frequency whereas up to 15.9 epigenetic loci at least doubled in frequency (Supplementary Figure 6). As expected, zebularine treatment did not affect genetic between-individual distance as revealed by PCoA and betadisper analyses (Figure 2A). Surprisingly, the overall between-individual distance for epigenetic loci was also not significantly reduced by the applied demethylation treatment (Figures 2B,C). Similarly, population-level descriptors of epigenetic diversity, i.e., the number of epigenetic loci present, the percentage of polymorphic epigenetic loci and the Shannon diversity index, were not significantly affected by the zebularine treatment (Supplementary Figures 4, 5). Additionally, epigenetic population-level variation correlated with genetic variation as revealed by Mantel tests (Table 2). This relationship, however, vanished for MSAP-m but not for MSAP-n loci when looking at zebularine-treated plants of the same populations (Table 2). Population-level genetic and epigenetic diversity did not significantly reflect phenotypic population-level variation (Table 2). This was also the case for epigenetic population-level variation when controlling for genetic diversity and spatial genetic neighborhoods (Table 2).

Our genome-scan approaches only revealed few outlier loci in S. canadensis populations that showed signs of adaptive significance. BAYESCAN detected only a single genetic outlier locus (label: loc185; size: 135 bp; marker sequence: ACG-CAT) but no epigenetic outlier loci (Figure 3A). BAYESCENV and LFMM, which relied on WorldClim 2.0 variables as environmental predictors derived from PCA (Supplementary Figure 8), detected up to 15 outlier loci in AFLP (BAYESCENV: 8; LFMM: 9), MSAP-m (BAYESCENV: 4; LFMM: 15), and MSAP-n datasets (BAYESCENV: 3; LFMM: 9). From this pool of putative outlier loci, both approaches jointly detected six genetic and epigenetic outlier loci (AFLP: n = 2; MSAP-m: n = 1; MSAP-n: n = 0; Figures 3A–C). One of these jointly detected genetic loci (label: loc58; size: 58 bp; marker sequence: AAC-CCT) was significantly less probable in zebularine-treated plants as revealed by LogMM, although with an increasing probability along PC1 in both control and zebularine-treated plants (Figure 4A and Table 3). Any other jointly detected outlier locus was not significantly affected by the zebularine treatment and did not significantly vary with climatic variation among source populations (Table 3 and Supplementary Table 3).

The RDA approach, which was applied only on epigenetic datasets (MSAP-m, MSAP-n) and relied on spatial genetic neighborhoods (MEMGENE axes; see Supplementary Method 3) as environmental predictors, detected eight epigenetic outlier loci (Figures 3B,C). From this pool, a subset of four putative outlier loci, that were jointly detected by additional approaches to reduce false positives, was further analyzed (Table 4). One of these epigenetic outlier loci (label: loc135; size: 135 bp; marker sequence: AGG-CGA) was significantly less probable with increasing values of the spatial genetic neighborhoods (MEMGENE1), however, this was evident for control plants only but not for zebularine-treated individuals (Figure 4B). All other jointly detected putative epigenetic outlier loci were not significantly affected by either the zebularine treatment or spatial genetic neighborhoods of source populations (Table 4 and Supplementary Table 4).

In contrast to clonal, apomictic or selfing species, outcrossing invasive species, such as S. canadensis, should be more affected in their genetic structuring at different spatial scales due to unequal gene flow and dispersal barriers (Ward, 2006). Spatial spread has been found to be often non-uniform, non-linear or to proceed from several locations simultaneously (Hastings et al., 2004). Based on Moran’s eigenvector maps, we found that genetic variation, but not epigenetic variation reflects the spatial distribution of S. canadensis plants. Epigenetic variation has been linked to rapid phenotypic evolution (Zhang et al., 2013), which might be beneficial for outcrossing invasive species, e.g., in ruderal sites which are characterized by spatial heterogeneity and unstable environmental conditions creating potential for epigenetic memory and priming (Turgut-Kara et al., 2020). Inheriting the parental epigenome, however, was also suggested to being less beneficial after long-distance dispersal into a new range and advantageous only in environments where gene expression states were already fine-tuned to cope with local conditions (Verhoeven and Preite, 2014). Linking epigenetic diversity of the outcrossing species S. canadensis to spatial variation in the invasive range in Central Europe, we found that long-term spread within the invasive range did not result in spatial epigenetic structuring. The results of our study indicate that cytosine methylation appears to be less significant in adaptation processes of outcrossing invasive plant species that exhibit high levels of admixture. Our findings support the idea that, in contrast to clonal or apomictic non-native plant species that might exhibit a high degree of epigenetic differentiation (Verhoeven and Preite, 2014), in outcrossing and admixed invasive plant species epigenetic differentiation might rarely exceed genetic variation and thus will likely not result in population structuring.

Spatial heterogeneity has been found to create and to increase time lags during invasion and to act as an environmental filter during the different invasion stages (Theoharides and Dukes, 2007). Rather than environmental selection on genotypic composition, the spatial genetic neighborhoods might also reflect that the analyzed S. canadensis populations vary in their stage of introduction or could be affected by bottleneck as well as by multiple introduction events at different spatial scales (Ward, 2006). Herrera et al. (2016) found that epigenetic between-individual similarity in Helleborus foetidus, a native plant species for which clonal propagation is exceptional, was greater than genetic similarity at the shortest distances. They pointed out that differences in local environmental features might play a similar or even greater role than spatial distance for epigenetic population structure, arguing for isolation-by-environment rather than IBD (Herrera et al., 2016). Our study focused on climatic variation and spatial autocorrelation along a latitudinal gradient, therefore, we did not include additional local site-specific features into our analyses. For example, both local soil properties and functional diversity have been found to affect the performance of S. canadensis in its invasive range in Central Europe (Czortek et al., 2020). Therefore, future studies on S. canadensis in Central Europe might more strongly emphasize local and site-specific conditions and their impact on epigenetic variation.

Ancestry estimation using the LEA framework indicated two epigenetic clusters, however, not structured along the latitudinal gradient under study. This is in line with Lele et al. (2018), who analyzed genetic and epigenetic differentiation in the non-clonal plant species Vitex negundo var. heterophylla in its native range and found that adaptation to heterogeneous habitat conditions was mainly genetically driven. The results of this study, however, also indicated a weak connection of epigenetic diversity and adaptive phenotypes. In our study, we did not find a significant relationship between epigenetic population-level diversity and trait variation, but there was low but significant epigenetic population-level differentiation as revealed by AMOVA. This epigenetic population-level differentiation may be linked to adaptation involving traits that were not measured in this study and that could be correlated to local environmental features. Our results, however, speak against an adaptive role of spatial epigenetics in S. canadensis in the invasive range. The three significant spatial genetic autocorrelation axes found in our study, however, might indicate that adaptive dynamics in S. canadensis, if present, take place at different spatial scales simultaneously and in a non-linear fashion.

Besides evident spatial genetic neighborhoods, S. canadensis populations were divided along the latitudinal gradient by two genetic clusters as revealed by the STRUCTURE and LEA approaches although with no IBD dynamics. This finding is surprising given that Eckert et al. (2021) found an overall lack of genetic population structure along this gradient based on seven microsatellite markers. In contrast to microsatellites, AFLPs are dominant markers with no information about heterozygosity, but have been found to be more efficient in individual-based population-assignment when there is weak population structuring (Campbell et al., 2003). Yet, population differentiation in S. canadensis based on AFLP analysis was weak and only a small amount of genetic variation could be explained by spatial patterns. While our results support the microsatellite analysis in Eckert et al. (2021), weak population differentiation is an unexpected result given the large latitudinal gradient and the long residence time of Solidago canadensis in Central Europe. Future approaches, however, might focus on more sophisticated molecular techniques such as SNPs (Manel et al., 2010) or spatial modeling based on large numbers of loci to map ecological adaptation (Fitzpatrick and Keller, 2015). The importance of including spatial heterogeneity into modeling plant invasion has been demonstrated in several studies (Herrera et al., 2016; Štajerová et al., 2017). Based on their results for Helleborus foetidus, Herrera et al. (2016) even argued that genetic and epigenetic differentiation might develop independently of each other leading to contrasting spatial patterns from spatially divergent selection. More studies are needed to assess whether this is the case for S. canadensis.

To the best of our knowledge, our study is the first to analyze heritable epigenetic spatial structuring in non-clonal offspring of an invasive plant species by quantifying effects of experimental cytosine demethylation. Although the demethylation agent zebularine shows dose-dependent downsides, such as delayed height and root growth (Cho et al., 2011; Eckert et al., 2021), it has been useful to detect stress-related transgenerational epigenetic variation in clonal and apomictic plant species (Richards et al., 2012; Preite et al., 2015). In our study, we found one epigenetic locus that was affected by the zebularine treatment, i.e., the presence of this locus varied depending on spatial variation. Variation in epigenetic patterns using the AFLP/MSAP technique has already been demonstrated in different plant species in their natural habitat. For example, Schulz et al. (2014) have detected one epigenetic locus in Viola elatior whose probability varied with the amount of photosynthetic active radiation present. Additionally, changes in methylation related to temperature have been shown for different plant species, e.g., Arabidopsis thaliana (Finnegan et al., 1998) and Festuca rubra (Münzbergová et al., 2019). The underlying plant material for our study, however, was only sampled once, at the end of the first year of the common-garden experiment. Therefore, we cannot exclude that there were changes in methylation profiles of control plants resulting from the common-garden environment itself, e.g., epigenetic shifts due to microclimatic conditions in the common-garden environment. In general, the majority of the analyzed 369 epigenetic loci (MSAP-m: 187 loci; MSAP-n: 182 loci) remained nearly unchanged, when zebularine effects were compared within populations. When there were changes in frequency, they were bidirectional: Up to 44.1% of epigenetic loci, i.e., MSAP-m and MSAP-n loci combined, showed a fourfold decrease, but up to 21.3% showed a fourfold increase in frequency. This non-directional change might be attributed to an increase in the level of transposable elements due to zebularine treatment, as found for Arabidopsis thaliana ecotypes (Konečná et al., 2021). With our approach, we could link the variation in two of these anonymous loci to climatic and spatial genetic variation, but more specific associations will need further analyses of the underlying mechanistic context, e.g., using next-generation approaches (Schield et al., 2016; Werner et al., 2020).

Although more epigenetic loci decreased than increased in mean frequency in zebularine-treated population subsets, the applied zebularine treatment did not result in a global change in the number of loci per population, percentage of polymorphic loci and Shannon diversity. Unfortunately, we did not sample at the beginning of the growing season, so the long-term stability of the applied zebularine treatment and, thus, persistent detectable changes remains to be tested. While we saw visible changes in plant phenotypes during the germination phase, the concentration applied during the germination phase might have been too low to cause persistent detectable changes. Hence, possible treatment differences, which may have been present at the beginning, may no longer have been visible in methylation profiles of later stages. In fact, plant DNA has been shown to be capable of compensating for experimental demethylation through repair mechanisms (Liu et al., 2015). If this is the case here, then this implies that experimental cytosine demethylation in S. canadensis can be compensated for within one growing season and, that cytosine methylation appears to be short-term only in this species.

Genetic and epigenetic population-level diversity in S. canadensis were strongly correlated. This is not surprising given that epigenetic variation is frequently found to be genetically controlled (Richards et al., 2017). In our study, however, we found that zebularine treatment decoupled genetic from epigenetic (MSAP-m) diversity. Still, there are plant species in which epigenetic variation may also exceed genetic variation or more strongly relate to a distinct environmental variable, e.g., grazing intensity in Plantago lanceolata populations in Germany (Gáspár et al., 2019). Although epigenetic variation might play a role in plastic responses to environmental challenges in the field in S. canadensis, the results of our study suggest that cytosine methylation most likely does not play a crucial role in transgenerational adaptive changes in its invasive range in Central Europe.

Solidago canadensis is often found in anthropogenically disturbed habitats (Weber, 2000). These habitats are often hotspots of continued introduction of non-natives and it seems that introduction rates for S. canadensis have not declined since this species has been introduced in Central Europe in the seventeenth century (Weber, 2000). Thus, variation in (epi)genotypes, i.e., variation in the presence of (epi)genetic loci, could also be due to multiple introductions and secondary spread pointing to unrelated genotypes. A similar scenario, has been assumed for the invasive hexaploid S. canadensis populations in China (Guo et al., 2016). We do not know the donor populations of introduced S. canadensis, so direct comparisons with populations in their native North American range were not possible. Future studies should do intensive sampling in the native range and analyze genome-wide nucleotide variation to infer possible source populations (e.g., Vallejo-Marín et al., 2021). It would also be interesting to expand the latitudinal gradient in the invaded area in future studies to include aspects of land use and urbanization related to genetic and epigenetic effects at different scales (e.g., Gáspár et al., 2019).

Multiple studies have demonstrated elevational and latitudinal phenotypic clines in fitness-related traits for S. canadensis in the invasive range in Central Europe, and some of these clines persisted when offspring were grown in a common-garden environment (Weber and Schmid, 1998; Moran et al., 2017; Eckert et al., 2021). These phenotypic clines should thus be reflected at the genetic level, but possibly likewise at the epigenetic level. Applying three genome-scan approaches and subsequent logistic mixed-effects models, we were able to detect one genetic and one epigenetic locus putatively under selection. One of these genetic loci was significantly associated with climatic variation along the latitudinal gradient and one epigenetic locus was significantly correlated to spatial genetic variation. A comparably low number of two genetic outlier loci with signatures of selection were jointly detected for Viola elatior (Schulz et al., 2014) when applying multiple genome-scan approaches, though compared to a much higher number of 39 detected epigenetic outliers. Conversely, Herrera and Bazaga (2010) found that only 10 out of 23 candidate loci were significantly associated with epigenetic differentiation in Viola cazorlensis pointing to a link between genetic and epigenetic divergence. The low number of only two candidate outlier loci detected in our study might be attributed to ongoing introductions that maintain admixture in S. canadensis and, thus, introduce unrelated (epi)genotypes (Verhoeven et al., 2011; Dlugosch et al., 2015). In this case, the outlier loci would be false positives despite their detection by multiple genome-scan approaches. It could also be a sign of a yet cryptic differentiation, which could be obscured by already pre-existing (epi)genetic variation from the invasion process of this species. Therefore, it remains to be tested whether this might be due to adaptive processes or genetic drift. Likewise, putative non-neutral epigenetic changes might have triggered the Baldwin effect, i.e., plastic responses to environmental cues induced via epigenetic variation might eventually be replaced by genetic change (Simpson, 1953; Bräutigam and Cronk, 2018). This process is yet hypothetical and still needs to be tested in plants (Bräutigam and Cronk, 2018).

Our study of S. canadensis in Europe revealed the presence of one genetic locus associated with climatic variation and one epigenetic locus associated with spatial variation along the studied latitudinal gradient. Interestingly, not all spatial patterns detected with Moran’s eigenvector maps and introduced into the genome screening using RDA resulted in the same, if any, putative loci with signatures of selection. This could suggest that within the analyzed latitudinal gradient, small-scale heterogeneous dynamics and non-linear spatial corridors (Hansen and Clevenger, 2005) between populations might affect allele frequency more strongly than isolation-by-distance processes along the gradient under study. For example, De Kort et al. (2020) used Fragaria vesca offspring collected as seeds from both a steep altitudinal (<2 km) and a wide spatial (>500 km) gradient and have grown individuals in a controlled environment to characterize DNA methylation signatures at different spatial scales. Their study demonstrated that epigenetic differentiation arises at varying spatial scales with CG methylation divergence more pronounced at the fine-scale altitudinal gradient. Since plant material for our analyses was collected at the end of the first growing season in the common-garden, we can only map one point in the lifetime of the plants. Many studies on epigenetic structuring among plant populations collected the samples directly in the field (see e.g., Herrera and Bazaga, 2010; Bewick and Schmitz, 2015; Gáspár et al., 2019; De Kort et al., 2020). Under field conditions, epigenetic structure, if present, could be directly affected by the local environment of the species. In our case, because we sampled in a common garden, the results would be due to transgenerational epigenetic inheritance. Still, it remains to be tested whether variation in the detected loci will be transgenerationally inherited and which underlying mechanistic processes these loci might trigger.

No genetic and epigenetic outlier loci were detected jointly by all three (for epigenetic loci four) genome-scan approaches. This could be because each of these approaches followed different strategies: BAYESCAN is able to detect stabilizing vs. diversifying selection if present (Foll and Gaggiotti, 2008), and BAYESCENV (de Villemereuil and Gaggiotti, 2015) and LFMM (Frichot et al., 2013) were developed to detect loci associated with environmental variables. Although the AFLP/MSAP approach is still useful for non-model organisms in detecting population structuring, it does not provide the regulatory context of the genomic sequence, i.e., linking gene-expression states to complex trait variation. Therefore, future approaches might overcome this limitation by focusing on more advanced methods of genome-wide association, such as QTL analysis (Long et al., 2011), bisulfite sequencing (Lu et al., 2020) or epiGBS (Schield et al., 2016; Werner et al., 2020).