Cistus laurifolius L. subsp. laurifolius is a shrub with large, white, hermaphrodite flowers, pollinated by generalist insects (e.g. beetles, bees and flies), flowering from May to June (Astuti et al. 2015). This taxon has relatively superficial roots and nutrient-poor leaves, mainly adapted to oligotrophic and degraded soils, as other Cistaceae (Moro et al. 1996). Cistus laurifolius typically occurs in open and dry habitats of non-coastal Mediterranean areas (Rizzotto 1979; Grossoni and Venturi 2009; Roma-Marzio et al. 2016b). It shows relatively small seed mass (about 1 mg/seed), a feature related to easier penetration and accumulation of seeds in the soil (Fenner 1985; Thanos et al. 1992). Seed dispersal is seemingly barochorous (Thanos et al. 1992) but, despite this feature, long-distance dispersal events have been documented in the genus (Rizzotto 1979).

Concerning the native geographic distribution, C. laurifolius occurs in the western (Morocco, Portugal, Spain, France) and eastern Mediterranean basin (north-eastern Greece and Turkey). An isolated population is found in Central Italy (Tuscany) (Warburg 1968; Rizzotto 1979; Fernández-Mazuecos and Vargas 2010, Dimopoulos et al. 2013).

The only Italian population of Cistus laurifolius L. subsp. laurifolius is located in Tuscany, near the village of Santa Brigida (Firenze). The mean annual temperature and mean annual rainfall of the area, measured between 1992 and 2010 by a thermopluviometric station located in Pontassieve (WGS84: 43.812324, 11.399167; 120 m a.s.l.) are 13.7°C and 856.2 mm, respectively (http://agrometeo.arsia.toscana.it/).

With the exception of a single individual surviving in the near proximity of the village of Santa Brigida, the population is fragmented into four sub-populations (Astuti et al. 2015; Roma-Marzio et al. 2016b) (Figure 1). Three of them (C5, D and FOR in Figure 1) are located in an open habitat (garrigue) dominated by Cistus salviifolius L., Erica arborea L., E. scoparia L. subsp. scoparia and Cytisus scoparius (L.) Link subsp. scoparius. The remaining one (MAS in Figure 1) is located in the underbrush of a Pinus nigra L. subsp. nigra plantation, together with sclerophyllous vegetation, mainly composed of Cistus salviifolius L., Quercus ilex L. and Rubus ulmifolius Schott (Roma-Marzio et al. 2016b).

This fragmentation has most probably been caused by human induced landscape transformation, as is suggested by the ongoing disappearance of the sub-population in the near proximity of the village of Santa Brigida and by the occurrence of buildings (e.g. farms) and a network of roads surrounding the population (Figure 1).

Figure 1.

Distribution map of the Italian population of C. laurifolius L. subsp. laurifolius. Yellow symbols indicate the location of each sub-population. Circle size is proportional to the area occupied by the sub-populations of Fornellaccio (FOR), Fontassenzio (D) and west of Fornellaccio (C5). The yellow star, not proportional to the area of occupancy, refers to the small sub-population of Masseto (MAS). The exclamation mark (bottom right) indicates the single individual surviving in the near proximity of the village of Santa Brigida (not investigated in this study). In the top left corner, the location of the study area (red square) in Italy is indicated.

To estimate the population size (number of individuals), its density and demographic structure, an evaluation of the area occupied by each sub-population was carried out. To this end, preliminary data published by Grossoni and Venturi (2009) were verified and adjusted by field surveys, using a GPS receiver. Thereafter, each sub-population’s area of occupancy was determined within a GIS environment (QGIS software v. 2.18; QGIS 2016), by delimiting a minimum convex polygon. All individuals were counted in the smaller sub-populations (C5, D, MAS), while, for the largest sub-population (FOR), all the individuals occurring in fifteen randomly sampled 10 × 10 m plots were counted. The number of individuals in FOR (I_FOR) was estimated according to the following formula:

where I_p is the total of individuals occurring in all the plots, A_FOR is the total area occupied by the sub-population FOR and A_p is the sum of the areas of each plot (1500 m²).

The age in shrub plants is hardly detectable and plant growth is highly dependent upon environmental and ecological parameters (Chapman 1986). Therefore, to describe the demographic structure of each sub-population, three practical “age-stage” classes have been arbitrarily selected, based on the presence/absence of flowers and the number of leaves. According to the terminology proposed by Gatsuk et al. (1980), the three classes were defined as following: class I) immature plants (lacking flowers, showing less than six leaves and with a plantlet-like appearance; class II) virginile plants (lacking flowers, but showing more than six leaves; class III) mature plants (showing flowers).

