Environmental Heterogeneity Explains the Genetic Structure of Continental and Mediterranean Populations of Fraxinus angustifolia Vahl

Tree species with wide distributions often exhibit different levels of genetic structuring correlated to their environment. However, understanding how environmental heterogeneity influences genetic variation is difficult because the effects of gene flow, drift and selection are confounded. We investigated the genetic variation and its ecological correlates in a wind-pollinated Mediterranean tree species, Fraxinus angustifolia Vahl, within a recognised glacial refugium in Croatia. We sampled 11 populations from environmentally divergent habitats within the Continental and Mediterranean biogeographical regions. We combined genetic data analyses based on nuclear microsatellite loci, multivariate statistics on environmental data and ecological niche modelling (ENM). We identified a geographic structure with a high genetic diversity and low differentiation in the Continental region, which contrasted with the significantly lower genetic diversity and higher population divergence in the Mediterranean region. The positive and significant correlation between environmental and genetic distances after controlling for geographic distance suggests an important influence of ecological divergence of the sites in shaping genetic variation. The ENM provided support for niche differentiation between the populations from the Continental and Mediterranean regions, suggesting that contemporary populations may represent two divergent ecotypes. Ecotype differentiation was also supported by multivariate environmental and genetic distance analyses. Our results suggest that despite extensive gene flow in continental areas, long-term stability of heterogeneous environments have likely promoted genetic divergence of ashes in this region and can explain the present-day genetic variation patterns of these ancient populations.


