Conserving on the edge: genetic variation and structure in northern populations of the endangered plant Dracocephalum ruyschiana L. (Lamiaceae)

Loss of biodiversity is accelerating, including the loss of genetic diversity. Conservation of small, isolated populations may be important, as they can provide valuable contributions to overall genetic variation and long-term viability of species. Furthermore, such populations may play an essential role in adaptation to new environments following changes in e.g. land-use and climate. Dracocephalum ruyschiana is a threatened plant species throughout its European distribution, but 25% of the European populations are situated within Norway. Therefore, the species has its own action plan in Norway, which includes demographic monitoring. However, this monitoring does not cover genetic variation nor is the selection of monitored populations based on genetic differentiation, therefore this fundamental level of biodiversity is overlooked. We analyzed 43 sites using 96 SNPs developed for D. ruyschiana, to investigate whether the monitored populations cover the genetic variation and differentiation found within the Norwegian distribution. The results show structuring and differentiation between populations and indicate that there are at least four distinct genetic groups, of which only two are covered extensively by current demographic monitoring. We suggest that two sites representing the two other genetic groups should be included in the national monitoring program to better conserve the genetic variation found in the Norwegian population of D. ruyschiana. Overall, our results highlight the importance of an integrated, interdisciplinary framework to better monitor and conserve biodiversity at several levels.


