Research Article
Research Article
Conservation biology of the last Italian population of Cistus laurifolius (Cistaceae): demographic structure, reproductive success and population genetics
expand article infoGiovanni Astuti, Francesco Roma-Marzio, Marco D'Antraccoli, Gianni Bedini, Angelino Carta, Federico Sebastiani§, Piero Bruschi|, Lorenzo Peruzzi
‡ University of Pisa, Pisa, Italy
§ IPSP, Consiglio Nazionale delle Ricerche, Firenze, Italy
| Università degli studi di Firenze, Firenze, Italy
Open Access


Isolated populations are usually subject to low fitness and reduced genetic diversity, both of which may negatively affect their survival and adaptive potential. Hence, these issues cannot be neglected when planning conservation actions for isolated populations. The Italian population of Cistus laurifolius subsp. laurifolius is extremely isolated. Furthermore, it is affected by fragmentation, being constituted by a single larger subpopulation, surrounded by three much smaller subpopulations, a few hundred metres to a few kilometres apart. In order to fill gaps in demographic and genetic knowledge concerning the Italian population, its area of occupancy, size, age-stage structure and phenology were investigated and its reproductive fitness, pollination strategies and genetic variability were assessed. The population was inferred as fully xenogamous and showed good reproductive performance. Despite this, its genetic variability was low and it showed relatively high levels of inbreeding depression (Fis), seemingly not affected by sub-population size. These results suggest that the Italian population recently suffered fragmentation and reduction in size. The low genetic diversity observed could be explained by the high percentage of mature individuals found in the population, possibly established before fragmentation. For these reasons, the Italian population of C. laurifolius subsp. laurifolius should be monitored and concrete actions aimed at its conservation planned.


Conservation biology, genetic diversity, reproductive biology, Italy


The knowledge of relationships amongst population size, age-stage structure, fitness and genetic diversity is of crucial importance in plant ecology and conservation, in order to gain an insight into population dynamics and evolutionary potentialities (Stockwell et al. 2003; Rodriguez-Perez 2005; Leimu et al. 2006).

Population size strongly affects local adaptation (Leimu et al. 2008), whereas population structure can influence breeding systems and, eventually, reproductive outcomes (Carta et al. 2016a). Population dynamics are usually age-stage-dependent and related to the survival potential of a species (Harper 1977).

Small outcrossing plant populations may suffer disadvantages due to the Allee effect, i.e. an individual fitness reduction caused by a decrease in population size/density (Forsyth 2003). Populations with a low density may be affected by an insufficient pollinator service and consequently express a reduced fruit and seed set (Kunin 1997; Forsyth 2003; Ashman et al. 2008).

According to Ghazoul (2005), the vulnerability of plant species to the Allee effect depends upon their mating system. Self-incompatible species are more likely to experience pollen limitation than self-compatible plants, as the former cannot compensate for reduced pollinator services through selfing (Aizen and Feinsinger 1994; Knight et al. 2005; Xia et al. 2013). In addition, when peripheral and isolated plant populations (PIPPs) (Abeli et al. 2009) display reduced genetic variability, they may fail to effectively cope with environmental changes, especially in the light of dramatic shifts imposed by global climate change (Willi and Hoffmann 2009, Hoffmann and Sgrò 2011). The Italian peninsula hosts several examples of such populations, highlighting the importance of case studies for planning proper conservation actions (Gargano et al. 2007, 2009; Rossi et al. 2009; Carta et al. 2016b; D’Antraccoli et al. 2016).

The laurel-leaved rock rose, Cistus laurifolius L. subsp. laurifolius, can be regarded as a notable case for studying relationships between genetic diversity and conservation in the context of PIPPs. This species shows a distribution scattered across the Mediterranean, the main populations being located in the Iberian peninsula and south France in the west and Anatolia in the east (Fernández-Mazuecos and Vargas 2010). In between, C. laurifolius survives in a single population in Italy, near Santa Brigida village (Tuscany) (Roma-Marzio et al. 2016a, 2016b and literature cited therein). A single eastwards migration event may have been the source of a recent dispersal of the species from putative refugial areas in the western Mediterranean (Fernández-Mazuecos and Vargas 2010). According to the latter authors, narrow ecological requirements (altitude and edaphic conditions) and low germination rates may explain the current distribution pattern of this species.

