Introduction

Species’ distribution ranges are commonly defined by ecological limits, which are most often determined by ecological gradients. As species approach their range limits, their populations typically become smaller and more fragmented (Bridle and Vines 2007). As a consequence, range edge (margins or verge) populations often feature decreased genetic diversity due to reduced gene flow as result of the more fragmented distribution or random genetic drift effects enhance genetic diversity due to neutral, negative, or positive mutations which can appear and increase in frequency over time (e.g., Ellstrand and Elam 1993). If gene flow patterns align with geographic distance (isolation by distance, IBD; Wright 1943), this affects the distribution of variation within and among range edge populations (Eckert et al. 2008; Sexton et al. 2009). If ecological gradients cause range limits, populations are generally more likely to maintain specialized genotypes that are well adapted to particular ecological conditions (e.g., Rehm et al. 2015). This is especially the case when populations that are geographically closer together are genetically more similar than populations that are further apart (IBD) and where genetic and environmental differences among populations (isolation by environment, IBE) are positively correlated. Especially in smaller and more fragmented range edge populations, positive, advantageous mutations are more likely to become fixed, as directional selection within and among range edge populations might be stronger than in central populations, where stabilizing selection through higher rates of gene flow tend to oppose effects of local selection and therefore limits adaptation (Hoffmann and Blows 1994; Lenormand 2002; Bridle and Vines 2007; Sexton et al. 2009). Under an IBD and IBE pattern, range edge populations eventually diversify and undergo niche evolution during adaptation to novel environments, or alternatively may depauperate where adaptation is prevented by small population size (Sexton et al. 2009).

For decades, evolutionary ecologists have investigated local adaptation across different systems and scales, however, rarely on micro-geographic scales. This is due to the assumption that high rates of gene flow prevent adaptive divergence at fine spatial grains (Richardson et al. 2014). However, several investigations on small geographic scales suggest that micro-geographic divergence is more widespread than commonly assumed (Kettlewell 1955; Antonovics and Bradshaw 1970; Steiner and Berrang 1990; Kavanagh et al. 2010; Willi and Hoffmann 2012; Krueger-Hadfield et al. 2013; Richardson and Urban 2013; Richardson et al. 2014).

Environmental association studies allow the investigation of the nature of local adaptation by identifying its leading causes (Savolainen et al. 2013; Gray et al. 2014). This is achieved by linking genetic variation to environmental variables (Manel et al. 2012); however, in non-model species with limited or absent genomic information, the identification of adaptive genetic variation can only be achieved indirectly. By investigating anonymous loci in numerous individuals, it is possible to detect outlier loci of ecological relevance that may be linked to adaptive genes (Haldane 1948; Endler 1986; Schmidt et al. 2008; Stinchcombe and Hoekstra 2008; Manel et al. 2012). Their distribution differs from alleles at neutral loci, and their correlation to environmental influences reveals potential indications for adaptations (Holderegger et al. 2010).

Environmental association studies have formerly concentrated on the correlation between genetic and geographic distances to analyze IBD, while more recent studies also incorporate multiple environmental variables with a main focus on environmental gradients (e.g., climatic-, elevation-, environmental- and habitat gradients) with climatic variables being the most commonly used (Gerber et al. 2004; Nahum et al. 2008; Nakazato et al. 2008; Manel et al. 2012; Jones et al. 2013; Gray et al. 2014; Harter et al. 2015). Environmental association studies by means of climatic variables generally comprise temperature and precipitation that over large geographic distances normally change simultaneously. Thus, it is often impossible to detangle one factor from another to determine the exact driving forces for putative adaptation (Linhart and Grant 1996). As a consequence, only a few studies could detect outlier loci in non-model plant species. In Campanula (Jones et al. 2013) and Cotinus (Lei et al. 2015), these outlier loci are clearly associated with precipitation, while Manel et al. (2012) in alpine plant species, Gray et al. (2014) in Andropogon, and Hübner et al. (2009) in Hordeum showed that precipitations in combination with temperature were the best environmental predictors.

Here, we investigated the genetic diversity and structure of range edge populations of the annual Geropogon hybridus (L.) Sch.Bip. along a steep precipitation gradient (450–300 mm) on a micro-geographic scale (45 km) without any significant temperature change. Using amplified fragment length polymorphism (AFLP) analysis (Vos et al. 1995), we asked the following questions: (1) Does genetic diversity and differentiation change gradually toward range edge? (2) Can we identify significant IBD or IBE pattern despite the micro-geographic scale of the study area? (3) Is there putative precipitation-related adaptation among the surveyed populations?

