Are there morphological and life‐history traits under climate‐dependent differential selection in S Tunesian Diplotaxis harra (Forssk.) Boiss. (Brassicaceae) populations?

Abstract Adaptation of morphological, physiological, or life‐history traits of a plant species to heterogeneous habitats through the process of natural selection is a paramount process in evolutionary biology. We have used a population genomic approach to disentangle selection‐based and demography‐based variation in morphological and life‐history traits in the crucifer Diplotaxis harra (Forssk.) Boiss. (Brassicaceae) encountered in populations along aridity gradients in S Tunisia. We have genotyped 182 individuals from 12 populations of the species ranging from coastal to semidesert habitats using amplified fragment length polymorphism (AFLP) fingerprinting and assessed a range of morphological and life‐history traits from their progeny cultivated under common‐garden conditions. Application of three different statistical approaches for searching AFLP loci under selection allowed us to characterize candidate loci, for which their association with the traits assessed was tested for statistical significance and correlation with climate data. As a key result of this study, we find that only the shape of cauline leaves seems to be under differential selection along the aridity gradient in S Tunisian populations of Diplotaxis harra, while for all other traits studied neutral biogeographical and/or random factors could not be excluded as explanation for the variation observed. The counter‐intuitive finding that plants from populations with more arid habitats produce broader leaves under optimal conditions of cultivation than those from more mesic habitats is interpreted as being ascribable to selection for a higher plasticity in this trait under more unpredictable semidesert conditions compared to the more predictable ones in coastal habitats.

population-level variation and its fitness consequences are assessed (Geber & Griffen, 2003). Local adaptation in functional traits through natural selection could be either studied by trait-fitness correlations (selection differentials and selection gradients) in natural populations (Endler, 1986;Kingsolver et al., 2001), by field experiments such as common-garden or reciprocal transplant experiments under assessment of trait and fitness, or by population genetic/genomic approaches (Savolainen, Lascoux, & Merilä, 2013). While in the two former method groups, the phenotypes are the starting point for selection studies and an assessment of fitness is necessary, the latter is less straightforward and starts with variation at the DNA level as a footprint of local adaptation, infers loci (so-called outlier or candidate loci) that are under selection, and tests for statistically significant correlations of those loci with functional traits (Storz, 2005;Vasemägi & Primmer, 2005).
Studies on local adaptation are especially informative when clinal variation in functional traits is studied along environmental gradients, and alternative explanations for phenotypic clines (e.g., population demography, phylogeographical patterns) are excluded (Keller et al., 2009;Keller & Taylor, 2008;Savolainen et al., 2013). Among the different clinally varying environmental variables, aridity has been in the focus of several studies on local adaptation in the last years, mainly as a consequence of questions raised in conjunction with expected climatic changes and their impact on the vulnerability of plant species (e.g., Rehfeldt, Wykoff, & Ying, 2001). In a common-garden experiment carried out in the Botanical Garden of Tel Aviv University (Israel), Petrů et al. (2006) found support for the expectation that desert populations have a much higher potential to respond plastically to differences in rainfall conditions. When irrigated, individuals of Biscutella didyma L.
(Brassicaceae) from arid environments produced almost three times as many seeds as plants from the mesic Mediterranean. Therefore, while the former was considered following a strategy of maximizing reproduction during good years, a significantly higher investment of the latter populations into vegetative growth was interpreted as adaptation to aboveground competition in more stable and environmentally predictable habitats. A recent population genomic study of Geropogon hybridus (L.) Sch.Bip. (Compositae) along a steep precipitation gradient in C Israel revealed natural selection at 11 of 123 polymorphic loci studied that showed signs of putative precipitation-related adaptation (Müller et al., 2017).
The coastal region of SE Tunisia exhibits steep gradients in precipitation, ranging from areas with mean annual precipitation values larger than 200 mm on the island of Djerba, over an approximately 20-km broad strip parallel to the Mediterranean coast with precipitation values between 150 and 200 mm, to semidesert areas west of the more inland Matmata Mts with precipitation values of below 100 mm (Frankenberg & Klaus, 1987). Due to their elevation reaching more than 500 m a.s.l., the Matmata Mts themselves are again favored with mean annual precipitation values reaching those of lowland and coastal Djerba. As a consequence of these steep gradients, changes among vegetation types are pronounced and range between Mediterranean coastal vegetation (Djerba), over different steppe types in the coastal strip, to plant associations substituting Juniperus phoenicea woodland (Matmata Mts) and the semidesertic Anthyllis sericea steppe type at the eastern boundary of the Sahara desert (Frankenberg & Klaus, 1987). This constellation provides a suitable setting for studying climate-dependent differential selection of morphological and life-history traits along an aridity gradient.
The crucifer Diplotaxis harra (Forssk.) Boiss. (Brassicaceae) is a southern Mediterranean annual to short-lived perennial herb with its distributional range also expanding to S Europe (Spain, Sicily) and SW Asia (Israel, Lebanon, Syria). Morphologically, it is easily characterized by its pendent siliquae with seeds in two rows. In North Africa, the species is widely distributed both in Mediterranean and in more arid, often sandy or gypseous areas with elevated salinity (Tlig, Gorai, & Neffati, 2008) and exhibits a broad range of ecological life-cycle types with different reproductive strategies, ranging from ephemeral over modular monocarpic to coppiced polycarpic behaviors in dependence of the precipitation regime (Hegazy, 2001). In SE Tunisia, Diplotaxis harra is found along the above described, steep precipitation gradients and is therefore a suitable object for studying local adaptation in morphological and life-history traits by means of a combined commongarden and population genomics approach. Following the strategy of Herrera and Bazaga (2008), who pinpointed floral traits being under adaptive divergence in the Spanish hawk moth-pollinated violet Viola cazorlensis Gand. (Violaceae), this study aims at the detection of amplified fragment length polymorphism (AFLP) loci under selection and asks whether these are correlated with phenotypic traits of interest, assuming that environment-dependent selection in these morphological and life-history traits has led to a selection signal in some of the AFLP loci via genetic hitchhiking.