Introduction
Understanding how environmental heterogeneity influences the distribution of genetic variation among natural populations along different spatial scales remains a central question in evolutionary biology and population genetics [1,2]. Genetic divergence of natural plant populations can be influenced by several evolutionary processes including gene flow, genetic drift, and natural selection [3]. If gene flow is locally restricted because of limited pollen and seed dispersal of the species, then the genetic differentiation of populations will show a pattern of isolation-bydistance (IBD) [4], which is considered to be the main force in the establishment of neutral genetic structure in plant populations.
Yet, greater genetic divergence than expected among populations inhabiting different environments has been used to suggest that contrasting ecological conditions may have a strong influence on the genetic differentiation of local populations [5,6]. Several studies have shown statistical associations between putatively neutral genetic variation and environmental variation in plant species, and such correlations may be interpreted as evidence of diversifying selection acting over the whole genome [3,[7][8][9]. In cases in which the genetic distance between populations correlates with their environmental distance, the pattern has been described as ''isolation by environmental distance'' (IBED) [10]. Nevertheless, the removal of geographic effects is necessary to detect the unique contribution of environmental gradients to genetic divergence because climatic differences are typically correlated with geographic distance [11]. A significant positive partial correlation after removing the effect of geographic distance suggests that genetic divergence is associated with environmental gradients and that natural selection may interact with neutral processes of gene flow and genetic drift [9][10][11].
Ecological niche modelling (ENM) allows the generation of biogeographical hypotheses and, when coupled with genetic data, provides new insights into the evolutionary history of animal and plant species [11][12][13][14][15]. For instance, ENM indicates that two population sets of unresolved taxonomic status with non-overlap-ping ecological niches and separated by a portion of unsuitable habitat may represent distinct evolutionary lineages ( [16,17], see example for geckos in Madagascar in [12]). As such, the application of ENM to delimit cryptic species has become an emergent field of ecological genetics [11,17,18]. In particular, the ability of ENM to test for a potential lack of spatial overlap in subpopulations of the same species allows the generation of different gene flow hypotheses that otherwise would be difficult to formulate. For example, if population differentiation with neutral markers is low but ENM simulations produce non-overlapping distributions, one can conclude that gene flow is still present despite the contrasting habitats associated with the populations. On the other hand, if genetic structuring of neutral markers reflects the simulated spatial distributions, rates of gene flow must be interrupted to allow for population differentiation. Although linking ENM with genetic data is increasing rapidly in phylogeographic studies at broad scales, the potential of the ENM approach and GIS-based environmental data has rarely been explored in conjunction with population genetic variation at finer scales. Studies of landscape genetics [1,19], which investigate how environmental factors influence the spatial distribution of intraspecific genetic variation, provide a promising framework for a better understanding of the microevolutionary processes involved in population genetic differentiation. While the majority of landscape genetics studies remain focused on animals [20], the integration of environmental features to explain the genetic variation of plants lags behind.
For this study, we selected the narrow-leaved ash (Fraxinus angustifolia Vahl, Oleaceae), a wind-pollinated tree species with a predominantly Mediterranean distribution, but found frequently in habitats of both Continental and Mediterranean biogeographical regions of Europe. It occurs naturally throughout Southern and Eastern Europe, from Portugal in the west to the Black Sea in the east. Although the phylogeography and genetic structure of the closely related species, the common ash (Fraxinus excelsior L.), and the hybridisation process between the two species have been well studied [21][22][23][24][25][26], no attempts have been made to understand the genetic structure of F. angustifolia populations. Unlike F. excelsior, which is widely distributed in mixed deciduous forests all over Europe, F. angustifolia is a habitat specialist associated with surface and ground waters, and thus, its dispersal ability is more restricted. In Central Europe, the Pannonian Basin and the Balkans, it occurs mainly in lowlands, riparian and floodplain forests along large rivers and their tributaries (Drava, Danube, and Morava), where it forms large and continuous populations. The distribution of this species in the Mediterranean region is patchy and reduced to smaller and more isolated populations on drier sites at higher altitudes or on wetland sites [27,28].
There are several different points of view regarding the taxonomic status of this species, but in general, the prevailing opinion is that F. angustifolia has three subspecies restricted by geographical regions [27,29]: ssp. angustifolia (in the western Mediterranean), ssp. oxycarpa (M. Bieb. ex Willd.) Franco and Rocha Afonso (in East Central and Southeastern Europe), and ssp. syriaca (Boiss.) Yalt. (in Turkey and eastwards to Iran). Fukarek [30] used morphological differences to further subdivide F. angustifolia in Croatia and the surrounding area of the Western Balkans. Continental populations along the rivers and floodplains of the Danube River Basin were named Fraxinus angustifolia ssp. pannonica (Fuk.) Soó and Simon. This putative new taxon is morphologically closer to ssp. oxycarpa, whereas Mediterranean populations along rivers and wetlands of the Adriatic River Basin probably belong to the typical ssp. angustifolia. Such a division differs from the more accepted general geographical classification, suggesting a possible increased within-region differentiation in the Western Balkan area.
Although located at a relatively high latitude and with a relatively small size, Croatia has been highlighted as having one of the most genetically divergent forests among 25 European forests investigated based on chloroplast genetic diversity [31]. Moreover, the area of the Western Balkans and Dinaric Alps served as an important refugium during the Pleistocene for the survival of many animal and plant species, including ash [23,[32][33][34]. Such large genetic divergence may have originated because of the comparatively higher environmental stability of this area during Quaternary climate oscillations, complex historical demographic events, and its geographical position. In fact, Croatia is situated along the contact line of three different biogeographical regions of Europe with contrasting climates: the Continental region including parts of the Pannonian lowlands, the mountainous Alpine region including parts of the Dinaric Alps and the Mediterranean region. Such environmental, landscape and historical diversity in a small area represents a valuable opportunity to investigate the influence of these factors on genetic variation and possible differentiation within species lineages.
In this study, we combined population genetic analyses of neutral loci, landscape genetic analysis using multivariate environmental data and ENM to examine genetic variation and its ecological correlates in F. angustifolia populations. We focused on answering these specific questions: 1) What are the levels of genetic diversity and divergence within and among natural populations of Fraxinus angustifolia distributed across the Croatian Continental and Mediterranean areas? 2) Does the neutral spatial genetic variation correlate with the environmental variation? and 3) Are areas predicted as suitable by ENM concordant with patterns of population genetic divergence?

Ethics statement
Collections of samples from protected areas were permitted by the authority of The Krka National Park and Lonjsko Polje Nature Park. For other locations no specific permits were required for the described field studies because sample collection did not involve endangered or protected plant species or privately-owned locations.