Materials and methods

Study species

Geropogon hybridus is a diploid annual herbaceous Asteraceae species, 10–40–(80) cm high. The flowering stem is erect, usually glabrous with narrowly linear, grass-like leaves. The pedicule is hollow or swollen below the capitula. The linear, long involucral bracts often exceed the capitulum. The zygomorphic flowers (from March to May) have pink to violet corollas, with dark purple anthers. The species is self-compatible. It is pollinated by a variety of generalist insect families that are abundant in the landscape. Achenes are dimorphic, narrowly fusiform, or cylindrical in outline. The outer achenes are smooth with a long awened scabridulous pappus, most likely adapted to endozoochorous dispersal (mammalians). The inner achenes are prominently ribbed with a long, plumose pappus consisting of 20 unequal bristles, most likely adapted to anemochorous dispersal (Bobrov and Tzvelev 2000). The species is widely distributed throughout the Mediterranean region and toward southeast Asia (Danin 2016, ICN 2009+), but occurrences are restricted to semidesert or more humid environments.

Study region and sampling

The study region is situated between the cities of Kirjat Gat and Be’er Sheva (31°24′00″–31°14′50″N, 34°48′30″–34°50′30″E) in the southern Judea Lowland, Israel (Fig. 1) and represents the southernmost range edge of G. hybridus in Israel (Giladi et al. 2011, Danin 2016). The northern part of the area is influenced by Mediterranean climate with an annual rainfall of 450 mm per year, whereas the southern part is situated in a much dryer, semiarid region with only 300 mm of annual rainfall (Ben-Gai et al. 1994; Goldreich 2003; Giladi et al. 2011). In contrast to the strong precipitation gradient, no significant temperature changes can be observed on this small micro-geographic scale (Goldreich 2003; Giladi et al. 2011).

Fig. 1
figure 1

