Centaurea hyssopifolia Vahl. (Compositae) is a dominant species of the dwarf shrub community on Iberian gypsum soils, restricted to gypsum massive soils of central Spain (Luzuriaga et al., 2006). The species is a medium-lived, small chamaephyte (20–60 cm) that exclusively grows on gypsum outcrops, i.e., it is a genuine gypsophile. It occurs on the upper piedmont of gypsum slopes, where a sparse shrub community of gypsophiles establishes. Long-term field observations suggest that individuals of C. hyssopifolia have an average longevity of 4–10 years. Flowering spans from May to July, with inflorescences (capitula) presenting purplish disk florets and a few white ray florets. The species is described as a mainly outcrosser with partial self-compatibility (Luzuriaga et al., 2006), which concurs with the mating system of other close congeners, and has limited dispersal mechanisms, with achenes with almost no pappus. Preliminary field data revealed that individually bagged inflorescences produced no viable seeds (Luzuriaga et al., unpublished). Field pollinator censuses showed that the study species is visited by a taxonomically diverse group of pollinators, including Bombyliidae species, several Apidae species as well as Lepidoptera and Coleoptera (unpublished data).

Populations were sampled in a fragmented gypsum landscape in the Tajo River Basin, near Chinchón, central Spain. Long-term, intensive agricultural activities have progressively fragmented the landscape, creating a mosaic of patches of natural vegetation remnants interspersed in a matrix of dry herbaceous croplands (cereal), dry arboreal croplands (olive trees), vineyards and Pinus halepensis plantations. Climate is semiarid Mediterranean with mean annual rainfall of 422 mm and mean temperature of 13.8°C. The study fragments laid within a 5 km × 6 km rectangular area (600–700 m a.s.l.), so climatic conditions across sites were very similar. We selected 10 populations from fragments that differed markedly in their size and degree of connectivity (see full details of fragment selection in Matesanz et al., 2015). We used high-resolution orthophotos of the area to identify and digitize remnant patches of gypsum vegetation (≈300 fragments). We measured the area of each fragment and calculated the minimum distance among surrounding fragments. To quantify the connectivity of each fragment (i.e., inverse of isolation), we used an index that accounts for the number of surrounding fragments weighed by their distance to the target fragment and their size (Tremlova and Muenzbergova, 2007). Specifically, connectivity was expressed as:

where C_i (unitless) is the connectivity of fragment i, n is the total number of natural fragments around the target fragment that are included within a 500 m radius circle, A_k is the area of fragment k, and d_i is the minimum distance between fragment i and k (see Supplementary Table S1 for distance among selected fragments). We used a conservative radius of 500 m because pollen movement among fragments located at larger distances is likely to be minimal, since generalist pollinators can forage pollen from many sources within a single fragment (Aizen et al., 2002; Fontaine et al., 2008). This index was selected because it takes into account not only the number of fragments of natural habitat surrounding the target fragment but also the among-fragment distance and their size, providing a complete measure of landscape connectivity. Based on these preliminary data, we surveyed 100 fragments in the field and finally selected 10 fragments of contrasting size and connectivity based on the following criteria: (i) the study species was abundant, with at least 150 individuals in each fragment, to evaluate the effects of fragment size per se, irrespective of population size and (ii) there was a well-developed biological soil crust, since its presence is indicative of unmanaged, undisturbed gypsum habitats. We selected fragments whose age (earlier date at which the fragment was identified in historical aerial photographs; Table 1) was much higher than the generation time of the species. The selected fragments create both a size and a connectivity gradient. Fragment size ranged from 1241 to 138849 m² (100-fold difference), and fragment connectivity ranged from 1.83 to 6.97 (Table 1). In order to test the independent effects of both components of habitat fragmentation –size and connectivity– as well as their interaction, we ensured that the gradients were not correlated (non-significant correlation between fragment size and connectivity R = 0.12, P = 0.76).

At the peak of the reproductive season in 2012, we randomly selected 12 reproductive individuals located within a 20 m × 20 m plot that was established in the south or south–east slope of each fragment, to minimize microclimatic differences among plots (Table 1). Plants were located ≥1 m from each other and were chosen to represent a random sample of the plant size distribution in each fragment. We collected 50–70 mature capitula from each mother plant (hereafter family). Capitula were dissected and viable seeds were cleaned and stored. In total, we sampled 120 families from 10 fragments. Since the sampled maternal families were open-pollinated, plants from the same family are considered half-siblings (Falconer and Mackay, 1996). Fresh leaves were also collected from 20 plants per fragment (including the maternal plants), air dried and stored in paper bags.

