Landscape genetic analysis suggests stronger effects of past than current landscape structure on genetic patterns of Primula veris

Recent changes in land use have led to substantial loss and fragmentation of semi‐natural grasslands. We assessed the relative effect of current and historical landscape composition, and landscape change on the genetic diversity and gene flow of a characteristic grassland plant.


| INTRODUC TI ON
Recent changes in land use have led to drastic loss and fragmentation of natural and semi-natural habitats (Picó & Van Groenendael, 2007).
In particular, land use has substantially changed towards more intensive agricultural practices, while traditional management has been abandoned (Picó & Van Groenendael, 2007;Prentice et al., 2006).
During the last century, the area of semi-natural grasslands has drastically decreased in Europe (Cousins et al., 2015;Hooftman & Bullock, 2012) due to lack of traditional hay-making and extensive grazing. Decrease in the area and connectivity of habitats, and altered habitat conditions often lead to loss of biodiversity (Hahn et al., 2013;, including the decrease in genetic diversity of wild plant populations (Leimu et al., 2010). Genetic diversity provides adaptive potential to species for coping with environmental change and should therefore be one of the focal aims of conservation management (Aavik & Helm, 2018). Furthermore, inclusion of knowledge about genetic diversity in restoration activities has a high potential to increase restoration success (Mijangos et al., 2015) and is crucial for recovering structurally and functionally resilient ecosystems (Moreno-Mateos et al., 2020).
Fragmentation and loss of habitats often impose serious negative effects on genetic diversity because they potentially lead to reduced gene flow between populations and increased inbreeding and genetic drift Leigh et al., 2019;Leimu et al., 2006;Prentice et al., 2006;Young et al., 1996). Increasing application of landscape genetic tools has advanced knowledge about the influence of environmental heterogeneity on genetic diversity and demonstrated that, in addition to geographic isolation, also the characteristics of the landscape around and between plant populations can strongly shape genetic patterns of plants (Emel et al., 2020;Holderegger et al., 2010;Lehmair et al., 2020;Manel et al., 2003).
This knowledge is also crucial for the recovery of resilient populations during ecosystem restoration (Moreno-Mateos et al., 2020).
As plants mostly depend on various abiotic and biotic vectors for dispersal of pollen and seed (Holderegger et al., 2010), landscape change affects plants mainly through influencing the movement of these vectors, either positively or negatively (Holderegger et al., 2010). For insect-pollinated grassland species, for example, woody elements may act as a barrier for gene flow by inhibiting the movement of pollinating insects (Aavik et al., 2014;Aavik et al., 2017;Aguilar et al., 2006;Hahn et al., 2013;Schmitt et al., 2000). Corridors in barrier elements can facilitate pollen dispersal (Tewksbury et al., 2002), which is particularly critical for maintaining gene flow in insect-pollinated plant species (Aavik et al., 2014(Aavik et al., , 2017. On the other hand, edge habitats, for example on the borders of woody landscape elements (Sõber et al., 2020) or between farmlands (Hass et al., 2018), can benefit pollinating insects by offering a more suitable microclimate (Bergman et al., 1996) and potential nesting and feeding areas (Neumüller et al., 2020) facilitating pollinator movement and thus gene flow across grasslands between edge habitats. Further, grazing cattle may facilitate the distribution of propagules of grassland plants, even of those not specifically adapted to animal seed dispersal (Holderegger et al., 2010). Therefore, grazing, especially rotational grazing, can be an important mechanism for maintaining genetic diversity via increased gene flow (DiLeo et al., 2017;Honnay et al., 2006;Jacquemyn et al., 2010;Plue et al., 2019; as well as species diversity (Pärtel et al., 2007). However, nowadays, due to changes in land use, cattle is rarely being moved across larger spatial scales, which may have also contributed to reduced seed dispersal.
Plant populations, which have experienced severe fragmentation, may still maintain considerably high genetic diversity despite landscape changes (Reisch et al., 2017). This phenomenon may be caused by so-called genetic extinction debt, due to which genetic diversity has not yet reacted to the changes of landscape structure but might be expected to do so in the future (Aguilar et al., 2008).
Indeed, several studies have found evidence of such time lags in the genetic patterns of plants Münzbergová et al., 2013;Reisch et al., 2017). Delayed responses to landscape change are more probable in long-lived plant species, while lagged responses in the genetic patterns of short-lived plants are rare (Aavik et al., 2017;Helm et al., 2009). Yet, there is still a gap in knowledge about the effects of change in landscape characteristics, such as the loss and transformation of semi-natural grasslands, on the gene flow between increasingly isolated grassland plant populations. A time lag in the response of genetic diversity and gene flow to habitat fragmentation may substantially influence conservation recommendations, but unfortunately, potentially delayed responses are often not accounted for despite the threat of further decay in genetic diversity of the species and eventually even extinction (Essl et al., 2015).
However, delayed responses could also be seen as a possibility in conservation to restore the historical habitats of species before the negative effects of habitat loss manifest.
Whether the response of genetic patterns of plants to current and historic landscape structure is detected may also depend on the applied method as well as the metrics of genetic diversity (Epps & Keyghobadi, 2015). In node-based analyses, the effect of landscape features (i.e. amount of different habitats) in a circular area at and around a study population on population-based genetic diversity is tested, using, for example, linear model or regression analyses (DiLeo & Wagner, 2016). With such on-site analyses, effects on genetic connectivity via potential immigration and establishment of pollen, seeds or juvenile plants can be tested. Link-based approaches, in contrast, focus on examining the effect of landscape features in a direct or indirect path between pairs of populations on K E Y W O R D S ddRAD, edge density, genetic diversity, grassland, landscape change, singlenucleotide polymorphisms their genetic differentiation (F ST ), using Mantel test or (multivariate) generalized linear model approaches with distance matrices (DiLeo & Wagner, 2016). Such between-site analyses provide the possibility to test for facilitating or inhibiting effects of the intervening landscape on gene flow and its vectors (e.g. pollinators) via pollen and seeds. The two approaches complement each other and should ideally be used in concert (DiLeo & Wagner, 2016). Yet, this has been rarely done so far (DiLeo & Wagner, 2016).
In the current study, we examined the genetic diversity of Primula veris populations occurring in semi-natural calcareous grasslands in Western Estonia that have experienced a dramatic loss of their area during the last century (Helm et al., 2006). In addition, we assessed genetic differentiation between the study populations as an indirect measure of gene flow. We examined the response of these node-and link-based genetic metrics to current and historical (pre-fragmentation; from the 1930s) landscape configuration and to landscape change. We proposed the following questions: (a) Does the amount of grassland in the landscape affect the genetic diversity of populations and gene flow between populations? (b) Does the amount of woody elements (forests and shrubs) in the landscape influence the genetic diversity of populations and gene flow between populations? (c) What is the role of grassland edge density on the genetic diversity of P. veris? (d) Do historical landscape characteristics have a stronger effect on current genetic diversity than current landscape features due to relatively long life span of P. veris, that is do patterns of genetic diversity exhibit a lagged response to landscape change? The knowledge obtained in the present study is crucial for interpreting measures of biodiversity in spatio-temporally dynamic landscapes and for planning conservation activities of threatened habitats.

| Study species
Primula veris L. (Primulaceae) is an herbaceous perennial plant with an approximate life span of up to 50 years (Ehrlén & Lehtilä, 2002).
Primula veris flowers in May in Estonia and is mostly found in dry, species-rich grasslands, in shrub or woodland ridges and edges and on calcareous cliffs. It is a light-demanding and drought-tolerant species (Brys & Jacquemyn, 2009). Primula veris is an obligate outcrossing plant depending on insects for effective pollen flow. The flowers are mostly pollinated by different species of Hymenoptera (e.g. bees), but also some species of Coleoptera (beetles) and Lepidoptera (butterflies). Pollen dispersal of P. veris is generally limited to a few metres from parental plants and seed dispersal is restricted to a few centimetres from maternal plants (Brys & Jacquemyn, 2009) but occasional longer dispersal distances might be reached by domestic and wild animals (Auffret & Plue, 2014;Plue et al., 2019;. The study species is heterostylous, that is it has two morphologically different flower types. Mostly, only inter-morph pollen exchange results in successful fertilization (Ganders, 1979).

| Study area
Study sites were located in semi-natural calcareous grasslandsalvars-in Western Estonia on the islands of Muhu and Saaremaa (58.6°, 23.25°; 58.42°, 22.5°; Figure 1a). Annual mean temperature in the area is 6.9℃, and precipitation is 600 mm (Estonian Weather Service, 2020). Alvars have characteristic shallow calcareous soils occurring on limestone bedrock. Estonian alvars are typically located near coastal areas in Northern and Western Estonia, making the climate of these grasslands more humid and with a smaller temperature range. Alvars have lost more than two thirds of their area during the last century because traditional grazing needed to sustain the grasslands has been mostly abandoned. As a consequence, most of the alvars have overgrown with trees and shrubs (Helm et al., 2006;Pärtel et al., 1999).  Table 1). The distance between sampled individuals within a population was at least one metre. Samples were stored on silica gel until DNA extraction.

| Genetic analyses
Samples were weighed to 25 mg of leaf material and ground for 2 min using two 2.3-mm metal beads in a Mixer Mill 301 (Retsch GmbH). DNA was extracted using the LGC sbeadex plant maxi kit (LGC) with some modifications. Extraction steps of binding, washing and elution were done on a KingFisher Flex Purification System (Thermo Fisher Scientific). For a more detailed description of DNA extraction, see Träger et al. (2021).
Extracted DNA was prepared for sequencing using double-digest restriction site-associated DNA sequencing (ddRADseq; Peterson et al., 2012). An extended description of sequencing library preparation is provided in Träger et al. (2021). Briefly, the ddRADseq method uses two restriction enzymes (in the present study, EcoRI and TaqI) to cut standardized DNA in a two-step process. Following a purification step, DNA fragments were ligated to corresponding adapters (48 EcoRI adapters and 2 TaqI adapters). Samples containing different EcoRI adapters, but the same TaqI adapter (48 samples), were Reads (sequenced DNA fragments) were bioinformatically analysed and filtered (see Appendix S1, for more info). Genotype information was extracted from the resulting VCF file using PGDSPIDER v2.1.1.3 (Lischer & Excoffier, 2012). Putatively adaptive SNPs were excluded using environmental association analysis (EAA) to detect SNPs associated to environmental factors related to habitat overgrowth and fragmentation and outlier tests to detect SNPs under potential diversifying or balancing selection (see Träger et al. (2021), for more information) to ensure that only putatively neutral SNPs were used in the further analyses. Population-based genetic diversity indices (unbiased expected and observed heterozygosity, uH e and H o , respectively, and percentage of polymorphic loci, %P) were calculated using GenAlex version 6.503 (Peakall & Smouse, 2005. Unbiased expected heterozygosity accounts for differences in population sizes and is shortened as expected heterozygosity (H e ) below. Inbreeding coefficients (F IS ) and genetic differentiation (F ST ) were calculated using the package "genepop" (Rousset, 2008) in R version 3.4.2 (R Core Team, 2017).

| Landscape data
For assessing the role of landscape characteristics on the genetic patterns of P. veris, we used historical and current data about the distribution of calcareous grasslands and woody landscape elements because overgrowth of grasslands with trees and shrubs has been the most profound change in these landscapes. We obtained historical landscape data (woody elements and grasslands) from digitalized maps of historical vegetation survey, which was carried out in the 1930s, that is when the distribution of Estonian alvar grasslands was at its maximum (Laasimer, 1965). Current grasslands were as- F I G U R E 1 Map of the study area and study design. a-Location of 19 study populations of Primula veris on the Estonian islands of Muhu and Saaremaa and the distribution of current and historical alvar grasslands in the region. b-An example of the node-based method to calculate the amount of current grassland area in circular buffers with radii of 500, 1,000 and 2000 m around the study populations of Primula veris. c-An example of the link-based method to calculate the amount of current grassland area around a straight line between two study populations of P. veris with radii of 200, 500 and 1,000 m

| Data analysis
Because predictor variables from different buffer radii used for the node-based approach were correlated with each other (Appendix S2), we made full models for each buffer radius separately (Appendix S3), using only uncorrelated variables in one model. We excluded woody elements change (r = 500 m, 1,000 m, 2,000 m), water amount (r = 1,000 m, 2,000 m) and historical grassland amount (r = 2,000 m) in models due to many strong correlations with other variables. TA B L E 1 Geographic coordinates, the number of sampled individuals before and after bioinformatic filtering, expected heterozygosity (H e ), observed heterozygosity (H o ), inbreeding coefficient (F IS ) and the percentage of polymorphic loci (%P) in 19 study populations of Primula veris on the islands of Muhu and Saaremaa, Estonia All statistical analyses were done in R version 3.4.2 (R Core Team, 2017). We conducted linear models to examine the effect of landscape variables on expected and observed heterozygosity and inbreeding coefficient of P. veris. We carried out generalized linear models to examine the response of the percentage of polymorphic loci applying a quasibinomial error distribution (package "lmerTest" (Kuznetsova et al., 2017)). We also included region (Muhu/Saaremaa) as a co-variable in the analysis to take into account potential landscape historical differences between regions (Appendix S4). In addition, we included the interactions between region and other landscape variables in the models. We also added population size in the models, but because it had no significant effect in any of the models, we present results without population size. We made several full models (Appendix S3) because of the correlations between landscape variables (Appendix S2  reached a plateau at about 27 km (Träger et al. 2021) indicating that at larger distances other (random) factors than geographic distance and/or landscape characteristics might be more important drivers of genetic differentiation. We chose the best models according to deviance information criterion (DIC). We used model selection based on DIC to compare different models to a null model, which included only geographic distance as an independent variable and population identity and region of both populations in a pair as random variables.

| Landscape analysis
Landscape variables calculated within different buffers (r = 500,  Note: Some landscape variables were correlated, so several initial models were made (Appendix S3). Genetic diversity indices (Gen div ind) used were expected heterozygosity (H Different measures of genetic diversity were on average higher in Muhu than in Saaremaa ( Table 2). The response of genetic diversity to landscape variables in the node-based approach depended on the scale (r = 500, 1,000 and 2,000 m) as well as the measure of genetic diversity. In the best models with the lowest AICc, historical edge density exhibited a most constant and significant positive effect on different variables of genetic diversity. In particular, both H o and H e increased with historical edge density within a 500-metre buffer radius around populations (Table 2; Figure 2). Furthermore, H o showed a significant positive response to historical edge density also in the models with landscape variables calculated within the buffers of 1,000 and 2,000 metres (Table 2) (Figure 3b). As the DIC of the best model differed less than 2 from four other models, a single best model may not be determined. However, the significant variables in these models still stayed the same, so we report the results of only the best model.

| The role of landscape composition on genetic diversity and gene flow
We found that higher (historical) proportion of alvar grasslands between the study populations has a negative effect on the genetic differentiation (F ST ) of P. veris populations, indicating more gene flow between populations with higher proportion of historical grasslands between them. Thus, high proportion of suitable habitats between grassland plant populations supports intra-specific gene flow by ensuring functional connectivity, that is the movement of seed and pollen (Auffret et al., 2017). In the study area, pollinator observations on "LIFE to alvars" (Helm, 2019) project sites revealed that the overgrowth of alvars with shrubs and trees was accompanied by a significant loss in the diversity and abundance of important pollinator groupsbutterflies and bumblebees (Prangel, 2017). It is thus very likely that Note: Variables in the model: F ST -genetic differentiation, dist-geographic distance, curr. grasscurrent grassland proportion, hist. grass-historical grassland proportion, curr. woody-current woody elements proportion, hist. woody-historical woody elements (forests and shrubs) proportion. All variables were calculated as proportional amounts in the buffers and rank transformed. DIC values for models containing only geographic distance as explanatory variable are independent of the buffer zone and have thus the same value. Model DIC values with ΔDIC < 2 within the buffer radius are highlighted in bold. The best model for each radius according to DIC is underlined.
*The overall best model according to DIC.

TA B L E 3
Deviance information criterion (DIC) values of models from the link-based analysis for three buffer zones (r = 200, 500 and 1,000 m) for the direct path between population pairs of Primula veris on the islands of Muhu and Saaremaa changes in landscape structure, such as the overgrowth of alvars and consequent loss of suitable habitats, can lead to deficiency of pollen vectors, which, in turn, is likely to have an overall negative effect on genetic diversity and gene flow. However, in semi-natural grasslands such as alvars, achieving functional connectivity may be also ensured by the movement of cattle or sheep as plant seed vectors between grassland complexes. With the recovery of cattle-grazing, semi-natural habitats might be preserved by preventing overgrowth of grasslands with shrubs on the one hand and by spreading plant seeds within and among grassland patches on the other (Rico, Holderegger, et al., 2014).
At the study sites, sheep and cattle were introduced (after sample collection) as the main management tools for ensuring the maintenance of "LIFE to alvars" project sites.
We found that the proportion of polymorphic loci (%P) and observed heterozygosity (H o ) of P. veris populations were generally lower in landscapes, which had experienced more drastic grassland loss. In addition, %P decreased with decreasing habitat amount (r = 500 m), whereas H o , surprisingly, increased with decreasing habitat amount (2,000 m). The area of grasslands has been found to have a positive effect on the genetic diversity in grassland plants (e.g. Aavik et al., 2019;Hahn et al., 2013;Toma et al., 2015). This is partly mediated by the positive effect of habitat amount on population size and thus the potential gene pool (DiLeo & Wagner, 2016;Leimu et al., 2006;Young et al., 1996). However, in our study, population sizes of P. veris might be still sufficiently large to not affect genetic diversity because we did not find a significant relationship between these two factors. The longevity of the study species should also be considered, because it might additionally buffer potential habitat and population size effects on the measures of genetic diversity (see section "Delayed response of genetic patterns to landscape change" below). Yet, in the long term, it can be expected that loss of habitat amount may lead to decreased population size that in turn could increase the negative influence of genetic drift (Balkenhol et al., 2015), which causes the random loss of alleles, as well as inbreeding (Young et al., 1996).  (Vellend & Geber, 2005) compared to homogeneous conditions. However, as the present study relies on neutral genetic markers, which should not respond to selection, this might not be the most applicable explanation. Second, edges of grasslands may be important habitats for pollinators, whose contribution to the gene exchange in insectpollinated plants may be substantial (e.g. Jabis et al., 2011). Indeed, pollinators have been shown to prefer edge habitats due to more feeding and nesting opportunities (Diaz-Forero et al., 2011Sõber et al., 2020) as well as more suitable microclimate (Bergman et al., 1996) and possibly to shrubs due to being potential landmarks (Cranmer et al., 2012). Additionally, in our study, current and historical grassland edge densities were strongly negatively correlated with current and historical amount of grasslands, respectively, meaning sites with more grasslands have less grassland edges per area unit. Thus, the negative influence of the amount of grasslands on genetic diversity observed in some of the models could likely mirror the preference of pollinators towards landscapes with more edge habitats. Therefore, alvar grasslands should still harbour some degree of shrub coverage (Helm, 2019), while the connectivity between sites should be maintained (Bergman et al., 2018;Steffan-Dewenter & Tscharntke, 1999) and overgrowth avoided.
The results revealed a positive effect of historical (but not current) amount of woody elements in the proximity of studied alvar grasslands on the percentage of polymorphic loci (%P) of P. veris. The effect of forest on the genetic diversity of insect-pollinated grassland plants is usually assumed to be negative as it can act as a barrier for pollinators and thus diminish gene flow (Aavik et al., 2014(Aavik et al., , 2017DiLeo et al., 2018). However, Hahn et al. (2013) did not detect the negative influence of forest on the genetic diversity and gene flow of Trifolium montanum. It has been hypothesized that forest would reduce the movement of pollinators and thereby hamper gene flow in insect-pollinated grassland plants when these woody landscape elements cover more than half of the landscape ). Yet, several pollinators are able to cross forests and shrubs up to certain distances (e.g. Zurbuchen et al., 2010) and some forest cover may be even beneficial to, for example, butterflies (Bergman et al., 2018). In addition, the relatively low historical amount of forest and shrubs may have acted similarly to grassland edges for pollinators providing feeding and nesting sites as well as better microclimate.

| Delayed response of genetic patterns to landscape change
We found support for patterns of genetic diversity of P. veris to exhibit a lagged response to landscape change, although some of the genetic variables responded to current landscape parameters.
Notably, historical grassland edge density had a significant effect Working in the same region as the study area, Aavik et al. (2019) found some evidence of lagged responses (at larger scale, r = 2,000 m) of the genetic diversity of T. montanum to landscape variables in alvar grasslands, but Helm et al. (2009)

| CON CLUS IONS
The present study is one of the few to use simultaneously both node-and link-based approaches for analysing the effects of landscape change on genetic patterns of plants. We demonstrate that despite the vast landscape change in Estonian alvars, genetic diversity as well as patterns of gene flow of Primula veris still reflects past landscape context. Therefore, relationships between current habitat amount and genetic patterns should be interpreted with caution.
Clarifying the impact of landscape change on the genetic diversity and differentiation of plants is thus crucial for guiding decisionmaking in nature conservation policy. Where possible, maps about historical distribution of habitats may help to determine delays in biodiversity patterns. Such delays may provide opportunities for successful restoration of semi-natural grasslands before they suffer from the consequences of substantially decreased genetic diversity.
Most of the effects of landscape variables on the measures of within-population genetic diversity of P. veris were apparent at the scales of 500 and 1,000 m, while the role of landscape was less evident at the scale of 2,000 m. Because these distances correspond F I G U R E 3 Relationship between pairwise genetic differentiation (F ST ) of Primula veris study populations in the corridor analysis and geographic distance (a) and historical proportion of grasslands (b). Eff. indicates the effect sizes of the best generalized mixed effect model (Table 4)

CO N FLI C T O F I NTE R E S T
The authors declare no conflicts of interest.

PEER R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/ddi.13357.

DATA AVA I L A B I L I T Y S TAT E M E N T
Sequence data used in this study will be made available at the