Map and sampling sites of the study region in Israel. In the magnified sections are the Isohyete (2011; worldclim.org), gray color depicts natural habitats and white color indicates agriculture fields. Cake diagrams at the right show the cluster of the STRUCTURE analysis (∆K = 2) for each population. In the left box, the distribution of Geropogon hybridus in Israel is shown [GBIF.org (29th July 2016) GBIF Occurrence Download http://doi.org/10.15468/dl.d2zwzr]. Map was prepared with ArcGIS Desktop (ArcGIS Desktop 10.2.2., Esri)

The area is characterized by a fragmented agroecosystem with scattered mosaic-like patches of natural vegetation, where insulated rock layers appear close to the topsoil (Svoray et al. 2007; Yaacobi et al. 2007; Gemeinholzer et al. 2012). On these patches, we sampled a total of 12 G. hybridus populations, with four populations each in the north (N1–N4), in the Center (C1–C4), and in the south (S1–S4) of the investigation area (Fig. 1; Table 1). The distance between the different regions are N to C ~12 km, C to S ~10 km, and N to S ~23 km (distances refer to the shortest distances between each region). Additionally, for some analyses, populations were grouped according to precipitation gradient (450 mm: N1–4; 400–350 mm: C1–4 + S1–2; 300 mm: S3–4). At each site, we randomly sampled 12 individuals; however, in some cases, evaluable data could only be obtained from less (Table 1; actual number of samples per population are provided). Samples of fresh and healthy leaves were randomly collected in the field and immediately dried in silica gel and stored at 4 °C until further processing. Total population size was estimated by counting the individuals in a plot of 150 × 150 m.

Table 1 Overview of sampled Geropogon hybridus populations and estimates of genetic diversity

Laboratory work and data analyses

DNA isolation and subsequent AFLP analysis were conducted as described in (Lauterbach et al. 2011). After an initial primer screening on eight samples (N1_1; N3_1; C1_1; C1_3; C3_1; C3_3; S2_1; S4_1) from the three subgroups each to exclude ascertainment bias (N, C and S) testing 12 different primer combinations, we chose three selective primer combinations that proved to be informative and provided clear bands which were sufficiently polymorphic to show variation within and among populations: EcoRI_AGC/MseI_CCA, EcoRI_AGC/MseI_CTC, and EcoRI_ACA/MseI_CGA. The laboratory setup was designed in a 96-well format. Thus, each primer could be analyzed in one run, assuming that AFLP produces clear and reproducible bands (Jones et al. 1997); therefore, we did not add replicates for error rate estimation. AFLP data matrix establishment was carried out using Genographer, version 2.1.4 (Banks and Benham 2008). Rare bands were included in the analysis (Online Resources 1 and 2).

Genetic diversity of the surveyed populations was estimated as percentage of polymorphic loci (PLP) and as Nei’s gene diversity (H e ) (Nei 1987) using AFLPsurv 1.0 (Vekemans 2002). To test for significant differences of diversity estimates between groups of populations (i.e., N vs. C vs. S), we applied a one-way analysis of variance (ANOVA) in R 3.0.2 (R Development Core Team Team 2013) using the Shapiro–Wilk test to check for normality (P = 0.41) and the Tukey HSD post hoc test to explore significant ANOVA results. The Bartlett test (P = 0.21) indicated no evidence for non-constancy of variance (heteroscedasticity) in our data.

Patterns of genetic population structure were visualized with a principal component analysis (PCA) using the R package ADEGENET version 1.4-2 (Jombart 2008). Additionally, we explored individuals’ genetic affiliation to genetic clusters using STRUCTURE version 2.3.3 (Pritchard et al. 2000). We explored individuals’ genetic affiliation to genetic clusters by applying the admixture model, 100,000 MCMC replicates, with a burn-in period of 50,020 and repeats per run for each chosen cluster number (i.e., K = 1–12), PLOIDY = 2 and RECESSIVEALLELES = 1. For all other settings, the default options were used. To identify the most likely K modal distribution, delta K (Evanno et al. 2005) was determined using STRUCTURE HARVESTER (Earl and von Holdt 2012). Corresponding graphs were constructed with DISTRUCT (Rosenberg 2004).

Genetic variation among groups of populations (ΦCT), among populations within groups (ΦSC), and within populations (ΦST) was partitioned with hierarchical analyses of molecular variance (AMOVA) using ARLEQUIN 3.5.1.2 (Excoffier and Lischer 2010). Additionally, pairwise ΦST values were estimated among populations. Significance levels were determined after 9999 permutations.

To evaluate patterns of isolation by distance (IBD) and isolation by environment (IBE), we tested for the relationships between pairwise ΦST values and geographic distance, and pairwise ΦST values and difference in annual precipitation, respectively. Pairwise and partial Mantel tests were conducted using the ‘vegan’ library in R (Oksanen et al. 2007) with 999 permutations. For IBD analyses, we used Euclidian geographic distances, and for IBE analyses, we constructed a matrix displaying pairwise precipitation distances (i.e., the difference of annual rainfall in mm between pairs of populations). For both types of analyses, we always tested (i) the complete data set of all populations and (ii) the two pairs of adjacent groups of populations (i.e., N–C and C–S, respectively).

To detect signatures of selection that could indicate putative adaptations to particular conditions along the gradient, we applied two differentiation-based genome scan methods: BAYESCAN 2.1 (Foll and Gaggiotti 2008) and DFDIST (Beaumont and Balding 2004) as included in the workbench MCHEZA (Antao and Beaumont 2011). BAYESCAN analyses were run with a burn-in of 50,000 iterations, a sample size of 10,000, and a thinning interval of 50, resulting in a total of 550,000 iterations. An additional burn-in was carried out by 20 short pilot runs of 5000 iterations. DFDIST analyses were conducted with 50,000 simulations, using the combined ‘Neutral mean FST’ and ‘Force mean FST’ algorithms and correcting results for multiple comparisons by setting the false discovery rate to 0.1. To consider a locus under selection, we chose a conservative approach. In BAYESCAN, only loci with a posterior probability over 0.99 and that were detected with a false discovery rate (FDR) of 0.01 were considered to be putatively under divergent selection. Comparably, in MCHEZA, only outliers that were identified at the 99% confidence interval were treated as putative adaptive loci. To test if adaptations occur either gradual along the gradient or are linked to a specific precipitation threshold (or any other biological significant line) for each genome scan method, we conducted two types of analyses: (i) a ‘global analysis’ including all populations, and (ii) ‘N–C specific’ and ‘C–S specific’ analyses including either only northern and central or central and southern populations.

To further substantiate the results of the two differentiation-based genome scans and to more directly test for signatures of precipitation-related adaptation, we applied the spatial analysis method (SAM; Joost et al. 2007), which uses multiple univariate logistic regression to test for correlations between environmental variables and binary molecular data. In our case, the site-specific mean values of mean annual precipitation in mm (Table 1) were assigned to the AFLP data. SAM uses the individual as reference unit, functions independently of any presumed population structure and is largely assumption free (Joost et al. 2007). Only if the two statistical tests implemented in SAM (likelihood ratio G and Wald test) reject the null hypothesis, a model is considered as significant (Joost et al. 2007). Also in SAM, the Bonferroni correction of the significance level for multiple comparisons was applied, which corresponds to a 99% confidence interval.

Results

Genetic diversity and differentiation

AFLP analyses with three different primer combinations and 91 analyzed individuals resulted in 123 unambiguously scorable polymorphic loci, ranging from 50 to 450 base pairs. The mean genetic diversity of the investigated Geropogon hybridus populations was H e  = 0.26, with genetic diversities ranging from H e  = 0.21 (S4) to H e  = 0.32 (C4, Table 1). The region specific mean values of diversity estimates (Table 1) varied from H eN = 0.26 to H eC = 0.29 with the lowest values in the south (H eS = 0.24). Genetic diversity among the three regions differed significantly: F = 4.33, P = 0.048. Post hoc tests showed a difference between central and south populations (P < 0.05). No differences in gene diversity estimates among precipitation specific groupings could be detected.

Principal component analysis depicted two large and distinct cluster that separated northern individuals from central and southern individuals along the first axes (Online Resource 3). Notably, one of the northern and three of the central individuals reflected contrary genotypic patterns, resembling those of the central or northern cluster, respectively. However, none of the southern individuals featured a northern-specific genotypic pattern. Overall, the first three principal components accounted for 31.4, 7.2, and 5.3% of the genetic variation.

The STRUCTURE analysis clearly confirmed the PCA results depicted two large and distinct clusters that separated northern individuals from central and southern individuals, detecting an optimal cluster number of ΔK = 2 (Fig. 2, Online Resource 4) that separated populations between the northern and central–south region (Fig. 1). Only few individuals of the northern and central site deviate from their respective group specific genetic pattern (i.e., in N2, C1 and C4).

Fig. 2
figure 2

Population genetic structure of 12 Geropogon hybridus populations as revealed by STRUCTURE analysis using the admixture model for ∆K = 2. Each individual is represented by a vertical bar, and fractional memberships in each of the clusters are indicated by colors

Analysis of molecular variance (AMOVA) resulted in a global ΦST of 0.29 (Table 2) and population pairwise ΦST values between 0.000 and 0.406 (Online Resource 5), indicating strong genetic differentiation between some of the populations (i.e., N1–4 vs. C1–4, and S1–4). Hierarchical AMOVA (Table 2) showed that 33.1% of genetic variance resided between study regions, whereas only 1.9% resided among populations within study regions. However, most variation was partitioned within populations (65%). When separating populations only in two regional groups, for ‘N versus C + S’ a total of 42% and for ‘N + C versus S’ only 13.5% of variance resided between groups (Online Resource 6), further depicting that the main population differentiation is located between northern and central sites. Simple Mantel tests revealed strong and significant IBD and IBE patterns (Table 3) for both, the ‘complete’ and the ‘north-central’ data set (Mantel r values between 0.76 and 0.86) but not for the ‘central–south’ data set (Mantel r = 0.29 and −0.02). To further unravel the correlated effects of IBD and IBE, we applied partial Mantel tests. These still revealed significant relationships for the complete data set (IBD controlling for precipitation: r = 0.54; IBE controlling for geographic distances: r = 0.37) but not for the two reduced data sets (Table 3).

Table 2 Results of simple and hierarchical analysis of molecular variance (AMOVA) performed by grouping the 12 Geropogon hybridus populations according to their geographic location along the surveyed precipitation gradient
Table 3 Results of simple and partial Mantel tests for pairwise ΦST with geographic and precipitation distance matrices of the surveyed Geropogon hybridus populations

Outlier detection

The two differentiation-based genome scan approaches detected 15 (DFDIST) and 12 (BAYESCAN) putatively adaptive outlier loci, or loci which might be in the same linkage group as loci under adaptation, for the complete data set, respectively. Of these, 11 were similarly identified by both approaches (Table 4, Online Resource 7). Whereas, for the ‘central–south’ data set, none of the two approaches detected any outliers, for the ‘north-central’ data set, at least BAYESCAN identified nine loci to be potentially under divergent selection. Almost, all of these ‘north-central’-specific outliers (8/9) were similarly detected by DFDIST and BAYESCAN in the complete data set (Online Resource 7), indicating that putative differences in the selection regime are located between the north and central sites and do not gradually occur over the complete north–south gradient.

Table 4 Numbers of outlier loci for populations of Geropogon hybridus as assessed by differentiation-based (BAYESCAN, DFDIST) and correlation-based (SAM) genome scan analyses

For correlation with mean annual precipitation data, SAM analysis revealed a total of 39 AFLP loci that showed a significant association at the 99% confidence interval. Interestingly, all 11 outliers that were similarly detected by the two differentiation-based genome scan approaches likewise were identified by SAM to be associated with the precipitation level (Table 4, Online Resource 7).

Discussion

Isolation by distance

Our population genetic investigation in fragmented range edge Geropogon hybridus populations along a steep precipitation gradient supports theory that range edge populations often feature reduced gene flow as result of the more fragmented distributions (Sexton et al. 2009). Other potential drivers explaining this pattern could be genetic bottlenecks (i.e., small population sizes), or selection differences. We detected potential IBD with the simple Mantel test (N–C–S r = 0.81; P = 0.001) and ΦST values (0.35; P < 0.001) indicative of a very strong differentiation along a short geographic distance. The findings are similar to the ones in Catananche lutea L., another annual Asteraceae species surveyed in the same investigation area (Gemeinholzer et al. 2012). This species co-occurs in the same vegetation along the same ecologic gradient under similar influences of fragmentation and also revealed variation to be greatest at largest distances, featuring an IBD scenario. Our analysis also supports the results of the summary analysis conducted by Sexton et al. (2014), who assessed gene flow with respect to environmental gradients in 110 published investigations. They found IBD to be the strongest pattern observed among plant studies.

IBD demands that the distribution of plots follows a linear model. The design of our sample setup in a natural environment of island like patches of natural vegetation in an agricultural surrounding takes this into account as far as possible. However, this strong overall population differentiation in our investigation was mainly driven by differences within northern populations, whereas central and southern populations alone showed no obvious IBD and only very weak population differentiation. Generally, significant IBD values can either be the result of continuous spatial variation in allele frequencies across the landscape, or they may derive from sharply divided groups where allele frequencies change rapidly at the point of division, due to environmental or dispersal limitations (Lindström et al. 2013) which may result in unaccounted overlapping effects of geographic barriers or potential bias connected to IBD (Meirmans 2012). Our results of the simple Mantel test (N–C) indicate significant barriers to gene flow and confirm restricted dissemination between northern and central populations, but not between central and southern populations, which is further corroborated by the STRUCTURE analyses. Meirmans (2012) argued, that by inducing allelic alterations in the sampling, the effect of geographic barriers may erroneously imply the occurrence of IBD (Silva et al. 2016) and IBD interpretations than, should be treated with caution.

Gene flow barriers may arise due to (i) factors restricting survival or persistence, or (ii) vector limitations in diaspore dispersal (pollen or seeds) as result of biotic or abiotic (physical or ecological) alterations. While dispersal limitation has been shown to limit species’ distributions at larger spatial scales (Eriksson 1998), its role in structuring communities at smaller scales has received less attention (Emery et al. 2009). No specific information about G. hybridus vector limitations in the investigation area is available, and the specific groups of insect pollinators are not known. However, composites are commonly visited by a wide variety of hymenoptera and coleoptera (Cheplick 1987; Bosch et al. 1997; Kaul et al. 2000) for many of which the distances on our micro-geographic mosaic-like vegetation patches are still in effective reach (e.g., Kwak et al. 1998; Morandin and Winston 2006). In addition, there are no seed dispersal distances for G. hybridus available; however, they are likely similar to the closely related and comparatively similar seed dispersed Hypochaeris radicata L. that has shown most seeds to spread within a distance of 100 m and rarely reaching distances up to 400 m (Soons et al. 2004). A population genetic-based maximization-of-explained-variance procedure resulted in a threshold distance of 3.5 km above which H. radicata populations were effectively genetically isolated (Mix et al. 2006). G. hybridus dispersal distances are probably quiet alike; however, isolation between fragmented habitat patches additionally aggravates the problem, not only because dispersal distance influences gene flow, but also due to the reduced likelihood that propagules reach suitable environments for establishment. The gradient-like structure of the project set-up with more or less evenly distributed mosaic-like populations in three regions across the investigation area was meant to account for that and clearly indicates that the statistically significant barrier of gene flow identified by both the simple and partial Mantel tests between N1–N4 and C1–S4 cannot solely be explained by geographic distances limiting dispersal. We found no explanation for the observed split into a northern and a southern gene pool based on soil type. Both locations, Galon in the north and Dvir in the south are characterized by grumosolic soil, with a higher content of loess in Dvir, whereas the central collecting point Lachish is situated on Redzina soil (Dan and Raz 1970).

In conclusion, other vector limitations or environmental influences must have contributed to the detected barrier between northern and central populations, either due to alterations in pollen or seed dispersal capacities or modifications in recruitment rates (Primack and Miao 1992; Nathan and Muller-Landau 2000; Clarke and Davison 2001; Givnish 2010). A possible hypothesis for the coherence of the central population Lachish with the southern population of Dvir might be the cultural history of Lachish. Situated on a strategic important hill, the area of Lachish was colonized around 6500 BC. Between 900 and 800 BC, Lachish developed into the second important garrison and residence town after Jerusalem (Utzschneider 2005: 46; Velikovsky 2009: 53). Located at the old caravan route between Syria and Egypt (Salvator 1879), Lachish was a center for trading with an own caravansary where most of the transport between Hebron and Gaza in south-east direction happened. In the 14th century, the large amounts of pilgrims from the Sinai region to Jerusalem also followed the old caravan route via Be’er Sheva to Hebron, laterthey made the pilgrimage by the old coast-way to Gaza and then north-east to Jerusalem (Robinson and Smith 1841: 736). In none of the records, we found any evidence for a transport route from the Lachish area to the north via Galon or Kiryat Gat. Thus, a possible explanation for the observed split into two genetic detectable entities might be an anthropogenic land use scenario, by which a dispersal limitation of the up to 3000 years isolation of the central and southern population (which for an annual plant is equal to 3000 generations) has led to a genetic drift at the edge of the distribution range. Caravanning from Lachish south to Be’er Sheva might have expanded the range of an edge distributed genotype to the southern area of Dvir, the opposite movement kept these populations in contact with the Lachish populations. Since no further route into the northern part of the species range existed, a further exchange between the larger northern populations and the smaller southern populations was hampered, leading to the now observable split within G. hybridus. In a further step, adaptation of the central and southern population to the low amount of precipitation might has occurred.

Theory predicts that IBD pattern often coincide with pattern indicative for niche evolution as result of strong selection and micro-geographic adaptation, or genetic impoverishment due to small population size and reduced gene flow from central populations (Ehrlich and Raven 1969; Nosil et al. 2005; Edelaar et al. 2008; Sexton et al. 2009; Urban 2011; Richardson and Urban 2013). IBD can especially be found in populations of fragmented habitats toward range edges (van Treuren et al. 1994; Young et al. 1996; Griffin and Barrett 2004). Eckert et al. (2008) reviewed 134 investigations to evaluate population structures across species’ geographic ranges by comparing central and range edge genetic population structures and found in 64% within-population genetic diversity declined toward range edges. In G. hybridus, however, nonsignificant lower H e values and equally high population sizes at range edge populations (Table 1) render genetic impoverishment unlikely and rather points to inter-population genetic divergence as result of increased random genetic drift, selection or adaptation (e.g., Young et al. 1996; Sork et al. 1999; Lowe et al. 2005; Aguilar et al. 2008). The mean genetic diversity of G. hybridusH e  = 0.262) over the whole study region was almost twice as high as in the comparable annual C. lutea (Gemeinholzer et al. 2012; Ø H e  = 0.136) in the same study area, but was comparatively low to investigations in fragmented solely out crossing H. radicata populations (Mix et al. 2006; Ø H e  = 0.88). Our values are similar to findings in other Asteraceae species across their whole distribution range, e.g., Leucochrysum albicans var. tricolor (DC.) Paul G.Wilson (Morgan et al. 2013, H e  = 0.13–0.39, with allozyme markers), Ixeridium dentatum subsp. nipponicum (Nakai) J.H.Pak & Kawano (Tanaka et al. 2014), Centaurea cineraria L. (Synonym: C. gymnocarpa Moris & De Not., Guarino et al. 2013; H e  = 0.027–0.567, with SSR markers) and Hieracium eriophorum St.-Amans (Frey et al. 2012), thus are comparatively high on this rather local scale. Similar high levels of genetic diversity were also found in Grevillea barklyana F.Muell. ex Benth (Hogbin et al. 1998) and Primula (P. veris L. and P. vulgaris Huds; van Rossum et al. 2004) in comparative analyses between central and range edge populations.