It was not possible to take into account the very small and ephemeral cotyledon phase (seedling), due to difficulties in detecting this stage in the soil and leaf litter and also because its presence could be easily affected by exceptional weather conditions (Alonso et al. 1992).

For each sub-population, reproductive fitness was evaluated by means of seed set (number of seeds/number of ovules) and seed mass. The number of ovules was averaged on ten randomly selected flowers in each sub-population (for a total of 40 flowers). The ovaries were dissected along their septa using a razor blade in order to count the ovules under a 60× magnification stereomicroscope.

Since no difference was found for mean ovule numbers amongst sub-populations (ANOVA, p > 0.05), then the mean ovule number at population level (94.72) was used as the reference to calculate the seed set for all sub-populations. In the case of single fruits showing a seed number exceeding the mean ovule number, a seed set of 100% was assigned by default. To calculate the seed number, 50 fruits (capsules) were randomly collected for each sub-population and the seed number was counted for each capsule. Aborted (i.e. showing a seed-shape, but lacking embryo) seeds were not taken into account. To evaluate the seed set, the data were averaged at the sub-population level as follows:

where n (=50) is the number of sampled fruits for each sub-population and x is the seed set calculated for each fruit.

To evaluate the seed mass, ten replicates, each consisting of a group of 50 randomly selected seeds, were weighed (± 0.001 mg accuracy), for a total of 500 seeds per sub-population.

To characterise the floral sex allocation and to infer the breeding system of the species at the population level, flowers were collected at an early developing stage (showing mature, but still not-dehiscing anthers), in order to evaluate the flower biomass (mean dry weight) and the mean number of pollen grains produced per flower (P).

According to the regression formula proposed by Herrera (1985), the dry flower weight was used as an indirect measure of the daily nectar secretion rate (Bosch 1992). Dry weight was calculated on 34 flowers (ten from the sub-population FOR, ten from C5, ten from MAS and four from D). After their collection, the flowers were dried at 60°C for 24h and then weighed (± 0.001 mg) (Kay and Picklum 2013).

The total pollen production per flower was estimated according to the dilution method proposed by Galloni et al. (2007). Pollen grains were estimated on seven flowers (two from the sub-population FOR, two from C5, two from MAS and one from D). Anthers were taken from fresh flowers and placed in vials with 1 ml of ethanol/glycerinated fuchsine-glycerol solution (3:1). For effective pollen release, vials were sonicated for 1 min at 14 kHz using a Sonoplus Ultrasonic Homogeniser GM 2070. A known volume (1 µL) of suspension was mounted on a microscope slide within 10 seconds from sonication, to guarantee the homogeneity of suspension.

Pollen grains were then counted using a light microscope (250× magnification) and recorded with the help of a manual cell counter. The number of pollen grains was finally multiplied by the dilution factor and then by the number of anthers to obtain the total number of pollen grains estimated for a whole flower. Finally, according to the values indicated by Cruden (1977), the P/O ratio was used to infer the breeding system of the species. For the number of ovules, the same mean value already used for the seed set was referred to, as explained in the previous section.

To evaluate the overall effect of three single predictors (logarithm of area of occupancy; sub-population density and % of adults) on the seed set, each of them was fitted with a single Generalised Linear Model (GLM), with a logit link function and a binomial error structure, followed by a likelihood test. The logistic regression was selected since the seed set is a binomial phenomenon.

The values of the area of occupancy were subjected to logarithmic transformation to reduce the large differences amongst sub-populations.

Differences in seed mass, flower mass and ovule number amongst sub-populations were tested by means of an ANOVA test, followed by Tukey’s pairwise comparisons, after checking normality and homoscedasticity of the data. Differences in seed set amongst sub-populations were tested by means of χ² test.