Study site and sampling
This study was conducted across the entire species' natural distribution range in Croatia (<56.538 km 2 land surface). Sampling was carried out in 11 natural populations of F. angustifolia ( Figure 1), with about 30 individuals sampled per population (total n = 345) ( Table 1). To avoid the sampling of close relatives, the minimum distance between sampled individuals was at least 50 m. Coordinates were recorded for each sampled tree using GPS. Trees were sampled from seven Continental and four Mediterranean natural stands, comprising the whole environmental gradient in Croatia in which the species occurs ( Figure 1). All populations from the Mediterranean region where stands had sufficient size to allow at least 30 individuals to be sampled were included.
About 9% of Croatian land area is frequently flooded [35], and natural floodplains of the Continental region represent the main habitat for the species. Floodplain forests in Croatia, dominated by narrow-leaved ash or mixed stands with oaks, are one of the most preserved in Europe. The majority of these forests stretch along the Sava, Kupa, Drava, and Danube rivers, all belonging to the Black Sea catchment (35.133 km 2 ). In this region, stands are  Table 1

Microsatellite analysis
Total DNA was extracted using DNeasy 96 Plant Kit (Qiagen) from up to ten mg of dried leaf tissue following the manufacturer's protocols. We used six polymorphic microsatellite markers that had been extensively used in Fraxinus spp. studies: Femsatl4, Femsatl11, Femsatl12, Femsatl16, Femsatl19 and M2-30 [36,37]. For Femsatl12, the redefined primers of Gerard [24] were used to avoid null alleles seen with the original primer set. The PCR conditions followed those used by Morand [22]. Fluorescent labelling of the forward primers enabled the detection of PCR products by capillary electrophoresis (ABI 3100), and allele sizes were scored using GeneMapper 4.0 (Applied Biosystems).
Population genetic diversity and structure GENEPOP 4.0 [38] was used to estimate the following genetic diversity parameters: average number of alleles per locus (N a ), observed heterozygosity (H O ), and expected heterozygosity (H E ). The program MICRO-CHECKER [39] was used to check for potential problems related to allele dropout and the presence of null alleles. The estimates of the null allele frequencies were based on the expectation-maximisation algorithm [40] and then calculated using FREENA [41]. The adjusted allele frequencies were used to recalculate the expected heterozygosity values (H Enull ). GENEPOP was also used to test for Hardy-Weinberg equilibrium (HWE) for each locus in each population and to test the loci for linkage disequilibrium. The probability tests were based on the Markov chain method [42]. The sequential Bonferroni adjustment [43] was applied to correct for the effect of multiple tests using SAS (SAS Ver. 9.1; SAS Institute Inc., Cary, NC, USA). FSTAT Ver. 2.9.3.2 [44] was used to calculate allelic richness (N ar ), which yields allele counts standardised to the minimum sample size, and to test the significance of the differences in average values of N ar , H O and H E between the Continental and Mediterranean populations. The number of private alleles (N pr ) per population was assessed by MICROSAT [45].
Pairwise genetic distances between populations (F ST ) and their significance were calculated using FSTAT. Pairwise F ST values were also estimated after correcting for the presence of null alleles using a method implemented in FREENA [41]. The overall population genetic structure was estimated for each locus and as a multilocus estimate with Wright's F-statistics using Weir and Cockerham's method [46] implemented in FSTAT. The analysis of molecular variance (hierarchical AMOVA) was performed to examine the partition of microsatellite variation between the Continental and the Mediterranean regions, among populations within the regions, and within populations using Arlequin Ver. 3.5.1.2 [47]. The variance components were tested by nonparametric randomisation tests using 10 000 permutations.