| General strategy
Following the argumentation of Taylor (2008) andKeller et al. (2009), for the study of natural selection of phenotypic traits in wild populations, it is important to discriminate effects due to selection from changes caused by prior evolutionary history and chance events. Adopting the analytical sequence suggested by Herrera and Bazaga (2008) in their study of natural selection in flower characteristics of Viola cazorlensis Gand. (Violaceae), we first performed AFLP fingerprinting in populations of Diplotaxis harra in SE Tunisia in order to discriminate putative loci subject to selection from neutral ones. This was done with two methods (MCHEZA, Antao & Beaumont, 2011;BayeScan, Foll & Gaggiotti, 2008) that use locus-wise measures of population differentiation (F ST values) to find loci with differentiation values significantly larger than average as candidates for as being under divergent selection ("outlier loci"). Additionally, a third approach was used (as implemented in Samβada; Joost et al., 2007;Stucki et al., 2017) that searches for AFLP loci under selection by testing the strength of relationship between the presence or absence of an allele in a single-individual genotype and environmental variables using logistic regression models. In a second step, we tested for geographical signals in the two nonoverlapping subsets of loci ("selected" vs. "neutral" ones) in order to detect isolation-by-distance (IBD) caused by demography and neutral history of populations (phylogeography) in the background of either natural selection along geographical gradients ("selected loci" with a significant IBD signal) or natural selection without a geographical motivation ("selected loci" without a significant IBD signal). Isolation-by-distance (IBD) was tested for using Mantel's test of correlation between a matrix containing geographical distances and a matrix with pairwise F ST values among populations. Finally, we searched for correlations between life-history and morphological traits and either of the two nonoverlapping subsets of loci ("selected" vs. "neutral" ones) by means of a Generalized Analysis of Molecular Variance (GAMOVA; Nievergelt, Libiger, & Schork, 2007). Here, we looked for traits being simultaneously correlated with the set of loci under selection and not correlated with the set of neutral loci because only in this constellation (scenarios 1.1 and 1.2 in Table 1) neutral variation as cause for phenotypic divergence among populations could be excluded and clinal or local selection could be assumed with confidence.