No information about the genetic diversity, demographic structure and reproductive traits is available for the Italian population of laurel-leaved rock rose. This population has been considered as a relict, resulting from fragmentation (Dansereau 1939; Rizzotto 1979; Fernández-Mazuecos and Vargas 2010), due to the stenoecy of this species (Rizzotto 1979). The dramatic reduction of the Italian population during the last centuries was also determined by human impact (Astuti et al. 2015; Roma-Marzio et al. 2016b). According to the latter authors, two out of five sub-populations are much depleted and, taking into consideration the threats recorded in the area, this taxon has been evaluated as Vulnerable at Regional level, according to the IUCN criteria (IUCN 2014).

The aim of this study is to accumulate information useful for the conservation of this species, including: i) area of occupancy of the Italian population, ii) number of immature, virginile and mature individuals, iii) phenology, iv) reproductive fitness, v) pollination strategies and vi) genetic variability. The obtained results will be a framework to design a conservation programme for the species.


Study species

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).

Study area

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 (

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.

Population size and structure

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 (IFOR) was estimated according to the following formula:

where Ip is the total of individuals occurring in all the plots, AFOR is the total area occupied by the sub-population FOR and Ap is the sum of the areas of each plot (1500 m2).

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).

Reproductive fitness

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.

Floral sex allocation

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.

Statistical analysis of demographic structure and reproductive traits

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 χ2 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

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 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.

Genetic diversity analysis

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 MgCl2, 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 (Ho), unbiased expected heterozygosity (He), the number of alleles at each SSR locus (NA), the effective number of alleles (NE) 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).


Demographic structure

The estimated number of individuals within the population was 9,962, occupying an area of 86,145 m2. 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).

Table 1.

Demographic structure of the Italian population of C. laurifolius subsp. laurifolius. In square brackets: % values of area of occupancy and individuals with respect to the whole population. In round brackets: % age-stage classes with respect to the whole number of individuals for each sub-population. * = mean value.

Subpopulation Area (m2) Density (individuals/m2) Mature individuals Virginile individuals Immature individuals Total individuals
FOR 77,215 [89.7%] 0.12 8,185 (86.4%) 1,132 (12%) 154 (1.6%) 9,471 [95.1%]
C5 6,680 [7.7%] 0.05 261 (74.4%) 45 (12.8%) 45 (12.8%) 351 [3.5%]
D 2,100 [2.4%] 0.04 69 (75.0%) 17 (18.5%) 6 (6.5%) 92 [0.9%]
MAS 150 [0.2%] 0.31 40 (85.1%) 6 (12.8%) 1 (2.1%) 47 [0.5%]
Population 86,145 0.11* 8,556 (85.9%) 1,200 (12.0%) 206 (2.1%) 9,962

Reproductive fitness and floral sex allocation

According to χ2 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.

Table 2.

Reproductive features of the Italian population of C. laurifolius subsp. laurifolius. For each sub-population, mean values and standard deviations are reported. Different letters indicate growing ranking significant differences amongst groups (ANOVA test for ovule number, seed number and seed mass; χ2 test for seed set).

Sub-population N. ovules
(N = 10)
N. seeds
(N = 50)
Seed set
(N = 50)
Mass of 50 seeds (mg)
(N = 10)
C5 95.9 ± 27.3a 55.28 ± 37.2b 58.36%b 42.43 ± 1.7a, b
D 100.4 ± 24.7a 36.84 ± 27.0a 38.89%a 43.28 ± 1.5a, b
FOR 102.6 ± 28.0a 82.90 ± 39.5c 87.52%c 44.24 ± 1.4b
MAS 80.0 ± 20.7a 56.96 ± 40.0b 60.14%b 41.15 ± 2.0a

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).