Neutral genetic variation was estimated using eight microsatellite markers transferred from other species of the genus Centaurea (Fréville et al., 2000; Marrs et al., 2006; Austin et al., 2011; see Supplementary Table S2). DNA was extracted from 60 mg of leaf tissue using the SpeedTools Plant DNA Extraction kit (Biotools, Madrid, Spain), following the manufacturer protocol. PCR conditions for each microsatellite marker are detailed in Supplementary Table S1. PCR products were verified in an agarose 1% gel stained with 1% of RedGel (Biotium, Fremont, CA, United States) and amplified in an ABI 3730 (Applied Biosystems, Madrid, Spain) in Unidad de Genómica (Universidad Complutense, Madrid, Spain). Scoring was performed using GeneMarker v. 2.61 (SoftGenetics, State College, PA, United States).

Five viable seeds per plant were randomly selected and weighed individually in a Mettler Toledo MX5 microbalance (1 μg precision; Mettler Toledo, Columbus, OH, United States), to assess family mean seed mass. In October 2012, simulating germination following natural dispersal, seeds of each maternal family were randomly selected and sown in common conditions in an experimental glasshouse at the CULTIVE facilities at URJC (690 m a.s.l). Seeds were sown in 1.4 L pots (10 cm × 10 cm × 17 cm) filled with locally collected topsoil and topped with a 4 cm layer of sieved gypsum substrate. This gypsum layer was collected from a single patch at the study site by removing the top layer of substrate to avoid the naturally occurring seed bank, and was added to improve germination. Ten pots per family were included in the experiment, and seven seeds were sown in each pot. During the germination period, pots were periodically watered to ensure permanent humidity, and temperature in the experiment was controlled to match the optimum reported for this species. The 10 pots with plants from the same family were randomly assigned to each of five blocks (two pots per maternal family; 240 pots per block) and arranged in a randomized complete block design. After 5 weeks, when no significant new emergence was observed, emerged seedlings were clipped to leave one plant per pot. Due to low germination of a few maternal families of some populations and several missing replicates, the final sample was 1099 plants (9–12 families/fragment × 10 fragments × 10 replicates). See Supplementary Table S3 for details on number of families and plants per fragment.

All plants were grown for ≈1 year under favorable conditions of water and nutrient availability, and several functional traits as well as survival were measured during this period to assess genetic variation and population differentiation (see below). After this period, we performed the plasticity experiment. We assessed plasticity to drought because water is the most limiting environmental factor in the Mediterranean; indeed, drought has been shown to be the main cause of seedling mortality in natural conditions in other coexisting gypsophiles (e.g., Escudero et al., 2000). In August 2013, simulating a typical Mediterranean summer drought, half of the remaining plants were assigned to a drought treatment and the other half was maintained in well-watered conditions (total N = 727 plants). In the drought treatment, soil moisture was slowly reduced by spacing out watering events (see Supplementary Figure S1). Soil moisture (expressed as % of field capacity for this substrate) was monitored weekly by weighing a random sample of 15 pots in each watering treatment. Plasticity to drought (details on measurements below) was evaluated when plants experienced a reduction of ≈50% of field capacity. After this, irrigation of plants in the drought treatment was ceased until all plants died.

Genetic diversity was high in all populations (Table 2). The eight microsatellite loci scored gave a total of 130 alleles in the 199 individuals of C. hyssopifolia successfully amplified, with an average of 16.25 alleles per locus. All loci were polymorphic with 4–33 alleles per locus. The average number of alleles observed per locus ranged from 6.30 in population 136 to 8.50 in population 133 (Table 2). Expected heterozygosity values were generally high and did not vary much among populations (from 0.61 to 0.67; Table 2). Observed heterozygosity ranged from 0.52 to 0.69 and was slightly lower than expected heterozygosity (Table 2). The inbreeding coefficient (F_IS) was very low, ranging from 0.02 to 0.09, and not significantly different from zero (Table 2). Population pairwise F_ST values were low yet significantly different from zero, with values ranging from 0.0113 to 0.104 (Supplementary Table S4).

All the measured plant traits (except for height, SLA and chlorophyll fluorescence) showed demonstrable genetic variance (significant Family term in Table 3), indicating that these traits have the potential to respond to selection. Narrow-sense heritability was variable, ranging from 0.119 (chlorophyll fluorescence, F_v/F_m) to 0.934 (emergence rate), and was significantly different from zero in most cases (entire confidence intervals above zero; Table 4). Similarly, the coefficient of genetic variance (CVA) and the percentage of phenotypic variance explained by family differences were lower for chlorophyll fluorescence and larger for emergence rate (Table 4). When evaluated for each population, we found no consistent pattern for the amount of genetic variance within populations. All populations had significant genetic variation for several traits (significant Family term, Supplementary Figure S2), but the percentage of variance explained by the factor Family varied widely across traits and populations.