| Plant material and cultivation
The 12 populations of Diplotaxis harra surveyed in this study were sampled in Tunisia in March 2009 covering a climatic gradient between coastal Mediterranean conditions with relatively high precipitation and steppe vegetation in the east (Djerba) and semidesert areas with pronounced aridity in the west (Tozeur; see Table S1). In total, 182 individuals were sampled with 12-19 individuals per population. For each individual, leaf material was collected and dried with silica-gel, and ripe fruits were stored in separate paper bags for the cultivation of half-sib progeny in the Botanical Garden of Regensburg University.
Seeds from the 182 plant individuals sampled were sown on 5 March 2015 into a 1:1:2 mixture of sand, compost, and Einheitserde ® Classic (Sinntal-Altengronau, Germany). After germination, five of the most vital half-sib seedlings of each pot were selected and pricked out into larger pots with a 1:1 mixture of sand and Einheitserde ® Classic between 16 March and 25 March 2015 for further cultivation. On 8 April 2015, the five pricked-out seedlings were reduced to a single one by removing the other four. While these first steps of cultivation were carried out in a glasshouse under a constant temperature regime (25°C), the further cultivation was carried out still under glass but under natural temperature conditions to avoid etiolation and rampant growth. During all steps of the experiment, we kept plants under identical light and water regimes, randomizing the pot positions every 3rd day. Due to germination failure in 16 seed samples from different populations, we propagated two or three half-sib seedlings from 19 of single mother plants, respectively, in order to keep sample sizes of populations in balance. Therefore, in total, 187 seedlings were raised and scored for the life-history traits and morphological characters described in the following.

| Assessment of life-history and morphological traits
Life-history traits and morphological characters scored for each plant under cultivation are summarized in Table S2. In terms of lifehistory traits, we recorded the number of days between sowing (March 25) and the opening of the first flower (variable V01) and the time between the starting of flowering and the formation of the first ripe fruits (V03). Plant architecture was described by scoring T A B L E 1 Hypothetical outcomes from GAMOVA (Generalized Analysis of Molecular Variance; Nievergelt et al., 2007) tests addressing the correlation of quantitative traits with selected vs. neutral loci and the results of Mantel (Mantel, 1967)  each plant for total height of the main shoot (V09) and the number of side shoots produced (V10). Additionally, the total number of leaves (V07), the number of flowers along the main shoot (V04), and the number of seeds per fruit (mean number of three siliques surveyed; V06) were assessed. For detailed morphometric analyses of leaves and flowers (described below), the sixth leaf produced along the main shoot was removed and dry-pressed in a herbarium press, while a representative flower of each plant was removed and photographed keeping constant its dorsiventral orientation and adding a millimeter scale for calibration. For the assessment of leaf hairiness, a 5-mm long stretch in the middle of each of the two lateral leaf margins was counted out for hairs, and the two values were averaged (V08).