Isolation by environment

High genetic diversity in combination with IBD at species’ distribution ranges may be indicative of diversifying areas where ecological niche evolution can occur. This may happen in the course of adaptations to environmental alterations (Sexton et al. 2009; Mallet et al. 2014). Especially, at strong climatic gradients, environmental variability can contribute to patterns of interspecific genetic variation and can act as ecological driver promoting population divergence (Huang et al. 2015). The positive IBE pattern in G. hybridus indicates genetic variation to be correlated with reduced precipitation. However, as shown by partial Mantel tests, both (‘potential,’ see “Discussion” above) IBD and IBE are present in the population genetic data (Maassen and Bakker 2001; Shafer and Wolf 2013). Such patterns may occur when IBD is really strong and environmental spatial autocorrelation is low, as geographic distance then functions as suppressor to the IBE correlation. By conducting partial Mantel tests controlling for either geography or precipitation, we could demonstrate that both factors significantly contributed to the population genetic data. This pattern changed however when we analyzed the populations at both sides of the hypothesized gene flow barrier separately. Toward range edge, we found low IBD and no IBE correlations, thus the vast amount of unaccounted variation there must largely be due to other non-spatially structured biological or unmeasured environmental variables, and/or random processes triggered by ecological drift and dispersal (Legendre et al. 2009; Huang et al. 2015). In contrast, the northern populations featured highly positive IBD and IBE values with clear indications that the enhanced geographic distance here is more strongly affecting the population genetic pattern than precipitation thus IBD may drive IBE.