Table 3.

Overall effect of the three predictors on the seed set, estimated by three single GLM analyses. log(area) = natural logarithm of the area occupied by each sub-population, III% = percentage of the class III stage-age individuals, SE = standard error, DE% = percentage of the deviance explained by each model.

Intercept Estimate SE p value DE%
log(area) -0.530 0.108 0.007 < 0.01 1.98
Density -0.008 3.031 0.148 < 0.01 3.20
III% -7.154 9.407 0.270 < 0.01 9.19

Population genetics

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 (NA) was 3.3 ± 0.5, the observed heterozygosity (HO) was 0.25 ± 0.05, whereas expected heterozygosity (HE) 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 FIS), with the exception of MAS, where the FIS value was slightly negative.

Table 4.

Genetic diversity parameters for the four Italian sub-populations of C. laurifolius. N = sample size; NA = number of alleles; NE = effective number of alleles; HO = observed heterozygosity; HE = unbiased expected heterozygosity; FIS = inbreeding coefficient; * = p < 0.05.

FOR (N = 75) C5 (N = 62) D (N = 42) MAS (N = 10)
Cislau1 3 1.03 0.03 0.03 -0.01 3 1.14 0.06 0.13 0.55* 1 1.00 0.00 0.00 0.00 1 1.0 0.00 0.00 -0.00
Cislau7 4 2.65 0.45 0.63 0.28* 2 1.83 0.40 0.46 0.11 3 2.98 0.52 0.67 0.21* 3 2.38 0.60 0.61 -0.03
Cislau11 3 1.10 0.10 0.09 -0.04 1 1.00 0.00 0.00 0.00 2 1.13 0.12 0.12 -0.06 2 1.34 0.30 0.27 -0.18
Cislau12 8 2.38 0.26 0.64 0.59* 8 4.32 0.65 0.78 0.16* 6 2.90 0.36 0.66 0.45* 4 1.70 0.50 0.44 -0.20
Cislau14 3 1.14 0.10 0.12 0.17 1 1.00 0.00 0.00 0.00 5 1.66 0.29 0.40 0.27* 3 1.87 0.30 0.49 0.35*
Mean 4.2 1.73 0.19 0.30 0.20* 3 1.86 0.22 0.27 0.27* 3.4 1.93 0.26 0.37 0.18* 2.6 1.6 0.34 0.32 -0.02

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 FST 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 FST was significantly different from zero, but very low (FST = 0.050; 95% CI 0.022-0.083), suggesting strong inter-subpopulation gene flow. Concerning this result, the mean value of NM (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.


Demography, reproductive fitness and genetic structure

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 (HO = 0.63; HE = 0.65) and for long-living perennials (HO = 0.63; HE = 0.68) (Nybom 2004), these results suggest that genetic diversity is low within the population as a whole (mean HO = 0.25 ± 0.06; mean HE = 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 FIS (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 FIS value of this small sub-population did not differ from Hardy-Weinberg expectations (Table 4). Although the FIS 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 FIS 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 (FST 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 FIS and population size (Honnay and Jacquemyn 2007). The absence of homozygotes for rare alleles could bias downwards the FIS value in small sub-populations (Young et al. 1999), whereas FIS 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.

Conservation approaches

As population size, breeding system and genetic structure of C. laurifolius in Italy were completely unknown, the results of the present study provided relevant new knowledge, crucial for designing a programme for species management and conservation. Despite habitat fragmentation seeming to have no effect on the reproductive fitness, it is argued that this species in Italy could be affected by an ongoing process of population size reduction, linked to inbreeding depression, loss of genetic variation and fixation of deleterious alleles. All these factors play a role in reducing the adaptive potential of a population (Delmas et al. 2014; Theodorou and Couvet 2015). This local reduction and impoverishment is congruent with a general range contraction and fragmentation, already highlighted in this species by Fernández-Mazuecos and Vargas (2010). Without active conservation actions and also taking into account low genetic diversity (Frankham et al. 2002), the Italian population of C. laurifolius will face increasing risks of extinction.

Gene flow amongst sub-populations may partially compensate for losses of genetic diversity. This could reduce the mating between relatives, avoiding the increase of homozygosity and inbreeding depression (Hardner and Potts 1995).

To cope with habitat fragmentation, often due to canopy closure, the following in situ conservation actions are needed: a) coppicing, to reduce competition and to provide adequate light intensity for the seedling growth (Jacquemyn et al. 2008) and b) filling the spatial gaps amongst the sub-populations as much as possible, by means of targeted translocations, to contribute in maintaining gene flow. For translocation activities, there is no particular need to prefer any sub-population as a source of material, given the general homogeneity of the genetic asset. Concerning ex situ conservation, seeds from each sub-population are available at Pisa Germplasm Bank (Bedini and Carta 2010; Hay and Probert 2013). Studies about seed viability and ecology (Baskin and Baskin 2014) and phenology of radicle emergence (Carta et al. 2014) in this species will provide basic knowledge to establish a nursery of immature and virginile plants, to be translocated in situ. This may alleviate the lack of generation turnover in the natural population. Finally, a population genetics study using the same markers, but applied to populations collected all across the distribution range, may provide further insights for the conservation of this species.


This work was funded by the “Progetto di Ricerca di Ateneo” (PRA) of the University of Pisa, under grant number PRA_2016_1.


Paola Furio, Marta Sfingi, Junior Lacerda and Romario Tabosa are deeply acknowledged for their help in field and lab activities.


  • Abeli T, Gentili R, Rossi G, Bedini G, Foggi B (2009) Can the IUCN criteria be effectively applied to peripheral isolated plant populations? Biodiversity and Conservation 18: 3877–3890.
  • Ågren J (1996) Population size, pollinator limitation and seed set in the self-incompatible herb Lythrum salicaria. Ecology 77: 1779–1790.
  • Aguilar R, Quesada M, Ashworth L, Herrerias-Diego Y, Lobo J (2008) Genetic consequences of habitat fragmentation in plant populations: susceptible signals in plant traits and methodological approaches. Molecular Ecology 17: 5177–5188.
  • Aizen MA, Feinsinger P (1994) Forest fragmentation, pollination, and plant reproduction in a Chaco dry forest, Argentina. Ecology 75(2): 330–351.
  • Aizen MA, Harder LD (2007) Expanding the limits of the pollen-limitation concept: effects of pollen quantity and quality. Ecology 88: 271–281.
  • Albaladejo RG, Sebastiani F, González-Martı́nez SC, González-Varo JP, Vendramin GG, Aparicio A (2010) Isolation of microsatellite markers for the common Mediterranean shrub Myrtus communis (Myrtaceae). American Journal of Botany 11: e23–e25.
  • Alonso I, Luis E, Tárrega R (1992) First phases of regeneration of Cistus laurifolius and Cistus ladanifer after burning and cutting in experimental plots. International Journal of Wildland Fire 2(1): 7–14.
  • Ashman TL, Knight TM, Steets JA, Amarasekare P, Burd M, Campbell DR, Dudash MR, Johnston MO, Mazer SJ, Mitchell RJ, Morgan MT, Wilson WG (2004) Pollen limitation of plant reproduction: ecological and evolutionary causes and consequences. Ecology 85: 2408–2421.
  • Astuti G, Roma-Marzio F, D’Antraccoli M, Gennai M, Villani M, Peruzzi L (2015) Schede per una Lista Rossa della Flora vascolare e crittogamica: Cistus laurifolius L. subsp. laurifolius. Informatore Botanico Italiano 47(2): 245–289.
  • Atiqur BR, Khan ML, Dayanandan S (2016) Genetic Structure and Diversity of Natural and Domesticated Populations of Citrus medica L. in the Eastern Himalayan Region of Northeast India. Ecology and Evolution 6: 3898–3911.
  • Baskin CC, Baskin JM (2014) Seeds. Ecology, biogeography and evolution of dormancy and germination. Academic Press, San Diego.
  • Brys R, Jacquemyn H, Endels P, Hermy M, De Blust G (2003) The relationship between reproductive success and demographic structure in remnant populations of Primula veris. Acta Oecologica 24: 247–253.
  • Carta A, Probert R, Moretti M, Peruzzi L, Bedini G (2014) Seed dormancy and germination in three Crocus ser. Verni species (Iridaceae): implications for evolution of dormancy within the genus. Plant Biology 16: 1065–1074.
  • Carta A, Savio L, Bedini G, Peruzzi L, Fisogni A, Galloni M (2016a) All in an afternoon: mixed breeding system in one-day lasting flowers of Hypericum elodes L. (Hypericaceae). Plant Biosystems 150(5): 1001–1009.
  • Carta A, Probert R, Puglia G, Peruzzi L, Bedini G (2016b) Local climate explains degree of seed dormancy in Hypericum elodes L. (Hypericaceae). Plant Biology 18: 76–82.
  • Chapman WE (1986) Age versus stage: an analysis of age- and size-specific mortality and reproduction in a population of Laminaria longicruris Pyl. Journal of Experimental Marine Biology and Ecology 85: 225–234.
  • Chapuis MP, Estoup A (2007) Microsatellite null alleles and estimation of population differentiation. Molecular Biology and Evolution 24(3): 621–631.
  • Colling G, Matthies D (2006) Effects of habitat deterioration on population dynamics and extinction risk of an endangered, long-lived perennial herb Scorzonera humilis. Journal of Ecology 94: 959–972.
  • D’Antraccoli M, Aiello F, Mirabile M, Perfetti A, Logli F, Bedini G, Peruzzi L (2016) Contro il vortice dell’estinzione: il caso della traslocazione di Symphytum tanaicense Steven (Boraginaceae) nel Parco Regionale di Migliarino-San Rossore-Massaciuccoli. Codice Armonico 2016: 145–154.
  • Dansereau PM (1939) Monographie du genre Cistus. Boissiera 4: 1–90.
  • Delmas C, Cheptou PO, Escarvage N, Pornon A (2014) High lifetime inbreeding depression counteracts reproductive assurance benefit of selfing in a mass flowering shrub. BMC Evolutionary Biology 14: 243.
  • Dempster AP, Laird NM, Rubin DB (1977) Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society: Series B 39: 1–38.
  • Dimopoulos P, Raus T, Bergmeier E, Constantinidis T, Iatrou G, Kokkini S, Strid A, Tzanoudakis D (2013) Vascular plants of Greece: An annotated checklist. Botanic Garden and Botanical Museum Berlin-Dahlem and Hellenic Botanical Society.
  • Duminil J, Fineschi S, Hampe A, Jordano P, Salvini D, Vendramin GG, Petit RJ (2007) Can population genetic structure be predicted from life-history traits? American Naturalist 169: 662–672.
  • Galloni M, Podda L, Vivarelli D, Cristofolini G (2007) Pollen presentation, pollen-ovule ratios, and other reproductive traits in Mediterranean legumes (Fam. Fabaceae - Subfam. Faboideae). Plant Systematics and Evolution 266: 147–164.
  • García D, Zamora R, Hódar JA, Gómez JM (1999) Age structure of Juniperus communis L. in the Iberian peninsula: conservation of remnant populations in Mediterranean mountains. Biological Conservation 87: 215–220.
  • Gargano D, Bellusci F, Pellegrino G, Palermo AM, Bernardo L, Musacchio A (2009) The conservation perspectives and value of small and isolated plant populations: preliminary clues from Gentianella crispata (Gentianaceae) at the western boundary of its range. Annales Botanici Fennici 46: 115–124.
  • Gargano D, Fenu G, Medagli P, Sciandrello S, Bernardo L (2007) The status of Sarcopoterium spinosum (Rosaceae) at the western periphery of its range: ecological constraints led to conservation concerns. Israel Journal of Plant Science 55: 1–13.
  • Gatsuk LE, Smirnova OV, Vorontzova LI, Zaugolnova LB, Ztiukova LA (1980) Age states of plants of various growth forms: a review. Journal of Ecology 68: 675–696.
  • Glémin S, Bazin E, Charlesworth D (2006) Impact of mating systems on patterns of sequence polymorphism in flowering plants. Proceedings of the Royal Society of London B 273: 3011–3019.
  • Grossoni P, Venturi E (2009) Il cisto laurino di Santa Brigida. Tipografia Poggiali S.A.S., Rufina.
  • Hamrick J, Godt M (1989) Allozyme diversity in plant species. In: Brown A, Clegg MT, Kahler AL, Weir BS (Eds) Plant population genetics, breeding and genetic resources.Sinauer, Sunderland, 43–63.
  • Hamrick JL, Godt MJV (1996) Effects of life history traits on genetic diversity on plant species. Philosophical transactions of the Royal Society of London 351: 1291–1298.
  • Hardner CM, Potts BM (1995) Inbreeding depression and changes in variation after selfing Eucalyptus globulus ssp. globulus. Silvae Genetica 44: 46–54.
  • Harper JL (1977) Population Biology of Plants. Academic Press, New York.
  • Herrera J (1985) Nectar secretion patterns in southern Spanish Mediterranean scrublands. Israel Journal of Botany 34(1): 47–58.
  • Ismail SA, Ghazoul J, Shaanker R, Kushalappa CG, Kettle J (2012) Does long-distance pollen dispersal preclude inbreeding in tropical trees? Fragmentation genetics of Dysoxylum malabaricum in an agro-forest landscape. Molecular Ecology 21: 5484–5496.
  • Jacquemyn H, Brys R, Honnay O, Hermy M (2008) Effects of coppicing on demographic structure, fruit and seed set in Orchis mascula. Basic and Applied Ecology 9: 392–400.
  • Jacquemyn H, Van Rossum F, Brys R, Endels P, Hermy M, Triest L, De Blust G (2003) Effects of agricultural land use and fragmentation on genetics, demography and population persistence of the rare Primula vulgaris, and implications for conservation. Belgian Journal of Botany 136: 5–22.
  • Jakobsson M, Rosenberg NA (2007) clumpp: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. Bioinformatics 23: 1801–1806.
  • Kettle CJ, Hollingsworth PM, Jaffré T, Moran B, Ennos RA (2007) Identifying the early genetic consequences of habitat degradation in a highly threatened tropical conifer, Araucaria nemorosa Laubenfels. Molecular Ecology 16: 3581–3591.
  • Knight TM, Steets JA, Vamosi JC, Mazer SJ, Burd M, Campbell DR, Dudash MR, Johnston MO, Mitchell RJ, Ashman TL (2005) Pollen limitation of plant reproduction: pattern and process. Annual Review of Ecology, Evolution, and Systematics 36: 467–497.
  • Kunin WE (1997) Population size and density effects in pollination: pollinator foraging and plant reproductive success in experimental array of Brassica kaber. Journal of Ecology 85(2): 225–234.
  • Moro MJ, Domingo F, Escarré A (1996) Organic matter and nitrogen cycles in a pine afforested catchment with a shrub layer of Adenocarpus decorticans and Cistus laurifolius in south-eastern Spain. Annals of Botany 78: 675–685.
  • Piry S, Luikart G, Cornuet JM (1999) BOTTLENECK: a computer program for detecting recent reductions in the effective population size using allele frequency data. Journal of Heredity 90: 502–503.
  • Pritchard JK, Stephens M, Donnelly P (2000) Inference of population structure using multilocus genotype data. Genetics 155: 945–959.
  • QGIG (2016) QGIS Development Team. QGIS Geographic Information System. Open Source Geospatial Foundation.
  • R Core Team (2016) R: A language and environment for statistical computing. R Foundation for Statistical Computing. Vienna, Austria. Available from:
  • Rodriguez-Perez J (2005) Breeding system, flower visitors and seedling survival of two endangered species of Helianthemum (Cistaceae). Annals of Botany 95(7): 1229–1236.
  • Roma-Marzio F, D’Antraccoli M, Astuti G, Peruzzi L (2016b) Riscoperta della stazione storica di Cistus laurifolius subsp. laurifolius (Cistaceae) in località Masseto (Pontassieve, Firenze). Atti della Società Toscana di Scienze Naturali, Memorie, Serie B 122(2015): 29–32.
  • Rosas F, Quesada M, Lobo JA, Sork VL (2011) Effects of habitat fragmentation on pollen flow and genetic diversity of the endangered tropical tree Swietenia humilis (Meliaceae). Biological Conservation 144: 3082–3088.
  • Rozen S, Skaletsky HJ (2000) Primer 3 on the WWW for general users and for biologist programmers. In: Krawetz S, Misener S (Eds) Bioinformatics methods and protocols: Methods in molecular biology. Humana Press, Totowa.–386
  • Severns P (2003) Inbreeding and small population size reduce seed set in a threatened and fragmented plant species, Lupinus sulphureus ssp. kincaidii (Fabaceae). Biological Conservation 110(2): 221–229.
  • Silvertown J, Franco M, Pisanty I, Mendoza A (1993) Comparative plant demography - relative importance of life-cycle components to the finite rate of increase in woody and herbaceous perennials. Journal of Ecology 81: 465–476.
  • Thanos CA, Georgidou K, Kadis C, Pantazi C (1992) Cistaceae: a plant family with hard seeds. Israel Journal of Plant Sciences 41: 251–263.
  • Warburg EF (1968) Cistus L. In: Tutin TG, Heywood VH, Burges NA, Moore DM, Valentine DH, Walters SM, Webb DA (Eds) (1976) Flora Europaea, vol. 2. Cambridge University Press.
  • Wright SI, Ness RW, Foxe JP, Barrett SCH (2008) Genomic consequences of outcrossing and selfing in plants. International Journal of Plant Sciences 169: 105–118.
  • Xia J, Lu J, Wang ZX, Hao BB, Wang HB, Liu GH (2013) Pollen limitation and Allee effect related to population size and sex ratio in the endangered Ottelia acuminata (Hydrocharitaceae): implications for conservation and reintroduction. Plant Biology 15: 376–383.
  • Zhang ZY, Chen YY, Li DZ (2005) Detection of low genetic variation in a critically endangered Chinese pine, Pinus squamata, using RAPD and ISSR markers. Biochemical Genetics 43: 239–249.

Supplementary materials

Supplementary material 1 

Characteristics of the 6 polymorphic microsatellites markers developed for Cistus laurifolius

Giovanni Astuti, Francesco Roma-Marzio, Marco D’Antraccoli, Gianni Bedini, Angelino Carta, Federico Sebastiani, Piero Bruschi, Lorenzo Peruzzi

Data type: molecular data

This dataset is made available under the Open Database License ( The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited.
Download file (13.24 kb)
Supplementary material 2 

Genetic diversity parameters for the 4 stands of C. laurifolius

Giovanni Astuti, Francesco Roma-Marzio, Marco D’Antraccoli, Gianni Bedini, Angelino Carta, Federico Sebastiani, Piero Bruschi, Lorenzo Peruzzi

Data type: molecular data

This dataset is made available under the Open Database License ( The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited.
Download file (14.04 kb)
Supplementary material 3 

Results of STRUCTURE analysis

Giovanni Astuti, Francesco Roma-Marzio, Marco D’Antraccoli, Gianni Bedini, Angelino Carta, Federico Sebastiani, Piero Bruschi, Lorenzo Peruzzi

Data type: molecular data

This dataset is made available under the Open Database License ( The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited.
Download file (29.59 kb)
login to comment