Population Genetic Diversity of Quercus ilex subsp . ballota ( Desf . ) Samp . Reveals Divergence in Recent and Evolutionary Migration Rates in the Spanish Dehesas

The Spanish dehesas have been severely affected by human activities that date to the prehistoric period and have suffered accelerated decline since the 1980s. Holm oak (Quercus ilex subsp. ballota (Desf.) Samp.) is a key component of this system, and its acorns provide an important food source for wildlife and domesticated livestock. Our earlier work showed structured variation in acorn morphology and biochemistry. Here, we used chloroplast and nuclear microsatellites to detect genetic structure among populations of Q. ilex from the major biogeographic regions of Andalusia. We found high levels of spatial differentiation with chloroplast DNA indicating little seed dispersal among populations. Spatial differentiation was weaker for nuclear DNA, presumably as a result of more widespread pollen dispersal and its larger effective population size. The Baetic Cordillera (Cádiz) population consistently appeared well separated from populations of the northern Sierra Morena, suggesting that the Guadalquivir Valley has played an important role in determining population divergence. This may be, in part, evolutionary, as suggested by chloroplast DNA, and, in part, a result of human-induced population isolation, as Q. ilex has been removed from the Guadalquivir Valley. Evolutionary gene flow rates were greater than contemporary rates, which were limited to unidirectional gene flow from Córdoba to other populations in the Sierra Morena and, surprisingly, to the southern population at Almería. The inconsistency between evolutionary and recent migration rates suggests an effect of anthropogenic activity over the last few generations of Q. ilex.