Several investigations confirmed high levels of divergence or local adaptation for populations distributed across climatic gradients, and particularly, steep environmental gradients can affect gene flow and local adaptation (Aitken et al. 2008; Leimu and Fischer 2008; Hereford 2009; Alberto et al. 2013; Savolainen et al. 2013; Rehm et al. 2015). Rundle and Nosil (2005) stated that divergent selections between distinct environments are the best understood drivers for speciation and population differentiation (see also Mallet et al. 2014). However, this mainly accounts for population genetic investigations across large geographic scales, while on micro-geographic scales, this is still relatively unexploited.

Outlier loci

By identifying outlier loci, we sought to determine how genetic differentiation is affected by sharp environmental gradients (Gray et al. 2014). We could identify 11 loci which putatively undergo diversifying selection and which moreover could be associated with precipitation levels. In the literature review about genomic heterogeneities, Nosil et al. (2009) found that many studies report 5–10% of loci to be outlier. This is in concordance with our analysis where 8.9% variation significantly correlates to the precipitation-related environmental gradient and which suggests that rainfall may put considerable spatially divergent pressure on the G. hybridus genotypes even on micro-geographic scales. However, unexpected was that potentially IBE associated outliers were found in the N populations instead of the marginal populations, indicative of the unexpected pattern that core populations may drive a potential IBE associated pattern.