| Geometrical morphometrics of leaves and flowers
Size and shape of dried leaves were assessed after scanning with a scale for calibration. We used ImageJ v1.49 (Schneider, Rasband, & Eliceiri, 2012) for the measurement of the absolute area of the leaf lamina in cm 2 (V02). For the comparison of leaf shapes, scanned-in leaf silhouettes were imported into the outline analysis package Momocs (Bonhomme, Picq, Gaucherel, & Claude, 2014) in R v3.2.1 (R Development Core Team 2015), and an elliptical Fourier decomposition (Giardina & Kuhl, 1977;Kuhl & Giardina, 1982) was conducted.
Following suggestions of Daegling and Jungers (2000), a preliminary analysis of all scanned leafs was carried out showing that the extraction of 15 harmonics was sufficient to describe 99% of the outline of even deeply dissected leaves. Subsequently, extracted Fourier harmonics were subjected to a principal component analysis (PCA) to gain PCA coordinates for each OTU as further three variables (V15-V17) for the statistical analyses.
In contrast to the Fourier approach for outline analysis of leaves, we followed a landmark/semilandmark approach for flower shapes as done in a comparable study in Erysimum (Brassicaceae) by Gómez, Perfectti, and Camacho (2006 μl Tris-EDTA buffer. Amplified fragment length polymorphism analyses generally adhered to the protocol described by Vos et al. (1995) and were carried out as described in Meister, Hubaishan, Kilian, and Oberprieler (2006)

| Analyses of population structure
We analyzed the binary matrix from AFLP fingerprinting by a principal coordinates analysis (PCoA) based on pairwise Bray-Curtis distances among the 181 operational taxonomic units (OTUs) using the program MVSP v3.1 (Kovach, 1999). We further used a Bayesian clustering approach as implemented in the program Structure v2.3.3 (Falush, Stephens, & Pritchard, 2007;Pritchard, 2010;Pritchard, Stephens, & Donnelly, 2000). The optimal number of genetic clusters (k) was calculated according to the method of Evanno, Regnaut, and Goudet (2005) using k values ranging from 1 to 9 with ten repetitions for each k. The data sets ran using the admixture model with 150,000 generations after 150,000 burn-in generations.

| Search for loci under selection
We have used three different approaches to search for AFLP markers that are either under direct selection themselves or closely linked to loci under selection. As the first two frequentist F ST outlier detection methods (i.e., MCHEZA, BayeScan) make the assumption of populations being in Hardy-Weinberg equilibrium and we were unaware about the realization of this prerequisite in the populations under study, we applied a third spatial analysis method (Samβada) that does not use a population genetic approach.

The program MCHEZA (Antao & Beaumont, 2011) is based upon the Dfdist kernel which searches for loci under selection
based on the principle that genetic differentiation among populations is expected to be higher for loci under divergent selection than for the rest of the genome (Beaumont & Balding, 2004;Beaumont & Nichols, 1996). The program calculates locus-wise allele frequencies and F ST values using the Bayesian method of Zhivotovsky (1999) and uses computer simulations to produce F ST values for modeled AFLP loci under neutral conditions. As suggested by the software developers, we used the "neutral mean F ST " and the "force mean F ST " options to calculate a "trimmed mean F ST ," which is then used as an estimate of the average "neutral" F ST uninfluenced by outlier loci.
Loci with a noticeable high F ST value were considered as being under divergent selection. The analysis was performed using 1,000,000 simulations, a 0.95 confidence interval (CI), and a false discovery rate (FDR) of 0.05 to correct for multiple testing.
For all other parameters, we used the default values.

2.
The program BayeScan v2.1 (Foll & Gaggiotti, 2008)  BayeScan uses a regression approach and a hierarchical Bayesian procedure to simultaneously estimate F ST values for every locus in each population. For this analysis, we started with 50 pilot runs of 5,000 iterations each and finally ran a total number of 250,000 iterations with a burn-in of 50,000 and a thinning interval of 20, resulting in a sample size of 10,000 iterations. All other parameters were set to default. As above, we used a FDR of 0.05 for outlier detection.

3.
A third method to detect adaptive loci is the spatial analysis method as implemented in the program Samβada (Joost et al., 2007;Stucki et al., 2017). In contrast to the two above-mentioned procedures, it does not analyze AFLP data in a population genetic framework but is based on logistic regression models testing the strength of relationships between the presence or absence of an allele/band in a single-individual genotype and environmental variables. As explained by Joost et al. (2007), this procedure has some advantages over the population genetic approaches by not being dependent on any assumption of inbreeding coefficients, therefore being applicable to sampling designs with only single individuals per location, and immediately returning direct relationships between environmental variables and candidate loci under selection coupled with these variables. Two statistical tests [a log-likelihood ratio test and the Wald test (Wald, 1943)] are implemented to evaluate the significance of the correlation between the presence/absence of an AFLP band at a locus and an environmental variable. We ran univariate models including the 20 climatic variables extracted for the 12 sampled populations (Table S3) and used a significance threshold of 0.05, which was Bonferroni corrected for multiple comparisons.
Aiming at a reduction in complexity, we additionally submitted the

| Correlation of selected and neutral loci with geography
After assignment of loci to one of the two nonoverlapping subsets of loci ("selected" vs. "neutral" ones) according to the results of the MCHEZA, BayeScan, and Samβada analyses, respectively, we performed tests for isolation-by-distance (IBD) in each of the resulting locus subsets by Mantel tests (Mantel, 1967)  repetitions.

| Variation in life-history and morphological traits
Boxplots of all life-history characters (i.e.,   Figure S4), in which k = 2 was found as the optimal cluster number, and population F I G U R E 1 Results of a principal component analysis (PCA) based on the extraction of 15 Fourier harmonics describing the outline shape of the sixth cauline leaf. For each of the three extracted principal components (PC 1 to PC 3) explaining more than 5% of the total variance, the variation of leaf shape is depicted with five representative leaf silhouettes along the axis concerned. While PC 1 explains 54.7% of the total variance, PC 2 with 12.4% and PC 3 with 9.5% are considerably less important than the first axis

| Analysis of bioclimatological data
The mentioned bipartition of sampled populations of this study found in the AFLP fingerprinting is mirrored in the climatological patterns.
When bioclimatic variables (bio01-bio20) for the 12 populations are subjected to a principal component analysis (PCA; Figure S5, Table   S4), the first axis PC1 (accounting for 55.7% of the total variance in the data set vs. only 26.9% on PC2 and 11.5% on PC3) clearly separates populations 11 (Tozeur) and 12 (Selja)  that these populations form the extremes in a gradient between relatively high precipitation, low temperature, and low seasonality/continentality at the former locations and relatively low precipitation, high temperature, and high seasonality/continentality at the latter ones.

| Loci under selection
When all 12 populations under study were considered in the search for AFLP loci under selection, MCHEZA simulations resulted in a trimmed mean of F ST = 0.0249 over all 732 loci (Table 2). Seventeen (2.3% of the total) fell outside the 95% and 14 (1.9% of the total) outside the 99.5% confidence interval constructed from the simulated distribution ( Figure S7). While all of these ( a highly significant correlation detectable between the absence/presence of a band and a single or more environmental (bioclimatic) factors. Twelve of these loci were from the set of 16 outlier loci found with the above-mentioned F ST -based methods.
When the search for loci under selection was based on the ten genetically more homogenous populations of SE Tunisia (except the isolated populations 11 from Tozeur and 12 from Selja), the results were much less bulky compared to the complete data set (Table 3): While MCHEZA found 11 loci falling outside a 99% confidence interval for neutral loci, five of these could be also corroborated as being under selection in the BayeScan analysis.
Of these, three loci (L039, L040, L459) were also found with band presence/absence being significantly correlated with at least one bioclimatic factor.

| Relationship between genetic and phenotypic divergence
Tables 4 and 5 summarize the results of the Generalized Analysis of Molecular Variance (GAMOVA) for the complete (12 populations) and T A B L E 2 Characterization of 26 outlier loci found as being under selection by one (italics), two (simple), or three (bold) of the three methods implemented in MCHEZA, BayeScan, and Samβada, respectively, in the complete data set (pop01-pop12 the reduced data set (10 populations), respectively. In the complete data set (Table 4) Following the reasoning behind our search strategy for lifehistory and morphological traits under selection (Table 1), V15 represents the only trait, for which neutral variation could be ruled out significantly. Regressing of population means of V15 on the bioclimatic variables with a significant contribution to the selection regime of the reduced data set (i.e., bio05, bio12, bio13, bio18, PC3; Table S3), however, revealed quite spurious relationships (Figure 4, left). The situation improved considerably after removing population 1 (Djerba) and revealed marked trends for the relationships of V15 (leaf shape) with bio05 (max. temperature of warmest month), bio12 (annual precipitation), bio13 (precipitation in wettest month), and bio18 (precipitation in warmest quarter; Figure 4, right). Deducing from these trends, narrow leaves were found, therefore, in plants from populations with lower temperatures in the warmest month (bio05) and higher annual and seasonal precipitation values (bio12, bio13, bio18).

| Correlation of selected and neutral loci with geography
Tables 4 Table S2; V11-V14: principal component coordinates PC 1 to PC 4 of flower morphometrics, see Figure 2; V15-V17: principal component coordinates PC 1 to PC 3 of leaf morphometrics, see Figure 1), significant correlations with outlier and neutral locus sets resulting from one (Samβada), two (MCHEZA, BayeScan), and three (MCHEZA ∩ BayeScan ∩ Samβada) analyses for loci under selection are given in bold. Additionally, for each locus set, the results of a Mantel (Mantel, 1967) test for correlation between pairwise geographical distances and pairwise F ST values among populations are given

Pseudo-F p-Value
Pseudo-F  Seldja" while many of her indications for the former subspecies (e.g., "Gafsa, Matmatas") rather correspond to localities sampled for the present study under population numbers 7, 8, and 10, it may appear self-evident that the genetic (and morphological) bipartition observed might be caused by sampling of the two subspecies and might indicate rather a taxonomic-biogeographical structure than a selectiondependent and taxonomy-independent one. There are, however, considerable arguments against this interpretation: (1) in contrast to the above-mentioned two sources, Maire (1965) in his treatment of the genus in the Flore d′Afrique du Nord only indicates the first subspecies, D. harra subsp. harra, for Tunisia, but not the second one (only found in Morocco and Algeria); (2) according to diagnostic differences between the two subspecies given by Fennane (1999), Maire (1965, and Pottier-Alapetite (1979) with leaves of subsp. harra being T A B L E 5 Results of the GAMOVA (Generalized Analysis of Molecular Variance; Nievergelt et al., 2007) tests conducted in the resulted data set (pop01-pop10) of the present study. For each of the 17 variables (V01-V10: morphological and life-history traits, see Table S2; V11-V14: principal component coordinates PC 1 to PC 4 of flower morphometrics, see Figure 2; V15-V17: principal component coordinates PC 1 to PC 3 of leaf morphometrics, see Figure 1), significant correlations with outlier and neutral locus sets resulting from one (Samβada), two (MCHEZA, BayeScan), and three (MCHEZA ∩ BayeScan ∩ Samβada) analyses for loci under selection are given in bold. Additionally, for each locus set, the results of a Mantel (Mantel, 1967)  F I G U R E 4 Linear regression of trait V15 (leaf shape of the sixth cauline leaf described by PC 1 of the morphometric analysis; see Figure 2) on four individual bioclimatic variables and the third axis (PC 3) of a principal component analysis (PCA) of all 20 bioclimatological variables (Table  S3). In the left column, all ten populations (pop1-pop10) of the reduced data set are used, while in the right column, pop1 (Djerba) has been omitted, leading to a higher coefficient of determination r 2 in the case of the four individual bioclimatic variables bio05, bio12, bio13, and bio18 crassifolia is described as being strictly perennial and subshrubby, a feature definitively not exhibited by the plants of the present study, which died at the end of the fruiting period even under optimal glasshouse/garden conditions. This is additionally corroborated by observations of Diplotaxis harra (subsp. harra) in Egypt done by Hegazy (2001) who describes that a part of populations of this species exhibits coppiced and perennial behavior. Therefore, we are quite confident that the results of our selection study are not invalidated by a background taxonomical structure and infraspecific biogeographical or neutral demographic history (phylogeographical) patterns. Nevertheless, in order to safeguard our results against any possible taxonomy-prone distortions (and possible biases caused by the artificial geographical sampling gap between populations 1-10 and populations 11-12), we did all of our analyses both in the full (populations 1-12) and in a reduced data set (populations 1-10).

| Detecting selection in SE Tunisian Diplotaxis harra populations
In both data sets, the complete one with populations 1-12 and the reduced one excluding populations 11 and 12, we found a good, yet not complete correspondence between the two compared to studies carried out at both regional or even continental scales (Manel et al., 2010). Another explanation, however, may be that our conservative approach for outlier detection may have led to an inflated type II error and an underestimation of the number of loci under selection in our study compared to other ones. However, following the argumentation of Henry and Russello (2013), we consider the applied cross-validation procedure a valid strategy for increasing confidence in identified outliers.
We found that inclusion of the two geographically, climatically, and genetically deviating populations from the west of Chott el Djerid was tremendously influencing the number of outlier loci both detected by the two F ST -based and the environment-allele association approaches. The Samβada results (Table 2) demonstrate that divergent selection at 19 of 26 outlier loci is caused by climatological variables, either by individual ones (e.g., bio15, precipitation seasonality), but mostly by conglomerates of many, covarying ones and their surrogates PC1 and PC3 of the principal component analysis (PCA). Interestingly, almost all climatological variables with high factor loadings on PC2 (i.e., bio01, bio03, bio06, bio11, bio14, and bio20; see Table S4), which is the principal component axis that  (Table 4), we cannot preclude completely that this larger-scale pattern may be influenced by demography and neutral history of populations (phylogeography) in the background of natural selection.
When the two deviating semidesert populations are excluded, environmental gradients get much shorter but exhibit a more continuous pattern of variation among the ten remaining habitats surveyed (as demonstrated by the PCA in Figure S6), allowing for a possibly less phylogeographically biased analysis of molecular and trait selection. This expectation is supported by Lotterhos & Whitlock (2015) who have hypothesized that local adaptation may be most informatively studied by outlier detection methods when ecologically divergent populations sampled are geographically adjacent (and therefore genetically similar for neutral genes). As expected, the search for candidate loci resulted in a smaller number of loci under selection, five being jointly detected by the two F ST -based methods and only three of these corroborated by the environment-allele association approach (Table 3). The climate variables responsible for this small-scale selection were identified as being connected to precipitation (bio12: annual precipitation; bio13: precipitation in wettest month; bio18: precipitation in the warmest quarter), temperature (bio05: maximum temperature in the warmest month), and the climatic surrogate variable of PC3 (being mainly under the influence of latitude, longitude, and bio15: precipitation seasonality); all of them fitting the steep gradient between Mediterranean and steppe vegetation described by Frankenberg and Klaus (1987) for the study area.

| Traits under climate-dependent selection in SE Tunisian Diplotaxis harra populations
Irrespective of the composition of the set of outlier loci found with the three detection methods applied, the GAMOVA analyses of the complete data set revealed statistically significant correlations between divergent selection patterns and two morphological traits (Table 4): V02 (accounting for the absolute size of the sixth leaf produced along the main shoot) and V15 (length-width ratio of the lamina of the sixth leaf). However, as both traits (along with three others, V04: number of flowers along main shoot; V07: number of leaves along main axis; V17: pointedness of lamina of sixth leaf) also show strong correlations with the set of neutrally evolving loci, phenotypic divergence among populations could be caused by both selection and neutral processes, and the latter could not be excluded with confidence as explanation (see Table 1). Extremely strong isolation-by-distance (IBD) signals for both outlier and neutral loci sets corroborate the interpretation of morphological traits showing variation rather caused by neutral biogeographical and/or random factors than by environmentally mediated divergent selection.
The picture is slightly different in the small-scale analysis with the reduced data set of populations 1-10 (Table 5). In the outlier loci set found with the two F ST -based methods, only V07 (number of leaves) was detected as being under divergent selection; however, as in the complete data set, this morphological character showed an equally strong correlation with the set of neutrally evolving markers and is interpreted here again as possibly caused by random or population history factors. When the number of outlier loci is reduced to three by additionally accounting for the results of the Samβada analysis; however, V15 (describing the shape of leaves in the Fourier outline analysis) displays a significant correlation with the outlier loci set, which is not paralleled by a correlation with the neutral loci set. As the outlier loci also exhibit a significant IBD signal, this indicates that leaf shape and especially the width-length ratio determining the first principal component of the Fourier analysis (PC1: V15) is under a climate-mediated divergent selection along a geographical cline (see Table 1). Following studies addressing the genetic background of leaf-length/width ratio (leaf index) in the crucifer Arabidopsis thaliana (Horiguchi, Fujikura, Ishikawa, & Tsukaya, 2006;Tsukaya, 2006), the leaf index is controlled by four regulating systems being the combination of genes influencing the polar cell elongation (cell shape) or cell proliferation (cell number) in either the leaf-length or the leaf-width direction. Therefore, one could legitimately assume that this important functional trait is under oligogenic control also in the related crucifer Diplotaxis harra and that the three outlier loci (L039, L040, and L459) found could be linked with at least one of these genes. Methods described by Paris et al. (2012) could be used in the future to gain nucleotide sequences of these AFLP loci of interest and test their identity with length/width ratio genes of Arabidopsis.
The correlational trends in leaf shape as a morphological trait under environmental selection, with broader leaves being observed in plants from semidesert populations with higher temperatures in the warmest month (bio05) and lower annual and seasonal precipitation values (bio12, bio13, bio18), appears being counterintuitive because xeromorphic leaves are usually smaller and narrower than leaves from more humid environments. This tendency has been observed, for example, in leaves of summer ecotypes vs. leaves of winter ecotypes in the Sicilian Diplotaxis erucoides (Schleser, Bernhardt, & Hurka, 1989). Additionally, among other traits, length-width ratio of leaves was found exhibiting considerable plasticity in Piriqueta caroliniana Urban (Turneraceae), with narrower leaves being produced under more arid conditions in successive years (Picotte, Rhode, & Cruzan, 2009) or even in successive seasons of the same year (Picotte et al. 2007). However, as we lack measurements from natural habitats and our experiment has not been carried out under natural conditions in common gardens with Mediterranean and/or semidesert climate but under glasshouse/ garden conditions with optimal supplies of water, light, and nutrients, the observed trends rather demonstrate that plants with provenances from semi-arid habitats either have broader leaves or are capable to produce broader leaves than those from more mesic conditions when confronted with optimal growth conditions and that this trait or the plasticity of this trait is selected for. As described by Nicotra et al. (2011: 542), the functional significance of leaf shape may lie in its influence on leaf thermal regulation and its hydraulic properties, the latter "likely to be of greater importance." While it is easily conceivable that more unpredictable precipitation at semi-arid inland habitats may select for a higher plasticity in leaf shape than the more constant conditions along the Mediterranean coast, allowing plants to adjust to highly temporally variable aspects of transpirational demand (Sultan, 2000), our present observations could only be the starting point for more comprehensive, hypothesis-driven experiments addressing the evolutionary significance of this functional trait and its plasticity. Common-garden cultivation under different watering regimes (genotype-by-environment interaction) should be the first experiments in that direction (Des Marais, Hernandez, & Juenger, 2013), followed by artificial selection experiments and experiments addressing the question whether leaf-shape plasticity is adaptive (Scheiner, 1993).

ACKNOWLEDGMENTS
The help of Mrs Claudia Gstöttl (Regensburg) and Dr Robert Vogt (Berlin) with collection of plant material (leaves, seed) of Diplotaxis harra during the 2009 excursion to Tunisia is gratefully acknowledged.
For the cultivation of plants in the Botanical Garden of Regensburg University, we would like to express our sincerest thanks to Volker Debus and his team.

CONFLICT OF INTEREST
None declared.