Introduction
Holm oak (Quercus ilex subsp.ballota (Desf.)Samp.) is the dominant tree species in natural forest ecosystems over large areas of the western Mediterranean Basin and is the distinctive element of the Spanish agrosilvopastoral ecosystem "dehesa".Both holm oak natural forest and the dehesa are seriously threatened by several factors, the most important of which are of anthropogenic origin: overexploitation, poor management practices, and fire, as well as diseases and adverse environmental conditions [1].These have contributed to a severe decline observed since the early 1980s [2,3].Lack of regeneration that has led to populations dominated by aging individuals is one of the greatest threats to these ecosystems and is likely to become more problematic under anticipated climate change, with high temperatures and severe drought episodes expected for southern Spain and other Mediterranean countries [4,5].
Forest restoration and reforestation as well as sustainable management were major objectives for the afforestation programs implemented in Spain at the end of the 20th century and the beginning of the 21st, promoted by the European Union which established a Community aid scheme for forestry measures in agriculture lands [6].Q. ilex has become a priority species to achieve these objectives.The ultimate goal requires the selection of "elite" genotypes in terms of higher germination rates, increased acorn production, and desirable quality traits related to nutritional values, as well as adaptation to adverse biotic and abiotic stresses.This selection should be directed by phenotypic, physiological, and molecular analyses.Our group is carrying out a project in which variability in a number of traits in Andalusian holm oak populations is being analyzed.At the phenotypic level, we have shown the existence of differences in tolerance to drought stress and resistance to Phytophthora cinnamomi Rands infection among populations from different geographical locations; the differential response has been characterized at the physiological (water status and photosynthesis) and proteomic levels [7][8][9].In an attempt to analyze variability and differences among provenances of holm oak from Andalusia, acorn morphometry and Near infrared reflectance spectroscopy (NIRS) chemical analysis as well as protein profiling in pollen and acorns were performed [10][11][12].We found high levels of variability within and between populations.Although spatial grouping of populations was evident, spatially intermediate populations clustered differently according to the analysis employed (NIRS or proteomics), raising the question as to what extent population genetic structure can explain spatial variations in these phenotypic traits.
Here, we use molecular data to investigate whether regional structure that might be associated with past and contemporary gene flow can be detected.Any observed genetic structure would then form the basis for strategies for the selection of phenotypic traits and, possibly, the selection of seed sources for planting.High genetic diversity is expected for holm oak considering its very long life; fecundity, allogamous, and anemophilous reproductive traits; and promiscuity and long history of introgression with other species [4,[12][13][14].We hypothesize that genetic structure among populations within Andalusia could be significant because of its importance as a glacial refugium.For many species, the current spatial structure of genetic diversity is determined, in part, by past climatic changes that have led to population shifts or by demographic changes [15].The Iberian Peninsula is considered an important glacial refugial region within Europe [16,17].Several phylogeographic studies have shown complex patterns of sub-refugia within the Iberian Peninsula [18], of which the Baetic Cordillera recurs as a sub-refugium for a wide range of taxa, including tree species [19][20][21][22].A chloroplast phylogeography of Q. ilex revealed a unique chlorotype from southern Andalusia, suggesting an ancient refugial location [23].However, a multi-marker genetic study of the Iberian Peninsula suggested widespread gene flow in Q. ilex, with relatively little differentiation [13].Paleoclimate modeling of southern Spain has indicated that it was likely composed of multiple upland climatic refugia within the Baetic Cordillera and the Sierra Morena that could have been sites for species with different habitat requirements [24].
In view of the possibility of divergence among populations of Q. ilex in Andalusia and the likely partial genetic basis of our previously observed chemical composition and protein profiles [25], we undertook a molecular genetic study of populations of Q. ilex from the major regions of Andalusia.We performed nuclear and chloroplastic microsatellite analyses on 94 holm oak individuals from five natural populations distributed across bioclimatic gradients in the Andalusian territory.Our objectives in this research were to (1) determine levels of divergence among populations across Andalusia, (2) estimate differences in migration at an evolutionary and contemporary scale as well as the role of population bottlenecks on genetic diversity among the regions within Andalusia, and (3) assess Forests 2018, 9, 337 3 of 17 whether differences in morphological and biochemical traits that we have previously reported are associated with genetic divergence and gene flow within the region.

Populations, Plant Material, and DNA Extraction
The present study was performed with five populations covering the major bioclimatic regions of Andalusia, southern Spain: Southeast (SAA: Almería), Southwest (BCA: Cádiz), Northeast (PCO: Córdoba), and Northwest (APS: Seville, CTH Huelva) (Figure 1 and Figure S3; Table 1).Andalusia is characterized by two west-east trending mountain ranges: the Sierra Morena in the north and the Baetic Cordillera in the south (the latter composed of a Mediterranean coastal range, the Cordillera Penibética, and an inner range, the Cordillera Subbética).The two major chains are separated by a broad central depression, the Guadalquivir Basin.Today, Q. ilex occupies low to mid-elevations in the two major mountain ranges, northern and southern populations being separated by the Guadalquivir Basin.Geographical coordinates, altitude, climatic data, and soil characteristics of each surveyed area are shown in Table 1.Twenty individuals were sampled from each locality; however, for two localities, poor DNA amplifications resulted in fewer samples being analyzed (see Table 1 for sample sizes).Trees were sampled across the landscape to minimize the likelihood of relatedness through seed, and no trees were sampled that were less than 10 m apart.The present study was performed with five populations covering the major bioclimatic regions of Andalusia, southern Spain: Southeast (SAA: Almería), Southwest (BCA: Cádiz), Northeast (PCO: Córdoba), and Northwest (APS: Seville, CTH Huelva) (Figure 1 and Figure S3; Table 1).Andalusia is characterized by two west-east trending mountain ranges: the Sierra Morena in the north and the Baetic Cordillera in the south (the latter composed of a Mediterranean coastal range, the Cordillera Penibética, and an inner range, the Cordillera Subbética).The two major chains are separated by a broad central depression, the Guadalquivir Basin.Today, Q. ilex occupies low to mid-elevations in the two major mountain ranges, northern and southern populations being separated by the Guadalquivir Basin.Geographical coordinates, altitude, climatic data, and soil characteristics of each surveyed area are shown in Table 1.Twenty individuals were sampled from each locality; however, for two localities, poor DNA amplifications resulted in fewer samples being analyzed (see Table 1 for sample sizes).Trees were sampled across the landscape to minimize the likelihood of relatedness through seed, and no trees were sampled that were less than 10 m apart.For DNA extraction, a few leaf samples were collected from young shoots from the upper part of each tree, transported to the laboratory on ice, abundantly washed with tap water, blot dried, frozen in liquid nitrogen, and stored at −80 • C. Genomic DNA was isolated following the protocol reported by Echevarria-Zomeño [26] and diluted to 10 ng µL −1 to carry out PCR amplifications.

Nuclear and Chloroplast Microsatellite Analysis
Ten nuclear microsatellites (SSR) and 10 chloroplast SSR markers (Table S1) previously developed from other species, including Q. macrocarpa Michx., Q. robur L., Q. petraea (Matt.)Lieb., and Castanea sativa Mill., were screened in the populations collected in this study.PCR was performed in a 20 µL volume, with the reaction mixture containing 1× PCR buffer (Biotools, Madrid, Spain), 2 mM MgCl 2 , 0.2 mM dNTPs, 0.5 µM of each primer, 0.5 units Taq DNA polymerase (Biotools, Madrid, Spain), and 30 ng of DNA.Amplifications were conducted as follows: 1 cycle of 3 min at 95 • C, followed by 35 cycles of 1 min at 94 • C, 45 s at the annealing temperature indicated in Table S1 for each primer, and 1 min at 72 • C, followed by a final incubation of 7 min at 72 • C. The PCR reactions were carried out in a 96-well block thermal cycler (Applied Biosystems, Madrid, Spain).PCR products were detected using an ABI PRISM 3130xl Genetic Analyzer and GeneMapper analysis software (Applied Biosystems, CA, USA).For capillary electrophoresis detection, forward SSR primers were labeled with the 5 fluorescence dyes PET, NED, VIC, and 6-FAM, and the size standard used in the sequencer was Gene Scan 500 Liz (Applied Biosystems, CA, USA).Because of irregular amplification, the nuclear microsatellite locus, QpZAG9, was dropped from the analyses.

Chloroplast DNA
Chloroplasts are inherited maternally without recombination in oaks; thus microsatellite variants were combined into haplotypes, whereby each haplotype was considered to be inherited as a single allele.We estimated within-population haplotype diversity ( cp H S ), overall haplotype diversity ( cp H T ), global among-population differentiation by allele identity ( cp G ST ), and global among-population differentiation considering the allele size ( cp R ST ) using PERMUT software (Petit et al., 2002 [21]).We also used PERMUT to test if the observed cp R ST value was significantly different from cp G ST .This test is used to detect phylogeographic patterns (i.e., whether mutation has contributed significantly to population differentiation).
To explore divergence among haplotypes and their sharing among populations that would indicate evolutionary levels of successful seed migration, we constructed a haplotype network using the median-joining (MJ) network algorithm [27].MJ networks were post-processed with a maximum-parsimony (MP) algorithm [28] to remove unnecessary linkages and median vectors.Network construction was performed in NETWORK 5.0.0.1 (available as freeware from http://www.fluxus-engineering.com/sharenet.htm).

Nuclear DNA
We tested for null alleles using MICRODROP [29].MICRODROP uses an expectationmaximization algorithm to obtain joint estimates of allele frequencies, drop-out rates caused by sample factors, and locus factors and inbreeding coefficients to correct for deficits in heterozygosity.The software can provide imputed datasets by drawing genotypes according to the posterior distribution of the model and replacing false homozygotes by heterozygotes.We performed MICRODROP imputations separately by population.We tested for deviations from Hardy-Weinberg (HW) equilibrium within each population by the inbreeding fixation index F IS with the software FSTAT, version 2.9.3.2 [30].We estimated allelic richness with the rarefaction method (R t ) and expected heterozygosity (H e ) with FSTAT to analyze the level of within-population genetic diversity.
Population bottlenecks result in a transient disequilibrium by reducing allelic diversity through preferential loss of rare alleles, before having an effect on heterozygosity.We used BOTTLENECK 5.1 [31] to test for a signal of a bottleneck in all sampled oak populations under the two-phase model (TPM) with the default 70% single-step mutations and 30% multiple-step mutations.The TPM is an intermediate between the infinite allele model and the stepwise mutation model and is likely closer to the true mutation model of most microsatellites [31].Statistical significance was tested with the one-tailed Wilcoxon signed-rank test to determine whether observed heterozygosity deviated from expectations at mutation-drift equilibrium.Estimations were based on 10,000 replications.Reductions in population size were also tested using the "mode-shift" indicator of the distortion of allele frequency classes' distributions.In non-bottlenecked populations, the mode-shift distribution should be approximately L-shaped with the greatest number of alleles being detected at low frequencies [31].
To infer population structure, we used a Bayesian model based approach performed in the program STRUCTURE, version 2.3.4 [32] and a model-free multivariate approach, discriminant analysis of principal components (DAPC) [33].For STRUCTURE, an admixture model with correlated frequencies was used without prior population information.Preliminary runs revealed a complex hierarchical structure.We inferred population division (number of K populations) by performing 20 independent runs of each K (K = 1 to K = 5) with a burn-in of 100,000 iterations and 1 million iterations of the Gibbs sampler.Log-likelihood of the data was recorded for each run, and the ad hoc statistic ∆K was calculated following Evanno et al. [34].Output from STRUCTURE was post-processed to obtain ∆K and plots for publication using the online program CLUMPAK [35].
For DAPC, discriminant analysis (DA) was used to identify clusters on the basis of data that were transformed using principal component analysis (PCA).The number of PCA axes that were retained for DA was determined by 30 replicates of cross-validation.DAPC was performed in R, using the Adegenet package [33].
We used Barrier 2.2 [36] to test for the most prominent breaks in nuclear microsatellite data.Barrier detects geographic barriers by testing the correlation between the geographic distance and genetic distance among sampling locations using the Monmonier algorithm.We performed the Barrier analysis with 100 distance matrices (shared alleles) from bootstrapped data obtained by MSA 4.05 [37].

Inferring Migration
We used BayesAss [38] to estimate contemporary migration among the five sampled populations.BayesAss uses a Bayesian multilocus approach in which populations are not constrained by stationary conditions, so that individual populations may be out of HW equilibrium.Bidirectional migration rates among all populations are inferred over the last generations by detecting transient disequilibrium among multilocus genotypes of migrants or recent descendants of migrants.In preliminary runs, we adjusted mixing parameters to obtain acceptance rates in the range 0.2-0.6 (∆a = 0.5, ∆m = 0.5, and ∆f = 0.7).Final runs (five runs) were performed with different starting seeds, run for 10,000,000 iterations with 1,000,000 burn-in, and averaged to provide estimates of migration rates.Convergence of the runs was assessed by examining plots of total log-likelihood versus iteration to verify that there were no peaks or troughs using TRACER 1.5 [39].
Because contemporary and evolutionary migration rates can provide complementary information on population demography, we also estimated evolutionary bidirectional migration among populations using Migrate-N.We evaluated three population models that took into account the biogeography of Andalusia and the clustering of our earlier results from morphometrics and acorn chemistry: Model 1: Full model with 5 population sizes and 20 bidirectional migration rates.Model 2a: Two-population model that followed our population groupings based on acorn morphometrics and fatty acid composition that showed separation between Sierra Morena populations (Huelva, Córdoba, and Seville) from Baetic populations (Cádiz and Almería) [11].This model included two population sizes and two migration rates.Model 2b: Two-population model that followed our groupings based on acorn protein profiles that indicated separation of western Sierra Morena (Huelva and Seville) and central and eastern populations (Córdoba, Cádiz, and Almería) [10].Model 2c: Two-population model inferred by STRUCTURE and Barrier results that suggested divergence of Cádiz from all other populations.Model 3: Three-population model.This model was based on biogeographic separation of the Sierra Morena populations (Huelva, Córdoba, and Seville) from Cordillera Penibética (Almería) and Cordillera Subbética (Cádiz) and included three population sizes and six migration rates.We used Bayesian inference (BI) to infer population θs (4N e µ, where N e is effective population size and µ is the mutation rate) and mutation-scaled migration rates (M = m/µ, where m is the migration rate per generation).The Brownian motion model was used to approximate a stepwise mutation model for microsatellites.Initial runs were performed to determine appropriate priors that were set for final runs as uniform (θ of 0.00-10.00,δ = 1.0 and M of 0.00-1000.00,δ = 100.00).Each run comprised a single long chain for which 1 × 10 6 genealogies were sampled with a burn-in of 1 × 10 4 and static heating with default temperatures.Initial parameter values were obtained from estimates of F ST .We estimated probabilities for each model by dividing its marginal likelihood by the sum of the marginal likelihoods of all models using the Bézier approximation.

Chloroplast DNA
We detected a total of 30 distinct haplotypes among the five regional populations (Figure 2).The total number of mutational steps among the haplotypes inferred from the shortest trees was 61, with most individuals from the same population clustering together with few mutational steps, with the exception of Córdoba.The greatest number of haplotypes [9] was found at Córdoba, whereas we detected only five haplotypes each at Seville and Huelva, three of which were shared in the two populations (Figure 2).Other than the sharing between Seville and Huelva, the remaining populations comprised unique sets of haplotypes.The shared haplotypes at Seville and Huelva suggest recent seed dispersal between the two populations that could be bidirectional.However, interestingly, two individuals (one each from Cádiz and Córdoba) clustered with the Huelva/Seville group (one mutational step for the Cádiz individual and five mutational steps for the Córdoba individual), suggesting unidirectional historical seed dispersal from Huelva/Seville.The number of inferred mutational steps among individuals within populations averaged 7 each at Cádiz (excluding the one individual that clustered with Huelva/Seville) and Almería, 14 in the combined Huelva/Seville cluster, and 21 at Córdoba.Seven mutated positions were inferred between the Huelva/Seville cluster and the closest Córdoba chlorotype (excluding the one individual that clustered with Huelva/Seville), whereas three and four mutated positions were inferred between Córdoba and the nearest haplotypes from Almería and Seville, respectively.
Global among-population differentiation considering the allele size ( cp R ST = 0.57; SE = 0.05) was significantly greater than among-population differentiation by allele identity ( cp G ST = 0.24; SE = 0.05), indicating that mutation had contributed significantly to population differentiation.

Nuclear DNA
Significant allelic dropout was detected by Microdrop for genotypes at Almería, Córdoba, and Seville.Of these, the Pearson correlation coefficient was significant across individuals at Córdoba and Seville and across loci at Almería.We imputed missing alleles for each population using the option in Microdrop in which both locus-and individual-specific factors are taken into account for imputed data sets.Inbreeding coefficients estimated for each population from the original data and the imputed data returned similar results; thus for the analyses reported here, we used the imputed dataset.
We detected a total of 101 alleles across populations at the nine microsatellite loci in Q. ilex.Allelic richness, on the basis of a minimum population size of 15, was least in Córdoba (5.5 alleles) and greatest in Huelva (7.3) (Table 2).Over all populations, we found negligible heterozygote excess (H O = 0.77; H E = 0.71); among populations, expected heterozygosity ranged from 0.68 to 0.72 but did not vary significantly (Table 2).

Past Demographic Change
BOTTLENECK detected significant values of heterozygote excess, consistent with a recent bottleneck using the Cornuet and Luikart [40] test, in the population at Córdoba (probability of the one-tailed test of heterozygote excess of <0.001).Inspection of the mode-shift distributions showed weak signs of bottleneck effects in all populations but most evidently at Cádiz and Córdoba (Figure S1).

Population Structure
The ad hoc statistic ∆K, which summarizes our results from 10 independent runs per K (from K = 1 to K = 5) in STRUCTURE, suggested that the microsatellite data on Q. ilex could be assigned to three clusters, as shown in Figure 3.Although the three clusters were admixed in individuals from the five populations, some population differentiation was clear.Cádiz and Almería were well differentiated from each other, whereas Huelva, Seville, and Córdoba were intermediately differentiated.

Past Demographic Change
BOTTLENECK detected significant values of heterozygote excess, consistent with a recent bottleneck using the Cornuet and Luikart [40] test, in the population at Córdoba (probability of the one-tailed test of heterozygote excess of <0.001).Inspection of the mode-shift distributions showed weak signs of bottleneck effects in all populations but most evidently at Cádiz and Córdoba (Figure S1).

Population Structure
The ad hoc statistic ΔK, which summarizes our results from 10 independent runs per K (from K = 1 to K = 5) in STRUCTURE, suggested that the microsatellite data on Q. ilex could be assigned to three clusters, as shown in Figure 3.Although the three clusters were admixed in individuals from the five populations, some population differentiation was clear.Cádiz and Almería were well differentiated from each other, whereas Huelva, Seville, and Córdoba were intermediately differentiated.For the DAPC analyses, cross-validation indicated that 30 principal component (PC) axes and 4 DA axes were optimal.In the space of the first two DA axes, Cádiz and Almería populations were divergent from one another and from the remaining populations that grouped together as partially overlapping clusters (Figure 4).For the DAPC analyses, cross-validation indicated that 30 principal component (PC) axes and 4 DA axes were optimal.In the space of the first two DA axes, Cádiz and Almería populations were divergent from one another and from the remaining populations that grouped together as partially overlapping clusters (Figure 4).A barrier of decreasing strength was detected by Barrier extending between Cádiz and Seville/Huelva, Cádiz and Córdoba, and Cádiz and Almería (Figure 1 and Figure S2).

Migration
Significant contemporary migration rates estimated by BayesAss were detected from Córdoba to Almería (95% CI: 0.19-0.33),Córdoba to Huelva (95% CI: 0.18-0.32),and Córdoba to Seville (95% CI: 0.06-0.18),for which mean rates are shown in bold in Table 3.All other population pairwise rates were not significantly different from zero, including the reciprocal rates into Córdoba, suggesting recent unidirectional gene flow.A barrier of decreasing strength was detected by Barrier extending between Cádiz and Seville/Huelva, Cádiz and Córdoba, and Cádiz and Almería (Figure 1 and Figure S2).

Migration
Significant contemporary migration rates estimated by BayesAss were detected from Córdoba to Almería (95% CI: 0.19-0.33),Córdoba to Huelva (95% CI: 0.18-0.32),and Córdoba to Seville (95% CI: 0.06-0.18),for which mean rates are shown in bold in Table 3.All other population pairwise rates were not significantly different from zero, including the reciprocal rates into Córdoba, suggesting recent unidirectional gene flow.Estimates of evolutionary among-population migration rates (number of migrants per generation N m ) were calculated from Migrate-N estimates of θ and M as N m = (θ × M)/4 (Table 3).The highest rate (N m = 105.5)was from Seville to Huelva, approximately 3 times greater than the next highest rates from Cádiz to Seville (N m = 37.2) and Seville to Almería (N m = 34.8).Migration rates were asymmetrical for several population pairs; notably, Córdoba was much more important as a source than as a sink population.
Comparisons of Bayes factors from Migrate-N for the different phylogeographic breaks supported a two-group model (Cádiz and (Almería, Huelva, Seville, Córdoba)).

Population Structure
Chloroplast DNA showed relatively strong spatial structure; only populations at Seville and Huelva shared haplotypes.Because chloroplast DNA is maternally inherited in oaks, it is dispersed only by seed; thus our results were consistent with the bulk of seed dispersal being limited locally or to adjacent habitat for Q. ilex [41][42][43][44][45], for other oak species [46][47][48], and for the acorn-bearing tanoak [49].We found no evidence of shared haplotypes among other population pairs, but we could infer at least two examples of probable long-distance dispersal that would have occurred for generations in the past.In these instances, one individual from each of Cádiz and Córdoba had haplotypes that clustered with the Huelva/Seville group, which would suggest past seed dispersal from Seville/Huelva to Córdoba and to Cádiz.Our data, using a suite of chloroplast microsatellite markers, detected fine-scale haplotype variation that had not been detected in earlier studies using less variable markers.For example, Lumaret et al. [23] reported only 25 haplotypes (on the basis of PCR-RFLP markers) throughout the entire range of Q. ilex and 3 within southeastern Spain, compared with the 30 that we detected in Andalusia.Our nuclear DNA showed a much weaker spatial structure (differentiation among populations) than the chloroplast DNA.This was expected because of the smaller effective population size of chloroplast DNA and because of wind pollination leading to much more widespread gene dispersal.Nevertheless, we found evidence of nuclear genetic structure among populations of Q. ilex in Andalusia.We used several approaches to investigate potential breaks.The STRUCTURE and DAPC analyses suggested weak genetic structure partitioned into three groups; Cádiz, Almería, and Sierra Morena populations (Seville, Huelva, and Córdoba).However, this partition was not supported by inferred contemporary (BayesAss) or evolutionary-scale migration rates (Migrate-N).The most consistent break suggested by the migration rate analyses was between Cádiz and the more northern populations of Huelva, Seville, and Córdoba.This suggests that the Guadalquivir Basin, where Q. ilex is absent, has presented a "recent" barrier to gene flow between the southern Cordillera Subbética and the Sierra Morena.The Guadalquivir Valley likely served as a refugium during the full glacial period, but early human deforestation of the valley may have broken continuity between the northern and southern populations [50,51].Divergence among plant populations of the southern Iberian mountains likely has diverse origins.Fortuna et al. [52], using a comparative network approach to study genetic variation among four woody species from southern Spain, found that each species had responded to ancient landscape processes rather than more recent events and in unique ways.For all species, divergence across the Guadalquivir Valley was significant, although Q. coccifera L. showed least population structure, likely related to unique life history traits.In some instances, the Gudalquivir divergence may be associated with closer affinities of Baetic populations to those of North Africa, as is the case for Erophaca baetica (L.) Boiss.[53] and the aster Hypochaeris radicata L. [54].
We found no consistent break within the Baetic Cordillera, despite the population in Almería being in the Mediterranean coastal range and the population of Cádiz being in the inner mountain range of the Cordillera Subbética.This could reflect connectivity through a low-medium mountain corridor until very recent times [55][56][57][58][59].Our Barrier analyses also found the strongest break to be between Cádiz and the Sierra Morena, with a second break that was non-significant between Cádiz and Almería.

Gene Flow
Fragmentation of populations in the recent past due to human activity may result in reduced levels of gene flow.The long-term effects of limited gene flow will be reduced genetic variability that is commonly assumed to lead to reduced population viability [60].BayesAss estimates recent migration rates (<5 generations) that would be in the time-frame of human activity in the Spanish dehesas [61,62].On the other hand, Migrate-N uses a coalescence approach that infers migration rates over the last 4N e generations [63], which, in Q. ilex, would likely cover the last tens of thousands of years and so would be affected mostly by events not associated with human activity.However, prehistoric human activities could have played a role in vegetation changes in the Guadalquivir Valley associated with fire, change of species, and scrubbing [64].To compare the two estimates of migration, we converted the mutation-scaled migration rate (M = m/µ) from Migrate-N by multiplying M by an estimated mutation rate µ of 5 × 10 −4 [65], following the approach described by Chiucchi and Gibbs [66].Our estimates of recent migration from BayesAss indicate low levels of migration among population pairs, except for unidirectional immigration from Córdoba to Almería, Huelva, and Seville.After transformation, evolutionary migration rates were about an order of magnitude greater than recent migration rates.Although this would be consistent with reduced gene flow among populations since human activity began in the dehesas, the results should be treated with caution.A change in the mutation rate has an important effect on the estimates of migration rates based on mutation-scaled migration rates.Assuming a mutation rate of 1 × 10 −4 would bring many of the evolutionary rates close to the contemporary rates.However, reduced contemporary migration may reflect, in part, the management regime of the dehesas that has been based on the existing forest and not on the establishment of new forests [67].
Contemporary patterns of migration among populations showed some similarities and some contrasts with patterns of evolutionary rates.Notably, we found a strong asymmetry between immigration rates into Córdoba and immigration into other populations from Córdoba for both contemporary and evolutionary migration.Flowering phenology could explain some of this asymmetry.In Q. ilex, most inflorescences are protandrous, and flowering phenology is earlier at warmer lower elevations [68].The Córdoba population has the most continental climate of all populations studied, and thus its winters are significantly colder [5,69].The Córdoba population would likely be more successful at pollinating receptive female flowers from lower elevations, such as Huelva and Seville, because of the overlapping of male and female phenology, whereas the reverse direction would lead to a greater dispersion between male and female flowering times.The most pronounced difference between contemporary and evolutionary migration rates was between Seville and Huelva.As expected, evolutionary rates were highest between this closest population pair, but contemporary rates were low in both directions.Without a doubt, these populations have been the most disturbed by human activity in historical times [70,71].

Association between Phenotypic Traits and Population Divergence
Does the phenotypic variability and the population grouping derived from this variability reflect the genetic structure of the Andalusian Q. ilex as determined by n-and cp-SSR?In a previous study on seed morphometry and fatty acid composition, 13 populations from 3 main geographical areas (southern, northeastern, and northwestern provenances) were grouped into 2 main clusters, with the Guadalquivir as a hypothetical frontier separating them [11].The first cluster corresponded to northern, mesic, low-altitude provenances, while the second corresponded to southern, xeric, high-altitude provenances.Protein profiles from seed and pollen, as determined by NIRS and one-dimensional gel electrophoresis, reinforced relationships among populations being related to geographic location and climatic zone [10].Geographically extreme provenances were clearly included in the northern (Sierra Morena) or southern (Baetic Cordillera) clusters, independently of the method employed [10][11][12].However, interestingly, depending on the methodology employed, some provenances were grouped in different clusters.This was the case for those provenances that were closest to the Guadalquivir Valley (southern margins of the Sierra Morena and northern margins of the Baetic range).It would therefore appear that the morphological and biochemical structuring are in fairly good agreement with the structuring obtained using our neutral genetic markers.This does not imply any causality, for which we would need to prove a genetic basis of the phenotypic traits and look for variation in functional genes that might explain any potential environmental adaptations.The phenotypic variation that we observed could be related to environmental variation, such as xeric stress, that is presumably only indirectly related to variability in microsatellite markers.Nevertheless, it is interesting that the genetic and phenotypic partitioning of populations showed some concordance that justifies more rigorous tests, such as comparative tests across aridity gradients in each of the major regions.
Our genetic data suggest moderately important variation that should be taken into consideration in any reforestation of Q. ilex in Andalusia.Reforestation north of the Guadalquivir Valley should be from seed sources of that region, particularly as there may be growth-related adaptive traits for the more mesic conditions that we have not yet tested.In the south, we consider that it would be prudent to use local sources for eastern and western reforestation as the genetic data suggests some divergence in this region.

Conclusions
We used genetic markers (chloroplast and nuclear microsatellites) to detect structure among populations of Q. ilex from the major biogeographic regions of Andalusia.Our results show surprisingly important levels of genetic structure, with chloroplast DNA indicating little gene exchange as a result of seed dispersal among populations.Our nuclear DNA analyses showed a weaker but significant genetic structure among populations, presumably as a result of more widespread pollen

Figure 1 .
Figure 1.Map of Andalusia (inset showing position of Andalusia in Spain) showing the distribution of holm oak forest and sampling locations (see Table 1 for population codes).Blue arrows indicate inferred seed dispersal events on basis of chloroplast DNA; green arrows show significant recent directional gene flow (see text).Black curved lines partition populations following Barrier analyses.
Figure 1.Map of Andalusia (inset showing position of Andalusia in Spain) showing the distribution of holm oak forest and sampling locations (see Table 1 for population codes).Blue arrows indicate inferred seed dispersal events on basis of chloroplast DNA; green arrows show significant recent directional gene flow (see text).Black curved lines partition populations following Barrier analyses.

Figure 1 .
Figure 1.Map of Andalusia (inset showing position of Andalusia in Spain) showing the distribution of holm oak forest and sampling locations (see Table 1 for population codes).Blue arrows indicate inferred seed dispersal events on basis of chloroplast DNA; green arrows show significant recent directional gene flow (see text).Black curved lines partition populations following Barrier analyses.

Forests 2018, 9 ,
x; doi: FOR PEER REVIEW www.mdpi.com/journal/forestsHuelva/Seville.The number of inferred mutational steps among individuals within populations averaged 7 each at Cádiz (excluding the one individual that clustered with Huelva/Seville) and Almería, 14 in the combined Huelva/Seville cluster, and 21 at Córdoba.Seven mutated positions were inferred between the Huelva/Seville cluster and the closest Córdoba chlorotype (excluding the one individual that clustered with Huelva/Seville), whereas three and four mutated positions were inferred between Córdoba and the nearest haplotypes from Almería and Seville, respectively.

Figure 2 .
Figure 2. Haplotype network based on 10 chloroplast microsatellite loci produced in NETWORK 5.0.0.1.Size of circles is proportional to number of individuals, and pie shapes are the proportion of individuals from each of the five populations.

Figure 2 .
Figure 2. Haplotype network based on 10 chloroplast microsatellite loci produced in NETWORK 5.0.0.1.Size of circles is proportional to number of individuals, and pie shapes are the proportion of individuals from each of the five populations.

Figure 3 .
Figure 3. Plot of STRUCTURE results with K = 3 of Q. ilex populations from Andalusia.Figure 3. Plot of STRUCTURE results with K = 3 of Q. ilex populations from Andalusia.

Figure 3 .
Figure 3. Plot of STRUCTURE results with K = 3 of Q. ilex populations from Andalusia.Figure 3. Plot of STRUCTURE results with K = 3 of Q. ilex populations from Andalusia.

Table 1
for population codes).Blue arrows indicate inferred seed dispersal events on basis of chloroplast DNA; green arrows show significant recent directional gene flow (see text).Black curved lines partition populations following Barrier analyses.

Table 1 .
Geographical and climatic data of the population of Quercus ilex subsp.ballota included in this study.

Table 2 .
Parameter estimates for chloroplast and nuclear microsatellite diversity in populations of Q. ilex from Andalusia, Spain.H e -expected heterozygosity; AR 15 -number of alleles on basis of a minimum population size of 15; F IS -inbreeding coefficient.Numbers in parentheses indicate standard errors of the estimates.

Table 3 .
Migration rates estimated as proportion of migrants in population from BayesAss (BA) and estimated as number of migrants per generation (N m ) from Migrate-N (M) with estimated migration rate in parentheses.Rates in bold significant at p = 0.05.Pairwise F ST shown in lower diagonal.