Association between genetic variation and environmental heterogeneity
Species presence data. We obtained 335 occurrence points from 11 sampled populations. Further locality records were obtained from the Flora Croatica Database (FCD, http://hirc. botanic.hr/fcd/, n = 60), the National Forest Inventory (n = 352, [48]), and personal communications (n = 55, see acknowledgements). In total, we compiled 802 high resolution species occurrence points.
Environmental data. Climate data for current conditions were obtained from the WorldClim database with a spatial resolution close to a square km [49]. First, the correlations among all 19 WorldClim bioclimatic variables and topographic variables for all presence points were calculated to exclude the highly correlated ones (r.0.75), whilst keeping the variables useful in predicting the distribution limits of trees, such as climatic averages and extremes [50].
Ten environmental variables were selected to describe the ecological characteristics of the sampled stands, for the Principal Component Analysis (PCA) and for the calculation of environmental distances. The eight bioclimatic variables included averages, extremes and seasonal variation in precipitation and temperature, and the two topographic variables altitude and terrain slope (see Table 2). Three additional layers were included in the construction of the ENMs as predictors: terrain aspect, a distance-to-water variable, and habitat type (see Table 3). Because of its circular nature, terrain aspect was recalculated and presented as two variables, northness and eastness [51]. Rasterised layers of distance-to-water and habitat types were obtained from the Croatian Wetlands and Habitat Map GIS databases maintained by the State Institute for Nature Protection (SINP; http://www. cro-nen.hr/map/index_en). All topographic variables were based on a 90-m spatial resolution digital elevation model (DEM) (Shuttle Radar Topography Mission; http://www2.jpl.nasa.gov/ srtm) and prepared in ArcGISH 9.3 (ESRI).
Correlation between genetic, geographic and environmental distances. To generate the environmental distance matrix, we performed a canonical discriminant analysis (CDA) based on the ten environmental variables (see Table 2) using PROC CANDISC in SAS. Squared Mahalanobis distances (D 2 ) between the populations were computed to obtain a matrix of environmental distances among the populations. Mahalanobis distances are analogous to Euclidian distances but also account for covariance among variables. Mantel tests [52] were used to examine the extent to which the neutral genetic structure can be explained by the environmental heterogeneity. We computed and tested the correlations between (1) the matrix of the natural logarithm of geographical distances (in km) between pairs of populations and the matrix of pairwise F ST /(1-F ST ) ratios and (2) the matrix of environmental distances (D 2 ) and the matrix of pairwise F ST /(1-F ST ) ratios. Only individuals that had information about all three parameters (genetic, geographic and environmental) were used for the correlation tests (n = 335). In addition, a three-way Mantel test was applied between the matrix of environmental distances and the matrix of pairwise F ST /(1-F ST ) ratios while accounting for geographical distances among populations. The significance level was assessed after 10,000 permutations as implemented in NTSYS-pc Ver. 2.02 [53]. Finally, the relationships between the populations based on both genetic distances and environmental distances (D 2 ) were visualised by constructing two neighbour-joining trees. Pairwise Nei's standard genetic distances [54] were calculated and an unrooted phylogenetic tree was constructed using the Neighbour-joining algorithm with 1 000 bootstrap replicates over microsatellite loci as implemented in the software PHYLIP Ver. 3.6b [55].

Ecological niche analyses
The total set of 802 occurrence points was used to further examine the levels of ecological niche divergence between the Landscape Genetics of Fraxinus angustifolia PLOS ONE | www.plosone.org populations from the Continental and Mediterranean biogeographical regions. Each occurrence point was assigned to a specific region as shown in Figure 1. First, we conducted the PCA based on ten environmental variables ( Table 2) that describe the ecology of all presence localities using PROC PRINCOMP in SAS. Second, we generated an overall ENM based on all 802 presence points. In addition to the ten variables used for the PCA, terrain aspect, distance-to-water, and categorical habitat type variable were added (Table 3). All environmental layers were resampled to be used at a 100 m6100 m spatial resolution.
We applied a maximum entropy presence-only modelling technique to estimate the ecological niche of the species using   [56,57]. This method has proven robust for presence-only data [58]. We performed ten replicate runs using cross-validation with default parameters and we used a logistic output from Maxent [57]. Model performance was evaluated using the area under the curve (AUC) of a receiver-operating characteristic (ROC) plot. To depict suitable habitat maps, we used the minimum training presence threshold [56]. Finally, the relative contributions of the environmental variables to the Maxent model were recorded.
To further assess whether populations from different regions occupy divergent ecological niche space, we simulated two separated ENMs using only occurrence data from either the Continental or Mediterranean region following the same modelling procedures. To determine whether any areas of overlap between the putative ecotypes exist, we summed the probabilities of occurrence from the two regional ENMs after applying the minimum training presence threshold to each. In this way, we evaluated whether the stands present in the predicted overlap zones were congruent with low values of pairwise genetic distances. All resulting models were visualised in ArcGISH 9.3.

Within population genetic diversity
A total of 176 alleles were observed across the six markers, with the number of alleles per locus ranging from nine (Femsatl16) to 48 (M2-30) and a mean value of 29.33. N ar ranged from 9.81 to 17.15 (Table 1). High levels of both within populations H O and H E were found (mean values over loci and populations were 0.683 and 0.699, respectively). There was no evidence of allele dropout in the data according to MICRO-CHECKER. Null alleles were suggested in seven out of 66 locus6population combinations. Estimated null allele frequencies using FREENA ranged from 0.052 (Femsatl12 in population Trnjani) to 0.146 (Femsatl16 in population Lonjsko Polje). The H E values increased slightly when recalculated using adjusted allele frequencies (Table 1), but no significant differences were observed between values of H E and H Enull in any of the analysed populations (Kruskal-Wallis test, P = 0.52-0.81). Therefore, all subsequent analyses were conducted using the original data set. Significant differences (P,0.05) in genetic diversity (mean N ar , H O and H E ) were found between the Mediterranean and Continental populations, with lower values observed in the Mediterranean stands. Moreover, 57 private alleles were identified in the Continental region, whereas there were only nine in the Mediterranean region (Table 1). No significant departures (P,0.01) from the HWE were observed at any loci in any population after applying sequential Bonferroni corrections. Finally, among a total of 165 tests for linkage disequilibrium between pairs of loci, no test was found significant after applying sequential Bonferroni corrections (P,0.01).

Genetic structure among populations and biogeographical regions
Although testing for HWE within each population showed no significant departures, the average multilocus inbreeding coefficient of the overall sample was slightly positive but significant (F IS = 0.024, P = 0.0025). Moreover, mean multilocus values of F IS in the Mediterranean region were almost four times higher and significant than that for the Continental region ( Table 4). The overall multilocus differentiation among all populations (F ST ) was 0.022. The within-region F ST , however, was higher in the Mediterranean populations (F ST = 0.027) than that among the Continental populations (F ST = 0.012), showing that the coastal populations were more structured.
Pairwise F ST values ranged from zero between Trnjani/ Ž upanja to 0.074 between Č akovec/Neretva population pairs (Table 5). No significant differences were observed between raw pairwise F ST and pairwise F ST corrected for null alleles (Kruskal-Wallis test, P = 0.79), suggesting that null alleles did not affect this analysis. Most population pairs from the Continental region had non-significant pairwise F ST values, with the exception of the Č akovec population. In contrast, most population pairs within the Mediterranean were significantly differentiated. Pairwise differentiation was also detected between the Mediterranean and Continental populations, with the exception of the Istrian population Mirna, which was not significantly differentiated from several Continental populations. The AMOVA analysis (Table 6) showed that most of the genetic diversity was attributable to differences among individuals within populations (97.36%). However, a small but highly significant percentage of variation was explained by differences among populations within regions (1.72%) and by differences between regions (0.92%), confirming the geographic structuring of populations.

Association between genetic and environmental variation
The analysed populations showed both significant levels of IBD (r = 0.385, P = 0.026) ( Figure 2A) and even higher correlation between genetic divergence and environmental distance (r = 0.549, P = 0.004) ( Figure 2B). The correlation between genetic and environmental distances remained significant (r = 0.426, P = 0.002) even after accounting for the effect of geographical distance in a three-way Mantel test ( Figure 2C). On the other hand, the removal of the effect of environmental variation in the partial Mantel test resulted in a non-significant correlation between genetic and geographic distances (r = 20.025, P = 0.438). Therefore, our populations show a clear ''isolation by environmental distance'' pattern, rather than IBD as such. Ecological niche analyses PCA analyses. The first principal component (PC1) explained 60.32% of the total variation and clearly separated the Continental and Mediterranean localities along temperature and precipitation gradients (Figure 3). Elevation and slope were highly positively correlated with PC2, explaining 12.70% of the total variation and reflecting a topographic gradient. In addition, PCA on environmental data revealed a notable environmental sub-structure within the Mediterranean region, where each of the populations clustered along a different river valley. The results show that Continental populations occupy habitats that are cooler with lower winter temperatures (below zero) and higher summer rainfall. In contrast, coastal populations are characterised by warmer habitats with higher winter temperatures and wetter winters. Generally, the range of variation in the examined environmental variables appears to be more pronounced between the Mediterranean localities than among the Continental localities.
Ecological niche modelling. The Maxent model performed well with high average training and test AUC values across ten replicate runs (Table 7) and was congruent with the currently known distribution of F. angustifolia in Croatia ( Figure 4A). The highest probabilities of occurrence were in the Continental region, in lowlands along large rivers (Sava and Drava) with more or less continuous distribution. In contrast, the predicted distribution was discontinuous in the Mediterranean region with several isolated areas of high suitability associated with shorter karst river valleys along the eastern Adriatic coast (Mirna, Krka, Zrmanja, and Neretva). The overall model did not predict suitable habitats in the Alpine region, confirming that that the species distribution is strongly associated with river valleys and wetland sites ( Table 3).
The overlap between the two regional modelled distributions was very low (Fig. 4B), suggesting strong regional niche differentiation between the two putative ecotypes. Despite high AUC values, each model alone predicted a highly reduced distribution of the species in comparison with the overall model. Contrary to our expectations, the variables differed in their contributions to the three distribution models (Table 3). Distanceto-water contributed most to the overall and Mediterranean ENM, whereas habitat type and elevation were most important predictors for the Continental ENM.
Neighbour-joining analysis. Trees based on either genetic distances ( Figure 5A) or environmental distances ( Figure 5B) were congruent in their major features, suggesting that environmental variation may promote genetic divergence of the studied populations. Moreover, a comparison of the genetic tree ( Figure 5A) with the overlap map of the regional ENMs ( Figure 4B) showed that the coastal population Mirna, which is situated intermediately in the genetic tree, is located in the overlap area. There are, however, some incongruencies. For example, the Jastrebarsko population belongs to the Continental region based on genetic markers but can be considered intermediate from an ecological point of view.

Discussion
Our combined analysis of genetic data, multivariate statistics on environmental data and ENM suggests that current genetic variation patterns in natural Fraxinus angustifolia populations in Croatia may be influenced by the local ecological conditions rather than by geographic distances only. We observed an overall pattern of significantly higher genetic diversity in the Continental region and low local differentiation that contrasts with the reduced genetic diversity and stronger structuring in the Mediterranean region. The extent of potential ecological niche overlap between the continental and coastal populations was low, suggesting that two ecologically distinct lineages of narrow-leaved ash may occur in the study area. ENM was in agreement with the genetic distance   Table 2  analysis, which found that most Continental and Mediterranean populations were differentiated.

Genetic diversity and structure
The genetic analysis results agree with the expectations of high polymorphism within populations and low genetic differentiation between populations, as observed in ashes and wind-pollinated trees in general [21,26,59,60]. However, the large and significant local homozygote excess (F IS .0.15) found in many European F. excelsior populations (22,26,59) was not observed in the present study. This characteristic of common ash was often attributed to the presence of null alleles, biparental inbreeding or the Wahlund effect. Unlike common ash populations studied in Europe [21,22,26,59] all of the populations studied herein were in HWE, indicating weak evidence for the presence of null alleles or of local inbreeding. Finally, the average multilocus inbreeding coefficient F IS , which provides information on the cumulative effect of inbreeding, was low but significant at the level of the overall sample (average F IS of 0.024), but at the regional level only for the Mediterranean stands (average F IS of 0.046). Although positive F IS values may reflect the presence of null alleles, which are commonly suspected for microsatellite loci, controlled crosses carried out by Morand [22] showed that most of the loci used in this study do not show null alleles. Moreover, there were no significant differences between the H E and F ST values calculated with and without corrections for null allele frequencies. A plausible explanation for this positive inbreeding would be a Wahlund effect at the level of the overall sample population, mainly due to the differentiation between subpopulations, even though subpopulations themselves are in HWE.
Local population sizes determined by suitable habitat availability may explain the contrasting genetic structures found within and among biogeographical regions. The observed higher genetic diversity in the Continental populations can possibly be maintained due to larger effective population sizes compared to their Mediterranean counterparts. Longitudinal distribution along rivers and floodplains typical for the Pannonian lowlands with no barriers to gene flow allows free pollen and seed dispersion between the populations, which explains the lack of significant genetic structure in this region based on pairwise F ST , even between most distant populations. On the other hand, Mediterranean populations are reduced to a few smaller suitable sites associated with Dinaric karst fields, short karst rivers and rare natural wetlands with no apparent above-ground connections due to the limestone base. In consequence, populations are more isolated from each other and also from the Continental part of the distribution, limiting dispersion and favouring the maintenance of an intra-regional genetic structure.

Genetic divergence and environmental heterogeneity
Large amounts of neutral population genetic variation were explained by environmental variation rather than by simple geographic distances, suggesting a strong role of environmental heterogeneity in the genetic divergence of populations [6][7][8][9]61]. First, pairwise F ST values and multilocus estimates of F-statistics suggest that each region has different evolutionary constraints. For example, some very close Mediterranean populations (such as Zrmanja and Krka) have significant pairwise F ST values although they are not geographically distant (only 40 km apart), suggesting more restricted gene flow in this region probably caused by isolation of suitable habitat and/or habitat differentiation. Environmental distinctiveness and habitat discontinuities among Mediterranean populations are apparent from the PCA and niche models, revealing a similar differentiation pattern compared to pairwise F ST . The Mediterranean is known for pronounced environmental heterogeneity over very short distances because of factors such as slope, exposure, distance from sea, and rock type. In contrast, most of the Continental populations that occur in a rather homogenous environment exhibit non-significant pairwise F ST values (except for the Č akovec population), even in distant populations (.240 km). Similar trends were observed in Taxus baccata [62] on a broader geographical scale, in which populations located in the stronger Mediterranean climate displayed higher pairwise differentiation within regions than those from the inland areas, suggesting that the geographical and environmental features can influence population divergence of different tree species in this area.
Second, we found a significant correlation between genetic and environmental variation. IBED patterns may be arising from a neutral process of temporally disrupted gene flow among individuals living in environmentally distinct habitats, leading to phenological differences [7,22]. Gene flow should homogenise neutral genetic variation in wind-pollinated tree species at short geographical distances, but habitat differentiation can act as a barrier to gene flow, causing environmental isolation and genetic differentiation of spatially close plant populations [5], as found in our study. Finally, a significant impact of genetic drift in these populations can be discounted as drift alone would create a random genetic structure that was not observed herein; instead, ecologically similar populations were also grouped genetically. In sum, we show that the observed genetic variation pattern is associated with environmental gradients. We recognize that the observed IBED pattern does not imply causality and this correlation might have other plausible or more complex interpretations which cannot be explicitly tested using our current data, like different population ages in two regions because of independent colonization events. Future tests with candidate genes for traits of interest could also clarify the possible role of natural selection in shaping the divergence of these populations, but such markers are only recently emerging for Fraxinus spp. [63]. Further exploration in our study species is currently underway.

Ecological niche models and genetic structure
Using an ENM approach, we tested whether the predicted distributions of the species corresponded with the patterns of population genetic structure. In particular, we searched for areas where the Continental and Mediterranean ENMs overlap, as they could highlight populations from different biogeographical regions with lower levels of pairwise genetic differentiation. The overall ENM detected a separation of the two putative ecotypes by the mountainous region, representing an unsuitable habitat for the survival of the species and a potential barrier to gene flow. The Continental and Mediterranean ENM barely overlapped, indicating a clear divergence in the ecological niche space occupied by the populations in each region. An overlap of independent ENMs suggests only one point of contact in the Istrian peninsula at the Mirna population site ( Figure 4B). This stand is indeed the only Mediterranean population with non-significant pairwise genetic distances towards most Continental populations and is located within an area of intermediate environmental conditions between the coast and the continent. Because one set of populations poorly predicts the distribution of the other set of populations and because their ecological niches almost do not overlap, our populations may represent two distinct evolutionary lineages, even with the low levels of genetic divergence [12,16]. Lack of niche overlap in wide-ranging tree species may also appear to be driven by differences in abiotic conditions in different regions (soil, elevation, climate) [64]. Because our populations inhabit regions with clearly divergent climatic regimes and F. angustifolia has a wide distribution, both latitudinally and altitudinally, this is a confounding factor that needs to be kept in mind.

Stability, migration crossroads and environmental heterogeneity in refugia
From our observations, we cannot exclude the fact that historical migration processes in this refugial area could have raised the observed pattern of genetic variation in studied populations. Mediterranean tree populations have persisted in the southern refugia without significant geographical movements due to long-term stable environmental conditions in the Mediterranean [65]. Hence, the lower genetic diversity, higher genetic differentiation and higher fixation rates observed in the Mediterranean populations could result from older populations with smaller historic population sizes persisting in environmentally more stable regions over longer periods [66]. Croatia and the wider area of the Dinaric Alps have already been identified as an important refugium for ash and other temperate tree species during the Pleistocene and stand on a contact zone of their different postglacial recolonisation routes [34,67,68]. Recent studies have also confirmed the existence of a 'refugia within refugia' pattern in these areas, where differentiation of distinct lineages on a small geographical scale has been observed (see [69] and references therein). Therefore, forests in this region are expected to harbour greater regional genetic diversity and uniqueness in comparison with the rest of their range [31,65].
F. angustifolia is a thermophilic tree species with distinct moisture requirements that could have survived in situ during the last glacial period in several river valley sites along the Dalmatian coast at lower to mid-altitudes until today. These humid but not too cold sites provided continuous moisture availability and shelter for Mediterranean tree species during the Adriatic Sea level drop in the LGM, leaving the Northern half of the Adriatic Sea basin exposed and unsuitable for the survival of moisture-dependent species [33]. Northern coastal populations in Istria and the rest of the Continental populations could have been recolonised by expansion from the Dinaric Alps or from refugia in North Italy and/or Balkan Peninsula [68], most likely via the North Adriatic or along the Danube river lowlands. At least three F. angustifolia haplotype lineages meet at the vicinity of the investigated populations (H01, H02 and H03, sensu [68,27], authors' personal observations), confirming that various migration events occurred in the past within the study area and suggesting that Croatian populations may have originated from various colonising routes that likely brought new diversity. Higher levels of genetic variation in the Continent could therefore be due to a gene flow among individuals from different glacial refugia in newly colonised regions [31] while southern coastal populations probably represent relict, genetically more divergent populations [65]. In fact, modelling of potential distribution of Quercus robur in Europe during the LGM [14] shows that both the Adriatic coastal and Continental lowland parts of Croatia were suitable for the survival of this ecologically similar species in situ. Assuming that F. angustifolia followed a similar distribution during the LGM, this species likely survived for long time in this area, allowing enough time for the differentiation of distinct populations to occur through a processes of ecological isolation.

Conclusions
Overall, our results suggest that long-term stability of heterogeneous environments at regional spatial scales may explain current levels of genetic diversity and population genetic divergence in narrow-leaved ash in these ancient refugia. Environmental differences between the regions may have led to the general subdivision into two ecotypes, with the pronounced environmental heterogeneity in the Mediterranean further promoting the genetic differentiation of the coastal populations. Thus, the local genetic structure in the narrow-leaved ash is more complex than a simple allopatry divergence model as the populations are not clear-cut differentiated but rather in a complex genetic cline, probably resulting from the environmental heterogeneity over the studied geographical area.