For all statistical tests, significance was accepted at p ≤ 0.01. All analyses were performed using R 3.3.1 software (R Core Team 2016).

Microsatellite markers have been developed according to Albadejo (2010). About 1500 clones from a non-enriched genomic library were sequenced. The sequencing reads were assembled with CodonCode Aligner in 1348 unique sequences. Sequences were checked for the occurrence of di-, tri- and tetra-nucleotide repeats with the online software Sputnik (available at http://wheat.pw.usda.gov/ITMI/EST-SSR/LaRota). Although microsatellite motifs were detected in 25 sequences (1.8%), ten were discarded as the microsatellite motifs were too short, nucleotide repeats were too close to the vector for primer design or the clones showed high sequence homology. Subsequently, 15 primer pairs were designed by using Primer3 software (Rozen and Skaletsky 2000). In a preliminary screening of the markers’ variability, 30 individuals from two sub-populations were examined. Out of the 15 primer pairs, nine were discarded, as they failed to amplify, produced multi-banding patterns, were monomorphic or showed too pronounced stuttering. The presence of null alleles in the remaining six markers (Suppl. material 1) was examined following the Expectation Maximisation (EM) algorithm (Dempster et al. 1977) using FREENA (Chapuis and Estoup 2007). The estimated frequency of null alleles ranged from 0.13 to 0.23 with the exception of cislau11 (frequency = -0.01). Due to the excess of null alleles (frequency = 0.78) the marker cislau5 was excluded from further analyses.

The total DNA from 189 plants, sampled from the four sub-populations, was isolated using the Qiagen DNeasy Plant Mini Kit (QIAGEN, Hilden, Germany) from 80–100 mg of leaf dry tissue. Amplifications were performed by polymerase chain reaction (PCR) in 10 µL volumes, containing 10-50 ng of template DNA, 1× reaction buffer (200 mM Tris-HCl, 500 mm KCl, pH 8.4; Invitrogen), 0.5 U of Taq polymerase (Invitrogen), 0.5 µL of 1% W-1 solution (Invitrogen), 2 mm of MgCl₂, 1 µm of each primer, 60 µm of dNTP mix.

Reactions were performed in a Gene Amp PCR system 9700 (PE Applied Biosystems), with the following programme: an initial denaturation step of 3 min at 94° C, followed by 10 touchdown cycles of 30 s at 94 °C, 40 s at 60 °C (1 °C lower per cycle) and 30 s at 72 °C and 25 cycles of 20 s at 94 °C, 20 s at 50 °C and 30 s at 72 °C with a final extension step of 8 min at 72 °C. A final extension of 6 min at 72 °C was performed in all programmes. Amplified fragments were run in an ABI 3130xl automatic sequencer (Applied Biosystems). Electropherograms were analysed using GeneMapper version 4.0 (Applied Biosystems). Linkage disequilibrium between loci and deviations from Hardy-Weinberg (HW) expectations were tested using Fisher’s exact tests based on Markov chain procedures in GENEPOP ver. 3.4 (Rousset 2008). Basic statistics were calculated using the software GENALEX 6.2 (Peakall and Smouse 2006), to determine: allele frequencies, mean observed heterozygosity (H_o), unbiased expected heterozygosity (H_e), the number of alleles at each SSR locus (N_A), the effective number of alleles (N_E) and the inbreeding coefficient (Fis). GENALEX 6.2 was also used to perform the analysis of molecular variance (AMOVA). Significance levels were determined using 1,000 permutations. To test whether any of the sampled sub-population experienced bottlenecks in the recent past, the Bottleneck programme version 1.2.02. (Piry et al. 1999) was used to compute the difference (averaged over loci) between actual heterozygosity and the heterozygosity expected for a population in mutation-drift equilibrium. The Bottleneck software allowed the testing of the bottleneck hypothesis, taking into account two possible models: a stepwise mutation model (SMM) and a two-phase mutation model (TPM). The significance of heterozygosity excess was determined using the Wilcoxon signed-rank test, which has been demonstrated to be the most accurate test in case of both low number of polymorphic loci (< 20) and small samples sizes (< 30) (Piry et al. 1999). The genetic structure of the four sub-populations was analysed by the Bayesian algorithm implemented in STRUCTURE v. 2.3.3 (Pritchard et al. 2000), which assigns individuals to a K number of genetically homogeneous groups, based on allele frequencies at each locus. For the analyses with STRUCTURE, a burn-in period of 50,000 and a posterior number of Markov Chain Monte Carlo (MCMC) of 100,000 permutations was used.

Fifteen replications (runs) were performed for each value of K ranging from K = 1 to K = 10. An admixture and allele frequencies correlated model was used. The most likely number of genetic clusters (K) was estimated following Evanno et al. (2005), which uses an ad hoc parameter (ΔK) to estimate the rate of change of likelihood values amongst successive K values. Ten runs for each simulation were averaged using algorithms found in CLUMPP (cluster matching and permutation programme; Jakobsson and Rosenberg 2007) and represented as bar graphs using DISTRUCT (Rosenberg 2004). The membership probability of each individual in every cluster was assessed by the value of Q and each individual was assigned to a specific cluster taking into account a threshold of Q > 0.75 (Atiqur et al. 2016).

The estimated number of individuals within the population was 9,962, occupying an area of 86,145 m². The largest sub-population FOR hosted the vast majority of the plants (Table 1). Overall, mature plants appeared in the greatest numbers, followed by virginiles and immatures. Percentages similar to the overall values were found in the sub-populations FOR and MAS, whereas C5 and D showed a higher number of immature and a lower number of mature individuals. In addition, D showed the highest percentage of virginile plants (Table 1).

According to χ² test, the seed set of FOR (mean value 87.5%) and D (mean value 38.9%) showed the highest and the lowest values, respectively (Table 2). Regarding the seed mass, differences were significant only between FOR and MAS, with the latter sub-population showing the lowest values (Table 2), yet not significantly different from C5 and D.

The mean P/O value calculated for the population was 5,138.72 ± 4,310. Concerning dry flower mass, no differences were found amongst sub-populations (p > 0.05): the mean value for the entire population was 102.05 ± 0.02 mg. According to Herrera’s regression formula, the daily nectar’s production was estimated as 5.18 ± 0.78 mg per day per flower.

Concerning factors affecting seed set, a significant positive effect was found of the sub-population’s area of occupancy, density and frequency of class III (mature plants) (Table 3).

A total of 66 alleles was detected for five loci (Suppl. material 2) in the 189 individual genotypes sampled in this analysis. The overall diversity in all sub-populations was very low (Table 4); the mean number of alleles per locus (N_A) was 3.3 ± 0.5, the observed heterozygosity (H_O) was 0.25 ± 0.05, whereas expected heterozygosity (H_E) was 0.32 ± 0.06. No sub-population showed a significant excess of heterozygotes compared with the equilibrium expectation based on data modelling (one-tailed Wilcoxon signed-rank test, p > 0.80 for all sub-populations under both the SMM and TPM). Instead, all sub-populations showed a significant excess of homozygotes (as indicated by a positive F_IS), with the exception of MAS, where the F_IS value was slightly negative.

Linkage disequilibrium (LD), the non-random association of the alleles at different loci, was analysed for all pairs of SSR markers within each sub-population and across the whole population. Only one locus (cislau12) showed a significant departure from equilibrium (5% level) in the C5 sub-population. The population structure determined by AMOVA showed that approximately 4% of the total variation was attributable to variation amongst sub-populations and 96% of the total variation was attributable to differences amongst individuals within sub-populations. All pairwise F_ST values differed significantly from zero (p < 0.05), except between C5 and MAS (p = 0.252) and D and MAS (p = 0.076). The estimate of overall F_ST was significantly different from zero, but very low (F_ST = 0.050; 95% CI 0.022-0.083), suggesting strong inter-subpopulation gene flow. Concerning this result, the mean value of N_M (number of migrants) was estimated to be 7.25. The optimum cluster number inferred from the STRUCTURE analysis (Suppl. material 3) was K = 3. Further analyses were performed based on K = 3, to investigate the composition of each individual and each sub-population with respect to the three inferred genetically homogeneous groups. In agreement with the low levels of genetic diversity, the genetic structure was also very weak. The proportion of membership of each inferred genetically homogeneous group was, in all cases, lower than 52%; the major component of genetic composition is attributable to the first K group (51%) for D, to the second K group for FOR (41%) and to the third K group for MAS (50%) (Figure 2). This result demonstrates that each gene pool shows a high degree of admixture, pointing towards extensive gene flow or common ancestry.

Figure 2.

Estimated genetic clustering (K = 3), obtained with STRUCTURE analysis of 189 individuals from the Italian population of Cistus laurifolius. Each individual is represented by a vertical line, which is partitioned into coloured segments, indicating the individual’s estimated membership fraction in K clusters. Different sub-populations are separated by a vertical black line.

The survival chance of the Italian population of the laurel-leaved rock rose is mostly dependent upon the fate of the largest sub-population (FOR), including the vast majority of the individuals. The demographic structure of the population, as well, clearly parallels that of the largest sub-population, where a high percentage of mature plants and a low percentage of virginile and immature plants was observed (Table 1). The sub-populations of intermediate size and low density, namely C5 and D, show higher percentages of virginile and immature plants compared to the largest and to the smallest sub-populations (FOR and MAS, respectively). Accordingly, C5 and D may provide an important source of recruitment, albeit the seed set of these two sub-populations was relatively low, possibly due to the Allee Effect. As a consequence, C5 and D sub-populations can be regarded as patches of a fragmented habitat, surrounding the largest sub-population (FOR). Taking into account the smallest and most isolated sub-population (MAS), its peculiar environmental conditions should be underlined. As reported by Roma-Marzio et al. (2016b), MAS is subjected to severe threats (e.g. canopy closure and alien species) which are significantly affecting its ecological asset. Only 40 mature plants are surviving in this sub-population and only one immature plant was observed during this survey. Despite the low number of individuals in MAS, the seed set value was not low, probably due to the high density of mature individuals, as highlighted by the GLM results. Although there are relatively high values of the seed set, the low germination rates of this species (Fernández-Mazuecos and Vargas 2010) could however represent an additional threat.

According to the categorisation of P/O values made by Cruden (1977), the studied population seems fully xenogamous. This is also confirmed by the estimated daily nectar production which is consistent with a xenogamous breeding system, according to the data presented by Bosch (1992). The mating system is reported to have a significant impact on the distribution of genetic variability in a plant population (Duminil et al. 2007) and this genetic study indicates that 96% of the detected genetic variation is maintained within each sub-population. This value complies with that which is expected for an outcrossing plant (Hamrick and Godt 1989). Given that outcross seems the only way to produce new individuals in this species, its capability to cope with environmental shifts is strictly dependent upon the reproductive fitness of mature plants. The number of mature plants and the seed set measured in the population are quite high (Tables 1 and 2). In the light of these data, it seems that the Italian population of C. laurifolius shows a good reproductive performance. However, further investigations on the mating system using molecular markers could provide useful information (Glémin et al. 2006).

Compared with values reported for other outcrossing species (H_O = 0.63; H_E = 0.65) and for long-living perennials (H_O = 0.63; H_E = 0.68) (Nybom 2004), these results suggest that genetic diversity is low within the population as a whole (mean H_O = 0.25 ± 0.06; mean H_E = 0.31 ± 0.03). Low genetic diversity is commonly observed in endangered and rare species, with a reduced number of individuals in their populations (Hughes et al. 2003; Zhang et al. 2005), albeit different molecular markers can yield different diversity estimates (Nybom 2004). Although there are no concrete experimental data supporting a bottleneck hypothesis, the low genetic diversity, coupled with a relatively high F_IS (averaged value over all sub-populations: 0.21 ± 0.08), suggest that the sub-populations could have suffered a reduction in size (Frankham 2005; Wright et al. 2009). Furthermore, the loss of habitat and its fragmentation (Ellstrand and Elam 1993), as well as multi-generation isolation (Keller and Waller 2002), are also known to have significant effects on gene flow and genetic diversity, as a result of both drift and increased inbreeding levels (Aguilar et al. 2008). In particular, outcrossing species containing most of their genetic variability within sub-populations, are subjected to genetic erosion through habitat fragmentation (Frankham 2005). Furthermore, the proportion of short-distance mating events has been shown to increase as a consequence of an increased fragmentation process (Ismail et al. 2012).

There was evidence for inbreeding in all the sub-populations with the exception of MAS: the F_IS value of this small sub-population did not differ from Hardy-Weinberg expectations (Table 4). Although the F_IS value of MAS could be partially influenced by the low number of sampled individuals (N = 10, however representing 25% of mature plants), this further confirms the absence of a relationship amongst the sub-population size and the inbreeding coefficient, as already shown by other authors (e.g. Leimu et al. 2006; Honnay and Jacquemyn 2007; Aguilar et al. 2008). Age-structure of a sub-population is known to affect the Hardy-Weinberg expectations (Jacquemyn et al. 2004), therefore estimates of the inbreeding coefficient such as F_IS should also be considered with caution. Indeed, this sampling strategy involved only adult plants. Although no information is available in literature about the life span of C. laurifolius, it is speculated that the oldest individuals possibly established before habitat fragmentation (Evanno et al. 2005). This could explain the observed absence of effects of the sub-population size on the inbreeding coefficient. In this sense, minor sub-populations hosting more virginile and immature plants may play a significant role with respect to the survival of the last Italian population of C. laurifolius, despite no genetic differences amongst (mature individuals of) sub-populations were found. The lack of genetic structure could be explained by two main factors acting together: (i) a recent fragmentation of the population (F_ST has been shown to respond slowly to the fragmentation process) (Landguth et al. 2010) and (ii) xenogamous species behaviour. These two factors might have maintained gene flow amongst sub-populations, and gene flow is a very important factor counteracting loss of genetic diversity and inbreeding depression. Moreover, the long-living cycle of the species could be further delaying the genetic differentiation amongst sub-populations (Hamrick and Godt 1996; Colling and Matthies 2006).

Gene flow may be insufficient to counteract the effects of drift, especially at low levels of population density (Kettle et al. 2003; Rosas et al. 2011). For example, Ismail et al. (2012) have shown that pollen dispersal occurs over large distances in highly fragmented agro-forest landscapes. However, as forests become more fragmented, inbreeding due to non-random mating amongst related individuals, as well as genetic drift, are likely to be exacerbated. The low inbreeding in MAS may reflect the mating system preceding the fragmentation when the population was more continuous, so that extensive gene flow could take place. On the other hand, several additional factors could affect the relationship between F_IS and population size (Honnay and Jacquemyn 2007). The absence of homozygotes for rare alleles could bias downwards the F_IS value in small sub-populations (Young et al. 1999), whereas F_IS in large sub-populations could be biased upwards as a result of the Wahlund effect.

Whether recent or not, habitat fragmentation is one of the most important factors invoked to justify the low seed set values found in many natural plant populations, mostly due to its influence on pollination and genetic erosion (Ågren 1996; Young et al. 1996; Severns 2003; Ghazoul 2005; Leimu et al. 2006; Ouborg et al. 2006; Aizen et al. 2007; Aguilar et al. 2008). Indeed, as demonstrated in other species (Honnay et al. 2005), a negative effect of small areas of occupancy on the seed set was also detected, as well as a positive effect of plant density. High seed set values in all sub-populations point towards a good reproductive performance, but no data are available for the effective establishment of seedlings. Probably, the low number and density of mature individuals in the sub-population D accounts for the lower seed set observed (Allee effect). On the other hand, the high level of seed set shown by the large sub-population FOR parallels a low percentage of immature and virginile plants. In long-living perennials, the population fitness can be related to growth and survival rates which may reduce seed production (Silvertown et al. 1993, 1996). The low presence of immature and virginile plants in the largest sub-population FOR may be potentially deleterious for the population as a whole, since the lack of recruitment could lead to genetic depletion. In this perspective, immature and virginile plants supplied by the sub-populations C5 and D may be crucial for the long-term survival of the Italian population of laurel-leaved rock rose. On the other hand, a population structure characterised mostly by aged plants, in the worst scenario, may lead to local extinctions (García et al. 1999; Brys et al. 2003; Jacquemyn et al. 2003). The smallest sub-population (MAS) may experience such a trend within a relatively short period.