Introduction
The loss of biodiversity is accelerating with negative consequences for populations, species, communities and ecosystems (Hughes et al. 2008;Ceballos et al. 2015). The main threat to biodiversity loss is land-use change (IPBES 2019), which causes population fragmentation that may lead to loss of genetic diversity (Aguilar et al. 2008). Genetic variation is the smallest unit of biodiversity, and thus crucial for conserving biodiversity at higher levels, such as species and ecosystems (Bruford et al. 2017). Conserving genetic diversity provides genotypes for selection in a changing environment (Barrett and Schluter 2008), contributes to faster recovery of populations after disturbance (Reusch et al. 2005), and increases the stability of ecosystem processes (Crutsinger et al. 2006;Des Roches et al. 2018).
An important international political document for halting biodiversity loss is the United Nations' Convention on Biological Diversity (CBD; www.cbd.int), which also gives an explicit goal to protect genetic diversity. The fundamental idea is that conservation should aim to conserve evolutionary processes, adaptive potential, and not just current species (Bowen 1999). However, very limited action has yet been taken on a global scale to fulfil this goal (Laikre 2010; Electronic supplementary material The online version of this article (https ://doi.org/10.1007/s1059 2-020-01281 -7) contains supplementary material, which is available to authorized users. abilities to adapt to changing environmental conditions (Keller and Waller 2002;Barrett and Schluter 2008). Moreover, small populations are more prone to genetic drift and inbreeding than larger populations, further lowering their genetic diversity (for example, Wagner et al. 2012). Regardless of levels of genetic variation, populations at the distribution margin may play an important role in the range shift of species following climate change (Razgour et al. 2013). Furthermore, such populations can also hold unique ecological and genetic variation (Osborne et al. 2012) and potentially enhance within-species genetic diversity (García-Ramos and Kirkpatrick 1997). Thus, conserving these populations may be especially important, particularly if the species is threatened across its entire distribution.
To conserve genetic diversity and guide practical conservation strategies, it is crucial to identify population-genetic structuring, such as genetic distinctiveness and diversity, and ideally to couple this information with knowledge on geographical barriers and divergent local ecological adaptation. In Norway, only a few rare and endangered plant species have been investigated genetically with the intention to guide conservation plans (Westergaard et al. 2011a, b;Birkeland et al. 2017;Westergaard et al. 2019).
Here we use the vascular plant 'northern dragonhead' Dracocephalum ruyschiana L. as a case to investigate how analyses of genetic variation can be used to guide decision makers and management when conserving biodiversity nationally. The populations of D. ruyschiana are rapidly declining in Europe: the species is listed under Appendix 1 of the Bern Convention (http://conve ntion s.coe.int/Treat y/ FR/Treat ies/Html/104-2.htm) as well as on many national and regional Red Lists (e.g. Gärdenfors 2005;Király 2007;Witkowski et al. 2003), including the Norwegian Red List, in which it is listed as vulnerable (VU) (Henriksen and Hilmo 2015). Habitat loss due to urbanization or decreasing management of agricultural land, is considered the main threat to the species in Norway (Solstad et al. 2015), as well as in several European countries (Käsermann and Moser 1999;Lazarevic et al. 2009). Despite the documented decline of D. ruyschiana, knowledge of genetic variability, structure and evolutionary history of this species is still lacking for Norwegian and global populations.
Dracocephalum ruyschiana (Lamiaceae, 2n = 2x = 14) is a long-lived, insect-pollinated herb with a continental, fragmented distribution in Eurasia (Lid and Lid 2005;Lazarevic et al. 2009). It is a typical steppe plant, and grows on exposed, well-drained, shallow calcareous soils in dry meadows and rocky outcrops (Lid and Lid 2005;Lazarevic et al. 2009). In Europe (excluding Russia), more than 25% of the remaining populations are found in Norway (Norwegian Directorate for Nature Management 2010). The species reaches its European northwestern limit in Norway, where it has a southeastern distribution through the country (Faegri andDanielsen 1996, Lazarevic et al. 2009 (Lid and Lid 2005). The flowering period is typically in the middle of June, and the flowers are blue, 2-2.5 cm long, and gathered in a spike-like inflorescence with 50-60 flowers (Fig. 1d). The flowers are mainly pollinated by long-tongued insects such as bumblebees Bombus ssp. (Milberg and Bertilsson 1997). The fruits are typical for the family Lamiaceae, resulting in four dry one-seeded nutlets, without obvious adaptations for long-distance dispersal. The species is most likely self-compatible (Milberg and Bertilsson 1997), as self-incompatibility systems are not known from the Lamiaceae (Owens and Ubera-Jiménez 1992;Allen and Hiscock 2008).
In Norway, D. ruyschiana is one of only three vascular plants regulated to have priority for conservation (Lovdata 2011). As a large part of the European population is located in Norway, it is crucial to conserve the species nationally. Thus, an action plan for conservation has been prepared (Norwegian Directorate for Nature Management 2010) and a monitoring program for Norwegian populations of D. ruyschiana developed (Evju et al. 2016), aiming to collect data on population size and structure to inform managers on status and trends for the species at the national level. The monitoring protocol has so far been tested on 18 populations in and around the city of Oslo (Evju pers. comm.). The selection of populations was not based on genetic variation or structuring in the species, but on geographic clustering (regions). Furthermore, the action plan includes no actions to conserve genetic diversity.
Here we use single nucleotide polymorphisms (SNPs) developed for D. ruyschiana (Kleven et al. 2019) to investigate genetic diversity and structure within Norway. In particular, we aim at identifying genetic groups within its Norwegian distribution to evaluate whether established national monitoring sites sufficiently cover the observed c Norwegian distribution of D. ruyschiana (blue dots, data from Norwegian Biodiversity Information Centre, artsdatabanken.no), the collected sites for this study (yellow dots) and sites that are currently included in demographic monitoring (red dots). The site ID indicated in the map follows Table 1

Specimen sampling, DNA extraction and SNP dataset
With collection permits from local authorities (see Acknowledgements), leaf material of D. ruyschiana was sampled from 43 sites across the Norwegian distribution range (Table 1, Fig. 1), with relatively sparser coverage of the distribution margins. Samples were taken in 2012-2014, prior to the development and establishment of the monitoring program in 2016 and 2017, respectively, therefore there is not complete overlap between our sites and the monitored sites. Two different sampling schemes were used when collecting leaves. In the first scheme, ten 1 m 2 -plots were randomly placed across a site, after first identifying all flowering clusters of D. ruyschiana. Leaves were then sampled from the plant closest to the center of the plot. In total, 20 sites were sampled using this scheme, and later several of these sites were included in the monitoring initiated in 2017. For the second sampling scheme in the remaining 23 sites, plants were sampled randomly. Where possible, 10 individuals were sampled from each site, however, in sites with fewer than 10 clusters of individuals (likely belonging to the same clone), leaves from one plant representing each cluster were collected.
Genomic DNA was extracted from the silica dried leaves using either the NucleoSpin Plant II extraction kit (Macherey-Nagel) or DNeasyPlant Mini Kit (Qiagen), following the manufacturer's protocols. The amount of extracted DNA was quantified on a Qubit 2.0 using the HS Assay kit (Thermo Fisher Scientific). All samples were genotyped at 96 presumably non-coding SNPs applying a SNP-typing assay recently developed specifically for the purpose of facilitating monitoring of the northern dragonhead (see Kleven et al. 2019 for details). Briefly, SNPs were genotyped on a 96.96 Dynamic Array using the Fluidigm EP1 instrument according to the manufacturer's protocol and scored using the Fluidigm SNP genotyping analysis software (https :// www.fluid igm.com/softw are). Deviation from Hardy-Weinberg equilibrium for each SNP in each population was estimated using Arlequin ver. 3.5.1.2 (Excoffier and Lischer 2010). As none of the 96 SNPs deviated significantly from Hardy-Weinberg equilibrium consistently across populations, we included all SNPs in the downstream analyses.

Genetic diversity and structure
Genetic diversities in each D. ruyschiana site were explored by calculating mean number of alleles (N A ), observed (H O ) and expected (H E ) heterozygosity, the inbreeding coefficient (F) and percentage of heterozygous loci in GenAlEx 6.5 software Smouse 2006, 2012). Furthermore, structuring of genetic variation within and between sites, and between the six geographic regions was estimated in analyses of molecular variances (AMOVA), and population differentiation at the regional level was estimated using the F ST -value. We performed one more analyses in GenAlEx 6.5: the assignment test by using the "leave one out" option to explore the origin of individuals based on their genotypes. We also tested whether genetic distance was correlated with geographical distance, using the Mantel test (mantel.randtest) with 999 permutations with the R package ADE4 v1.7-13 (Dray and Dufour 2007) Genetic structuring was explored using MULTISPATI-PCA (Dray et al. 2008) performed with the R package ADE4 v1.7-13 (Dray and Dufour 2007) and SPDEP v1.1-2 (Bivand et al. 2008). This MULTISPATI-PCA is an extension of the principal component analysis that also takes the spatial covariability among variables into consideration. To further explore the number of genetically homogeneous groups (K) and overall structuring in the dataset, we ran genetic cluster algorithms in STRU CTU RE v2.3.4 (Pritchard et al. 2000;Falush et al. 2003Falush et al. , 2007Hubisz et al. 2009). In the initial analysis, the admixture and correlated allele frequencies models were applied with 50,000 runs as burn-in followed by 150,000 runs. The analysis was run 10 times for each K from 1 to 15. The probability of obtaining the genotype data X given K (Mean LnP(K)) increased with the number of Ks, however, the increase was only marginal from K = 10 (not shown). Using the ΔK statistics of Evanno et al. (2005), the most likely number of groups from the preliminary run was 2. However, using the default setting of α may lead to inaccurate estimates of K and assignments probabilities when population sizes are uneven and under the assumption that populations are descendants of recent ancestral populations (Wang 2017). Thus, we reran the analyses using 50,000 runs as burn-in followed by 150,000 runs, but with adjusted α-values based on the initial analyses (α = 1/assumed optimal K; α = 1/10 and α = 1/2), applying the admixture and correlated frequencies models.
The number of likely genetic groups was explored using calculations of Mean LnP(K) (Pritchard et al. 2000) and ΔK (Evanno et al. 2005) in Structure Harvester (Earl and vonHoldt 2012). To compare Structure results at multiple values of K, we aligned and visualized bar plots using the Clumpak (Cluster Markov Packager across K) web server (Kopelman et al. 2015).

Results
A dataset of 355 individuals collected from 43 sites was genotyped for 96 SNPs. The number of individuals sampled from each site varied between two and 13 (Table 2). Eight sites sampled five or fewer individuals, due to small total population sizes, 11 sites had between six and nine individuals, 23 sites had ten individuals, while one had 13. The level of missing SNP data was 0.003%. Genetic variation is summarized in Table 2 Table 2). The AMOVA results (Table 3) showed that most of the genetic variation was found within sites (87%), while 8% and 5% of variation was found among sites and geographical regions, respectively. The total F ST value was 0.134 (p < 0.001), while all values estimated between regions were lower than 0.10 ( Table 4). The lowest F ST value was found between Oslofjorden and Randsfjorden (F ST = 0.027) and the highest between Valdres and Hedmark (F ST = 0.081). The Mantel test revealed a positive correlation between genetic distance and geographical distance (R = 0.56, p = 0.001, Fig. 2).
The estimated likely number of K was not consistent across analyses, however, the different runs yielded the same results across sites, with only small variations in resolution among genetic clusters. Figure 3 shows the major mode of runs in Structure of K = 2-4 (10/10 runs) and K = 11 (9/10 runs) using the admixture and correlated frequencies models and α = 0.1. The Mean LnP(K) increased as the number of K increased and had not reached an optimum at K = 15 (Online Resource 1). The ΔK estimates indicated that K = 2 was the most likely number of genetic groups for most analyses (Online Resource 1). A combined interpretation of the results obtained for increasing K values showed a hierarchical resolution of genetically homogeneous groups (cf. Meirmans 2015). Individuals from sites in Oslofjorden and site H43 Solberg clustered together, while the rest of the individuals belonged to another genetic group at K = 2 (Fig. 3). In the latter group, more individuals were admixed between the two groups than were individuals in the former group. When K = 3, site T17 Vik splits off in a separate group (Fig. 3). There were no admixed individuals in this site and it stayed homogeneous as K increased. The nearby sites consisted of individuals admixed with the "purple" Vik-cluster and mainly the "blue" genetic cluster, but most of the sites split off to separate genetic groups as K increased (see K = 11 in Fig. 3). At K = 4 (Fig. 3), sites R24 Aslaksrud, H41 Bergseng and H43 Solberg separated into a fourth homogeneous group. These sites split into separate groups as K increased (see Fig. 3, K = 11). Other sites also split off in distinct groups as K increased, and there was a slight tendency that the sites in the Oslofjord region separated in two; an eastern group and a western group (see K = 11 in Fig. 3).
The Structure results were supported by the MULTIS-PATI-PCA result, showing that individuals from sites in the Oslo region group together, while individuals from sites in Randsfjorden group together (Fig. 4). Furthermore, individuals from Tyrifjorden group together, with site T17 Vik at a distance from the other sites. Site H43 Solberg was intermediate of sites from Oslofjorden and Randsfjorden regions, and forms its own group. This is in accordance with the Structure results for K = 4 (Fig. 3). The genetic variation observed among sites from Oslofjorden and Randsfjorden is thus covered in the monitoring program; these two regions represent the "blue" and "orange" genetic clusters identified by Structure (Fig. 3).
The admixture of genetic groups within sites recognized by Structure was also evident from the results of the assignment test. The test showed that most individuals originated within their sampled region, but within the regions there seemed to be a higher exchange of individuals between sites (Online Resource 2). In sites with no or low numbers of admixed individuals, such as T17 Vik, R24 Aslaksrud, and H43 Solberg, no individuals were assigned to other sites.

Discussion
Our results indicate the existence of at least four distinct genetic groups of Dracocephalum ruyschiana in Norway, of which only two groups are covered by the current monitoring program. To better conserve the genetic variation present in Norway, and thus a substantial portion of the European genetic variation, we suggest that the monitoring program should include sites with genetic groups not yet covered and downscale coverage of other sites.

Genetic diversity and structuring
Genetic variation in Norwegian populations of D. ruyschiana was estimated by expected heterozygosity (total H E = 0.27) and the levels found were relatively similar among populations (ranging from 0.20 to 0.32). Compared with other threatened plants (e.g. Tero et al. 2005;Dostálek et al. 2010; but see Minasiewicz et al. 2018), the genetic variation in the studied D. ruyschiana sites seem to be moderate. As different genetic markers have been applied in these studies, the interpretation of genetic variation should be treated cautiously. There was no sign of inbreeding, as the inbreeding coefficient was close to 0 for most sites. Thus, outcrossing seems to be maintained by pollinating insects (Milberg and Bertilsson 1997). Sites with few individuals have high negative inbreeding coefficients, but are not considered as this is likely an artefact of low sample size. Most of the genetic variation was found within sites (87%), indicating low differentiation between sites (F ST = 0.134). This is comparable to what was found in the close relative D. austriacum in the Czech and Slovak Republics; high genetic diversity within populations (80%) and relatively low differentiation among populations both within and between regions (Dostálek et al. 2010). In D. ruyschiana, isolation by distance was observed which may bias the Structure results (see Perez et al. 2018 and references therein). However, the spatial multivariate analysis (MULTISPATI-PCA) takes spatial data into account and the result obtained is very similar to the Structure result. Therefore, we use the two analyses combined to explore the genetic structuring of Norwegian D. ruyschiana. The PCAresult show that sites in the Oslofjorden region, Randsfjorden region and Tyrifjorden region, form three separate clusters. This finding is supported by the Structure analyses (K = 3) showing that sites in the Oslofjorden region belong to the same genetic group ('orange' group, Fig. 3), including site H43 Solberg (Hedmark region) and only a few individuals from other sites. All other sites form a separate group at K = 2 ('blue' group, Fig. 3), but at K = 3, a 'purple' group is revealed in Tyrifjorden region. As the number of K increases, more sites splits off in new genetic groups. At K = 4, H43 Solberg splits off in a 'green' group and also forms its own cluster in the PCA. Hence, genetic variation is geographically structured within the Norwegian distribution of D. ruyschiana.
The Structure results indicate that there is admixture between genetic groups, as two or more genetic groups are presented in single individuals (Fig. 3). Furthermore, individuals from two genetic groups are often found within one site. The assignment test also estimated that several individuals originated from other sites. However, between regions there seems to be limited exchange of individuals. This pattern may reflect gene flow within regions, which is likely as this is an outbreeding species dependent on insects such as Bombus species, for pollination (Milberg and Bertilsson (Dostálek et al. 2010). This could explain the observed isolation by distance. However, F STvalues between regions are very low and might reflect a historically larger meta-population. Land-use change and degradation have lately led to a more fragmented distribution of D. ruyschiana habitat (Evju et al. 2014), with an assessed 30-50% reduction in population size (Henriksen and Hilmo 2015). If the current distribution range of D. ruyschiana becomes even more fragmentated, a population decline can be expected, potentially leading to a decline in genetic variation and increased genetic isolation (Rubidge et al. 2012). A combination of isolated populations and decreasing numbers of pollinators (Hallmann et al. 2017) could further restrict future outcrossing within and between populations of D. ruyschiana. Reduced outcrossing rates may affect fitness, such as seed production (Milberg and Bertilsson 1997;Castro et al. 2015).
Review studies show a general positive relationship between population size, genetic variation and fitness, meaning larger populations tend to have higher levels of genetic variation and higher fitness (Leimu et al. 2006). Indeed, higher seed production was found in larger, and more  Fig. 3 Genetic structuring in Dracocephalum ruyschiana was investigated using Structure 2.3.4. The barplots show how 355 individuals from 43 sites using 96 SNPs were assigned to different genetic groups (K), when K = 2, 3, 4, and 11 (indicated to the left of the barplots), using admixture and correlated frequencies models and α = 0.1. For K = 2-4, the barplots summarize 10/10 runs, while for K = 10 the major mode 9/10 runs are shown. The site ID follows Table 1 Fig. 4 The PCA shows a spatial multivariate analysis of Dracocephalum ruyschiana. The boxes each represent one site and the site ID follow Table 1 and Fig. 3. Short arrows reveal a local spatial similarity, while a long arrows reveal a spatial discrepancy. The colours indicate the genetic group identified by Structure at K = 4 genetically diverse populations of D. austriacum, suggesting that genetic diversity might directly affect plant fitness and populations within this genus (Dostálek et al. 2010, but see Plenk et al. 2019). This finding could also be expected for the close relative D. ruyschiana, though further studies will be necessary to confirm this. Mimura et al. (2017) highlight the importance of conserving genetic variation and suggest monitoring of intraspecific genetic variation at global, regional, and local scales. Dracocephalum ruyschiana is threatened both at regional scales (for example Henriksen and Hilmo 2015) and continental scales (Europe; Bern Convention list I). This is the first genetic study of the species, and the first aimed at guiding monitoring and conservation plans. To prioritize among populations in conservation genetics, the delimitation of management units (MUs), defined as populations with different allele frequencies, has been suggested as a suitable conservation approach (Moritz 1994). Based on an interpretation of genetic structure and diversity in Norwegian populations and regions of D. ruyschiana, we presently recognize at least four management units corresponding to the groups recognized at K = 4 ("orange", "blue", "purple", and "green") by Structure (Fig. 3). The Structure result showed that as the number of K increases, more sites seem to split off in separate genetic groups. Nevertheless, the four groups cover all regions and genetic extremes (Fig. 4) studied so far. However, spatial genetic structuring may also be present even at a local scale (Minasiewicz et al. 2018;Tero et al. 2005) and regions such as the valleys stretching north-westwards, represented by few sites and/or individuals, should be further explored for undetected genetic variation in D. ruyschiana.

Implications for conservation in Norway
Current demographic monitoring of D. ruyschiana is concentrated on sites in the regions Oslofjorden and Randsfjorden, which together cover the two main genetic groups recognized ("orange" and "blue", respectively). The southernmost sites in Tyrifjorden seem to be the most genetically distinct sites in this region, with site T17 Vik being the most distinctive (Fig. 4); however, none of these sites are covered by ongoing monitoring. Hence, as a minimum, site T17 Vik should be included in the current monitoring program. Sites from the Hedmark region are also distinct in the analyses, and inclusion of site H43 Solberg in the monitoring program would increase the coverage of genetic variation monitored in Norway. At the same time, monitoring is unduly intense in some areas and some of the sites could be excluded based on genetic similarities, especially from Randsfjorden (Fig. 4). Criteria for exclusion should be based on representativity of habitat types and population states.
The SNP markers used were developed to guide monitoring of genetic variation within Norway and are based on individuals collected throughout the Norwegian distribution (Kleven et al. 2019). If genetic monitoring of the species is implemented in Norway, the markers are well suited for the purpose, as they can easily be applied to additional sites and/or used to monitor the species through time as the results are reproduceable. However, applying other genomic tools could likely aid in identifying locally adaptive genetic variation, revealing other management units (Funk et al. 2012). Genetic analyses of D. austriacum using AFLP markers, demonstrated that neutral and adaptive genetic variation was not correlated (Bonin et al. 2007). Thus, considering both neutral and adaptive genetic variation may be necessary to best conserve the evolutionary potential of a species (Razgour et al. 2018;Mable 2019). Furthermore, increasing our understanding of genetic structure and adaptive variation in D. ruyschiana on both regional and global scales, requires the application of modern genomic tools as well as phylogeographic analyses on a global dataset.
In 2009, Norway updated its laws on management of biodiversity (The Biodiversity Act (Naturmangfoldloven 2009)). The new act aims to give the authorities the possibility to halt the loss of biodiversity, including genetic diversity, through protection and sustainable use of natural resources. Dracocephalum ruyschiana is one species protected by this law. However, studies show a positive correlation between human population density and the loss of rare species (Thompson and Jones 1999), and as D. ruyschiana's core distribution overlaps with the most densely populated areas of Norway, the species is still under pressure of habitat loss and fragmentation (Evju et al. 2014). Loss of populations may have negative consequences for future communities and ecosystems (Des Roches et al. 2018). Therefore, active measures must be applied to prevent further fragmentation and loss of populations of D. ruyschiana. To better protect the species from becoming locally extinct and conserve its evolutionary potential throughout its distribution range, ecologists, geneticists, and conservation managers should work together to develop an integrated, interdisciplinary framework to better inform monitoring programs and conservation actions (Flanagan et al. 2017;Razgour et al. 2018).
proofreading the article. We are grateful to two anonymous referees for valuable comments on a previous version of this manuscript.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.