Due to the applied anonymous marker technique AFLP, it is still an open question if the loci—or some of them—are indeed indicative of natural selection, or other adaptive and non-adaptive processes leading to non-random gene flow, or rather if they account for false-positive observations. Even though AFLP markers are frequently used for outlier detection, it is well known that the anonymous marker technique may cause statistical bias due to homoplastic and non-independent characters, especially for smaller fragments (Vekemans 2002; Bonin et al. 2007; Caballero et al. 2008). As AFLP markers typically are <500 bp, they most likely fall in noncoding regions; thus, markers with potential outlier behavior will be in close linkage with the potential genes under selection rather than inside the gene sequences itself, or they may act as regulatory elements (Stinchcombe and Hoekstra 2008; Butlin 2010; Nunes et al. 2012). AFLP genome scans of potential outlier loci mostly detect noncoding fragments, often with repetitive or transposable elements, or the characters derive from areas of low recombination which are often the result of non-adaptive processes (Minder and Widmer 2008; Wood et al. 2008; Nunes et al. 2012; Paris and Despres 2012; Cruickshank and Hahn 2014). In such cases, the outlier regions could simply reflect demography, spatially heterogeneous selection (Guttman et al. 2009), background selection in areas containing genomic features (i.e., centromeres), or sequence assembly artifacts (Shafer et al. 2015) and thus are indicative for standing variation, sometimes referred to as soft selective sweeps due to contemporary disruptive natural selection. On the other hand, outliers could be the results of new variants (hard selective sweeps), e.g., by chromosomal regions with reduced recombination (Strasburg et al. 2012). Linkage as result of background or purifying selection or genetic hitchhiking can additionally affect differentiation which is known to be stronger in coding than in noncoding regions and might account for false-positive or false-negative observations (Lotterhos and Whitlock 2015). The influence of genetic hitchhiking is dependent on multiple evolutionary parameters, such as recombination rate, population size, the strength and mode of selection, loci number, and analytical methods (Hermisson and Pennings 2005; Oleksyk et al. 2010; Pritchard and Di Rienzo 2010; Tiffin and Ross-Ibarra 2014; Lotterhos and Whitlock 2015). In our analysis, several indications point to selective sweeps derived from standing variation, due to the comparatively high mutation rates on this micro-geographic scale, and as theory predicts standing variation to support mutations to sweep rapidly to fixation (Hermisson and Pennings 2005). Selective sweeps result in comparatively low variation at linked sites and high divergence between populations of different selective environments (Strasburg et al. 2012), which is in concordance to the pattern found here in the range edge populations of G. hybridus. To minimize false-positive observations, researchers often rely on very high thresholds which increase the false-negative rate but strongly reduce the probability to wrongly identify loci in linkage with selective sweeps (Thompson et al. 2003).