There was also significant genetic variation among populations, i.e., population differentiation, for most traits (significant Population term, Table 3 and see Supplementary Table S5 for population means). However, differences among families within populations were higher than differences among populations (except for RGR), i.e., differences among families explained a larger proportion of phenotypic variance (range 3–23.4%) than differences among populations (range 1.3–18.9%), and this was very consistent among traits. For all traits, however, the great majority of the variance was explained by differences among individuals (range 68.2–95.5%; Figure 1).

The drought treatment exerted a significant effect on all traits except for stomatal conductance (significant Water Treatment term, Table 5). Plants from all populations significantly decreased leaf size (78% reduction) and number (28%), SLA (24%), plant height (16%) and chlorophyll fluorescence (5%) in the drought treatment (Figure 2). This plastic response was similar across populations (non-significant Population × Treatment interaction, Table 5), i.e., there was no population differentiation for plasticity. Drought also had a significant effect on plant survival: plants in dry conditions died faster than those that were under well-watered conditions (log-rank test statistic = -5.20, P < 0.0001, df = 1; Supplementary Figure S3). Again, this response was similar for all populations: there were no significant differences among populations in survival in the drought treatment (χ² = 6.87, P = 0.65; df = 9).

There was no significant effect of fragment size on the genetic diversity indices per population (Supplementary Table S6). Similarly, there was no significant effect of fragment size and connectivity, or their interaction, on functional traits, i.e., the observed phenotypic differentiation among populations was not related to the fragment characteristics of each population (Table 3). Furthermore, we found no effect of fragment size, fragment connectivity or their interaction on the amount of genetic variation (% of phenotypic variance explained by family differences) in each population, or in the expression of phenotypic plasticity, i.e., populations from small/isolated fragments did not show lower genetic variation or less plasticity than populations from large/connected ones (Supplementary Tables S7, S8). However, there was a significant difference in survival among populations that was related to fragment connectivity. Plants from populations from more connected fragments had higher survival than those from populations from isolated fragments (χ² = 32.15, P < 0.001, df = 9; Figure 3).

We found no significant correlation between neutral and quantitative genetic variation (% of phenotypic variance averaged for all traits) and population fitness (survival; r = 0.173, P = 0.633; Supplementary Table S9). Similarly, we found no relationship between phenotypic plasticity (averaged for all traits) and survival (r = 0.004, P = 0.992) or between plasticity and genetic variation (r = -0.324, P = 0.368). See Supplementary Table S7 for tests at the trait level.

Except for two traits, no significant seed mass effects were detected in functional traits, population differentiation or the expression of plasticity (not significant seed mass term, Tables 3, 5).

Populations of the study species from a fragmented landscape showed high neutral genetic variation (comparable to other congeneric species; see, e.g., Fréville et al., 2001; Marrs et al., 2008), and revealed significant genetic variation in life-history, morphological and physiological traits. Furthermore, plants from all populations responded plastically to water stress –the most limiting environmental factor in this system– in patterns consistent with adaptive functional plasticity. Despite large similarities among populations in the amount of neutral and quantitative genetic variation and on the expression of phenotypic plasticity, populations significantly differed in a major fitness component, survival.

Several studies have shown a positive relationship between genetic variation – either neutral and/or quantitative – and plant fitness components (e.g., Hansson and Westerberg, 2002; Leimu et al., 2010). For instance, in populations with low genetic variation, increased inbreeding and allele loss may lead to an increased proportion of homozygotes and the expression of deleterious alleles, which may eventually reduce plant survival (Young et al., 1996; Leimu et al., 2010). Therefore, populations with high genetic variation are expected to also have high fitness. Similarly, it is also an accepted paradigm that (adaptive) phenotypic plasticity has a positive impact on plant fitness (Sultan, 1995; Pigliucci and Schlichting, 1996; Matesanz and Valladares, 2014), and a few studies have experimentally confirmed this expectation (e.g., Lazaro-Nogal et al., 2015). However, in our study, populations showed similar levels of neutral and quantitative genetic variation, and broadly similar patterns of functional plasticity, but significantly differed in their survival. In other words, the population’s genetic variation and evolutionary potential of important functional traits (and their plasticity) were not related to population fitness. Altogether, and contrary to expectations, these results suggest that genetic variation –either molecular or quantitative– and phenotypic plasticity may be in some instances poor predictors of population survival, and are in agreement with other studies where correlations between quantitative genetic variation and various components of plant fitness were not found (Lauterbach et al., 2011; Walisch et al., 2015). More generally, our findings call for caution on the use of standing genetic variation and/or plasticity, as is usually done in biological conservation contexts, to forecast the future performance of populations.