Several other outlier tests have been developed beyond the outlier tests applied here (Beaumont and Nichols 1996; Vitalis et al. 2001; Beaumont and Balding 2004; Foll and Gaggiotti 2008; Excoffier et al. 2009; Bonhomme et al. 2010; Duforet-Frebourg et al. 2014) with several of them also accounting for genotype and environment associations by considering relations between allele frequencies and environmental parameters (Joost et al. 2007; Günther and Coop 2013; Gautier 2015; Lotterhos and Whitlock 2015). However, these analyses often suffer from high false-positive rates as they do not explicitly control for population structure in the test statistics (Sternberg et al. 2009; Meirmans 2012; De Mita et al. 2013; Lotterhos and Whitlock 2015). This might also be the case in our analysis, as populations in a landscape typically exhibit isolation by distance, and spatial autocorrelation in allele frequencies can cause false-positive associations between gene frequencies and the environment by chance. Natural populations would adapt to spatially heterogeneous environments at different spatial scales; therefore, geographically close populations in different selective environments would probably provide less false-positive hits than from distant populations where neutral demographic history might influence the pattern. Lotterhos and Whitlock (2015) state that the informative value of outlier-detections by associating genotype and environment associations is higher when populations physically near to each other (and therefore genetically similar for neutral genes) and yet in different environments (and therefore differentiated by selection) are being compared, like in our investigation.

Conclusions

The population genetic analysis of the annual Asteraceae Geropogon hybridus along a strong environmental gradient in a mosaic-like fragmented habitat revealed no significant population genetic diversity decline toward range edges. Contrary to expectations, population differentiation did not change gradually from the Mediterranean-influenced climate zone to the central and southern populations toward range edge in a more semiarid climatic region along a raster like sampling scheme. IBD and IBE were significant, despite the micro-geographic scale of the study area and partial Mantel tests controlling for either geography or precipitation showed that IBD and IBE significantly contribute to the population genetic divergent pattern over the whole investigation area. This pattern diminished when the hypothesized gene flow barrier was taken into account. The non-verge populations (N1–N4) featured highly positive IBD and IBE values with clear indications that the enhanced geographic distance here has more strongly affected the population genetic pattern than precipitation. Toward range edge (C1–S4), we found low IBD and no IBE correlations; thus, variation here must largely be due to other non-spatially structured biological or unmeasured environmental variables, and/or drift and dispersal.

By conducting environmental association studies, we could detect 11 outlier loci toward range edges which potentially evolved under selection, and which clearly indicate intraspecific population adaptive responses correlated to ecological drivers, in this case, precipitation. The results suggest that environmental factors can play prominent roles in population divergence, genetic drift, and directional selection, even on micro-geographic scales. This is especially the case along strong environmental gradients at species range edges when gene flow barriers and mosaic-like structures of fragmented habitats hamper dispersal. The results provide indications about potential niche evolution during adaptation. This highlights the significance of gene flow barriers, especially along fragmented range edge populations toward species’ ecological limits.