The high genetic variation observed for most functional traits indicate that this gypsum specialist has substantial evolutionary potential, a first requirement to adapt to shifting environmental conditions. Our estimates of narrow-sense heritabilities largely varied among traits (0.119–0.934) but were moderate (<0.5) for 8 of the 10 traits measured (average h²: 0.44), and are in line with those reported for other species from fragmented landscapes (Mannouris and Byers, 2013; Weber and Kolb, 2014; Ye et al., 2014; Walisch et al., 2015). Lack of significant genetic variation for key traits such as plant height and chlorophyll fluorescence may be due to several factors, including low phenotypic variation due to the species rosette growth form and the lability of ecophysiological traits (see e.g., Agrawal et al., 2008). Only for two traits (emergence rate and leaf area) did we obtain very high values. It has been argued that in open-pollinated plants from mainly self-incompatible plants, a small proportion of the offspring can share the same father, leading to families being composed of both full- and half-siblings, which could bias the estimation of heritability (Falconer and Mackay, 1996). However, the generally high observed heterozygosity observed in these populations together with the fact that only two of 10 traits showed such a pattern suggests that high relatedness among siblings may not be playing a significant role in our calculations, and that outcrossing is indeed the main reproductive mechanism in our populations. Maternal effects can also lead to higher phenotypic similarity among siblings, particularly during early-life stages (Roach and Wulff, 1987). In our study, most of the functional traits were measured more than 6 months after germination, and seed mass effects were negligible in most cases. This suggests that maternal effects are likely to be minor in our study, and supports that the phenotypic differences observed have a genetic basis. We note, however, that unmeasured maternal effects may still have been present on traits expressed very early in life.

Alongside the high within-population genetic variation, plants from all populations responded adaptively to water stress. In addition to the passive response that reflects resource limitation, which included reduced plant size and photosynthetic performance, plants from all populations responded actively to water shortage (a reduction of 50% of soil field capacity) by reducing leaf size as well as the structure of the leaf tissue (SLA). The latter morphological and physiological adjustments are well-understood responses to moisture-limited conditions (Sultan and Bazzaz, 1993; Matesanz and Valladares, 2014). These results contribute to the limited current understanding on the ability of restricted edaphic endemics to express adaptive plasticity to varying environmental conditions (Damschen et al., 2012).

We did not find any effect of habitat fragmentation, i.e., fragment size and connectivity, on neutral and/or quantitative genetic variation of the study populations. All populations, irrespective of the size and degree of isolation of their source fragment, showed similar levels of neutral genetic variation and had significant quantitative variation – estimated by genetic differences among families– for several functional traits. In other words, populations from small and isolated fragments contained similar amounts of neutral and quantitative genetic variation, i.e., evolutionary potential, than populations from large and connected ones. Similarly, there was no effect of fragment size or connectivity on the measured functional traits or their plasticity.

Despite the lack of effects of habitat fragmentation on trait means, genetic variation and plasticity, we found a significant effect of fragment connectivity on a major fitness component, i.e., survival. Plants from populations from more connected fragments had higher survival in well-watered conditions, which can have important consequences for the long-term viability of populations in isolated sites. Low plant fitness in plants from isolated fragments can be the result of high inbreeding (Lienert, 2004; Kolb, 2005), the fixation of deleterious alleles (Young et al., 1996), seed provisioning (Stanton, 1984) and/or maternal effects (Roach and Wulff, 1987). Inbreeding and deleterious genetic effects seem unlikely in our case due to the low inbreeding coefficients and the high genetic variation observed. Seed provisioning and other maternal effects may have an impact on offspring survival. Matesanz et al. (2015) found lower seed mass on plants from small and isolated fragments. However, the almost negligible effect of seed mass on functional traits, genetic variation and plasticity together with the fact that significant mortality only started when plants were older than one year makes it unlikely that these factors are responsible for differences in survival among populations. Irrespective of the underlying mechanism, our results contribute to the mounting evidence that both fragment size and connectivity have an effect on fitness of C. hyssopifolia, either on its survival or reproductive output (Matesanz et al., 2009; Pias et al., 2010; Matesanz et al., 2015).

Lack of fragmentation effects on genetic variation, functional traits and plasticity is unexpected for several reasons. First, a comprehensive meta-analysis showed that self-incompatible species are more susceptible to the detrimental effects of habitat fragmentation than self-compatible ones because of their dependence on pollinator mutualisms (Aguilar et al., 2006). Second, recent studies reported high vulnerability (e.g., lower survival and reproductive output) of C. hyssopifolia populations to different components of habitat fragmentation (Matesanz et al., 2009, 2015; Pias et al., 2010). In particular, a field assessment of the fitness effects of fragment size and connectivity performed in the same study site showed lower reproductive output (i.e., fewer viable seeds per capitulum and lower seed set) of plants from isolated fragments (Matesanz et al., 2015), which could lead to lower genetic variation if fewer offspring are contributed to the next generation, and thus effective population sizes decrease over time. Finally, it could be expected that the anatomical features of C. hyssopifolia seeds, with a short pappus compared to seed weight, would limit their dispersal ability (Sheldon and Burrows, 1973), in turn decreasing gene flow and increasing the effects of genetic drift (Ellstrand and Elam, 1993; Young et al., 1996; Leimu et al., 2006).

Contrary to these expectations, our results therefore suggest that gene flow –either via pollen or seeds– among Centaurea populations is sufficient to maintain similar levels of genetic variation, even in the more isolated fragments. This is consistent with the low population differentiation observed in functional traits, which was in most cases lower than the relatively high within-population variation (Figure 1). Unrestricted gene flow is also supported by generally low inbreeding coefficients and low population differentiation estimated with microsatellite markers, indicating that in these populations, conspecific pollen from various sources may be readily available and breeding among relatives may be rare. Altogether, these results suggest that, at the fine scale, the sampled fragments of this species might functionally constitute a metapopulation with similarly connected nodes of different sizes and isolation levels (Manel et al., 2003). Alternatively, although the generation time of the study species is much shorter than the time since the fragments originated, we cannot rule out the possibility that the effects of fragmentation on neutral or quantitative genetic variation will be observed in the long-term (Debinski and Holt, 2000; Ellmer et al., 2011; Mona et al., 2014). Indeed, a lag effect of habitat fragmentation has been observed in several plant and animal species (e.g., Reitalu et al., 2009). Furthermore, various mechanisms other than gene flow, such as overdominance, balancing selection or mutation rates can be responsible for the similar genetic variation observed across populations. A few studies have reported a significant relationship between population size and quantitative variation (Willi et al., 2006). For instance, Weber and Kolb (2014) found higher heritabilities of phenological and reproductive traits in larger populations of the perennial herb Phyteuma spicatum. In contrast, our results concur with other studies where no relationship between genetic variation and population and/or fragment size was detected (see e.g., Ellmer et al., 2011; Walisch et al., 2015).

Similarly, there was no effect of fragment size or connectivity on the measured functional traits or their plasticity. Size, growth rate and leaf traits were similar in plants from all populations. Phenotypic differentiation among populations can be the result of differential selective pressures in each population, the effect of neutral processes such as genetic drift and inbreeding, or a combination of both (Merilä and Crnokrak, 2001; Holderegger et al., 2006). For instance, Walisch et al. (2015), in a study of the rare Saxifraga rosacea subsp. sponhemica from fragmented populations, found that most population trait means were significantly related to climate gradients, suggesting adaptive genetic differentiation in quantitative traits (see also Weber and Kolb, 2014; Ye et al., 2014). Likewise, Kittelson et al. (2015) found that inbreeding significantly affected the expression of 12 functional traits in populations of the prairie forb Echinacea angustifolia, and Mannouris and Byers (2013) found a significant effect of genetic drift load on populations of Chamaecrista fasciculata from small fragments. In our study, lack of population differentiation in functional traits suggests a low impact of genetic drift and inbreeding on plant traits, and/or similar selective pressures acting on our populations due to similar microclimatic conditions (see also Ellmer et al., 2011). We cannot discard, however, the possibility that populations have adaptively diverged to their local environment in a different set of unmeasured functional traits, or that inbreeding effects will be expressed in more stressful field conditions (Keller and Waller, 2002).

Lack of population differentiation in plasticity patterns may be due to past selection on specific norms of reaction, or to absence of within-population genetic variation for plasticity to evolve (Matesanz et al., 2010). Irrespective of the underlying mechanism, parallel populations’ norms of reaction suggest that fragment size and connectivity had no effect on the expression of phenotypic plasticity to varying water conditions. Our results concur with a study of Carlina vulgaris from populations of different size and degree of isolation that found no evidence that small or isolated populations were less plastic than large or connected ones (Berg et al., 2005 but see Kery et al., 2000; Paschke et al., 2005). More generally, these results contribute to our understanding on the yet poorly explored effects of habitat fragmentation on plasticity.