When ecological marginality is not geographically peripheral: exploring genetic predictions of the centre-periphery hypothesis in the endemic plant Lilium pomponium

Study species, data collection and sampling

Lilium pomponium L. (Liliaceae) is a perennial herb (geophyte) endemic to the Maritime and Ligurian Alps—a region that represents a major regional hotspot of endemism and plant diversity in the Mediterranean Basin (Thompson, 2020). It grows in rocky grassland and shrubland from 100 to 2,000 m of altitude in an area of roughly 7.000 km² (Noble & Diadema, 2011). The species has a self-incompatible outcrossing breeding system, with a poor capacity for selfing (Casazza et al., 2018). Species occurrences were obtained from field surveys and the Conservatoire Botanique National Méditerranéen SILENE data base (http://www.silene.eu/index.php?cont=accueil) and LiBiOss (Regione Liguria; http://www.cartografiarl.regione.liguria.it/Biodiv/Biodiv.aspx). Overall, a final dataset of 881 occurrence records was used to assess the global distribution of the species and to determine central and peripheral sites and those that may be environmentally marginal. A total of 197 individuals were sampled from 20 populations that cover the entire geographical and altitudinal range of the species (Fig. 1A and Table 1). Field collection was authorized by the Conservatoire Botanique National Méditerranéen de Porquerolles and the Regione Liguria (Decree n. 859 of 22/08/2018).

DNA extraction and genotyping

Genomic DNA was isolated from fresh leaves according to Kobayashi et al. (1998). Frozen leaves samples were ground in a mixer 150 mill “TissueLyser” (Quiagen-Retsch, Hilden, Germany). Total DNA was extracted using NucleoSpin Plant II 151 Kit (Macherey & Nagel, Germany). DNA concentrations were measured using a photometer (Biophotometer, Eppendorf, Germany). The AFLP method was performed as describe by Vos et al. (1995) using laboratory equipment previously employed by Vidaller (2018).

AFLP genotyping was based on the protocol described by Vos et al. (1995). A total of 100 ng of DNA were digested using the restriction enzymes Eco RI and Tru 9I (Fisher Scientific, France) for 3 h at 37 °C and then for 3h at 65 °C in a total volume of 25 μl (15 μL + 10μl of DNA). Digestion products were ligated to 0.5 μL Eco and 25 μL Mse adaptors for 3 h at 37 °C and treated with T4 DNA Ligase and 0.1 μL of 100 mM ATP to a final volume of 25 μL (5 μL + 20 μl of restriction products). Ligation products were diluted eight times and pre-selective PCR amplification was performed using EcoR1+A, Mse+C primers and Taq DNA polymerase in a 44.5 μL volume. The profile of pre-amplification thermocycle was 94 °C for 2 min, followed by 20 cycles at 94 °C for 45 s, 56 °C for 45 s, 72 °C for 1 min and 72 °C for 10 min. For the selective amplification, three primer combinations were chosen for PCR: ASII: EcoR1-AGG/MseI-CGG, ASIII: EcoR1-AGC/MseI-CAG, ASVII: EcoR1-AGC/MseI-CTG dyed with 6-FAM fluorescence at 5′ Eco end (Eurofins Genomics, Ebersberg, Germany).

One Hundred times diluted pre-amplification products were used to perform selective amplification in a final volume of 20 μL (15 μL + 5 μL of diluted pre-amplification products). For the selective amplification thermocycle we used 94 °C for 2 min, 10 cycles of 94 °C for 30 s, 65 °C for 30 s (step −1 °C per cycle), 72 °C for 1 min, followed by 22 cycles at 94 °C for 30 s, 56 °C for 30 s, 72 °C for 1 min, 72 °C for 5 min and 4 °C for 2h. The fragment length produced by the amplification was separated and quantified by electrophoresis using an ABI 3730xl DNA analyzer (Applied Biosystems, Foster City, California, USA) with GS600 LIZ size marker. Peaks were scored in Peak Scanner V 1.0 (Applied Biosystems, Foster City, California, USA) as present (1), absent (0), or “no data” (NA). We utilized Raw Geno 2.0 (Arrigo et al., 2009) to select the fragments longer than 100 bp and smaller than 200 bp. Maximum binning between peaks was set at 1.75 and minimum at 1.5. The error rate was calculated as mismatches ratio over matches in 32 replicated individuals (Bonin et al., 2004).

Estimation of geographical and environmental center-periphery gradient

We used different measures to describe the geographical and environmental center-periphery gradients. To identify geographically central and peripheral populations we calculated the geographical center of distribution by averaging the coordinates of the occurrences and the distributional limits by calculating the minimum convex hull polygon that included all the known occurrences of L. pomponium. Following this, geographical peripheral populations were identified by calculating for each occurrence its distance from the geographical center, its shortest distance to the edge of the convex hull and the ratio between these two distance measures.

To assess environmental marginality nineteen bioclimatic variables representative of the period 1979–2013 were downloaded from the Chelsa climate database website (Karger et al., 2017a, 2017b; www.chelsa-climate.org) at 30-s (c. 1 km) spatial resolution. We extracted values of bioclimatic variables from cells where the species occurs and defined the climatic space occupied by the species on the basis of the first two axes of a principal component analysis (PCA) using all cell values. Environmental marginality was assessed by calculating the Euclidean distance from the centroid of the climatic space and the Mahalanobis distance from the mean of the climatic space.

Data analysis

Genetic diversity was measured for each population by calculating Nei’s diversity index (h) and Shannon’s information index (I) using POPGENE v. 1.32 (Yeh, Yang & Boyle, 1999). Because the various measures of population differentiation quantify different aspects of population structure, we estimated genetic differentiation among populations by using three different measures: the fixation index F_ST, estimated by using ARLEQUIN 3.11 (Excoffier & Lischer, 2010) and Rho (Ronfort et al., 1998) and Jost’s D (Jost, 2008) pairwise distance matrix using GenoDive 3.0 (Meirmans, 2020). To estimate the genetic differentiation among-population, we calculated the average of pairwise F_ST, Rho and D values for each population with all other populations. We calculated Pearson product moment correlations between geographical and environmental center-periphery distance measures and between genetic parameters and the center-periphery distance measures.

Geographical and environmental gradients

The first two axes of PCA of climatic data explained 79.71% of the overall climatic variance (29.03% and 50.68%, respectively in Fig. S1). The first PCA axis is positively related to temperature-related parameters and negatively related to precipitation during the warm and dry season. The second axis is positively related to precipitation during the wet and cool seasons and negatively related to temperature seasonality indices. Sites that are marginal in climatic space occurred throughout the entire distributional range, as witnessed by the occurrence of sites that are far from the climatic centroid both near the geographical center and in northern and southern geographically peripheral sites (dark red areas in Figs. 1C and 1E; Table 1).

The distance from the geographical center was higher in eastern and western peripheral sites than in northern and southern peripheral sites (dark red areas in Fig. 1B). As a result, sites farthest from the distributional edges occurred in a long and thin belt close to the center of distributional range (dark red areas in Fig. 1D). Likewise, areas with a low ratio between distance from the edge and the center of distribution occurred mainly at the eastern and western geographical periphery (light red areas in Fig. 1F). In line with these results, the distance from the climatic centroid was not correlated (r = 0.09 and r = 0.05 according to Euclidean and Mahalanobis distance, respectively) with distance from the geographical center (Figs. 1, 2A and 2D; Table 1) and was only poorly negatively correlated with distance from distributional edge and the edge/centre distance ratio (Figs. 1, 2B, 2C, 2E and 2F; Table 1).

Genetic diversity

The three AFLP primer pairs produced 134 polymorphic fragments for 197 individuals (Table 2). The AFLP error rate based on 32 replicates was 2.02%. Genetic diversity indices per population ranged from 0.11 to 0.30 with a mean of 0.22 (±0.04) and from 0.19 to 0.43 with a mean of 0.34 (±0.06) according to Nei’s diversity index and Shannon information index, respectively (Table 3). Average F_ST between each population and all remaining populations ranged from 0.034 to 0.152, the average Jost’s D ranged from 0.032 to 0.108 and the average Rho ranged from 0.008 to 0.321 (Table 3). Populations with low and high genetic diversity were scattered throughout the distributional range (Figs. 3A and 3B). In particular, populations with lowest values occurred near the geographical center (i.e., P05), as well as near both the southern (i.e., P01) and the northern (i.e., P20) periphery of the distribution. Similarly, populations with high values occurred both near the geographical center (i.e., P07 and P11) and near the periphery (i.e., P14 and P16).

Populations with low and high F_ST and Jost’s D values were scattered throughout the distributional range (Figs. 3C and 3D). In contrast, populations with high Rho values occurred mainly near the geographical periphery (i.e., P01, P13, P16 and P17 in Fig. 3E) while populations with low Rho values occurred near the geographical center (i.e., P04, P06 and P07). In line with this pattern, intra-population genetic diversity was not significantly correlated with geographical distance measures (Table 4). Two measures of average inter-population differentiation (i.e., F_ST and Jost’s D) were not correlated with measures of geographical marginality while Rho index was significantly correlated with measures of geographical marginality, increasing toward the geographical periphery (Table 5). Genetic diversity indices were significantly correlated with both distances from climatic centroid (Table 4) while none of the genetic differentiation indices per population were significantly correlated with environmental distances (Table 5).

Lack of concordance among geographical and environmental gradients

An implicit assumption of the CPH is a symmetrical, monotonic deterioration of environmental condition from the center of the distribution toward the periphery (Brown, 1984). However, in mountainous environments, such as the Maritime and Ligurian Alps where L. pomponium occurs, environmental factors change over very short distances because of high topographic complexity (Körner, 2003; Casazza et al., 2008; Thompson, 2020), imposing environmentally marginal conditions in scattered and diffuse parts of a given species’ range. Moreover, topographic heterogeneity may decouple site climatic conditions from regional climatic values, reducing our capacity to detect causal factors related to demographic process and environmental variation (Hannah et al., 2014; Patsiou et al., 2014).

In L. pomponium climatically marginal conditions were found both at the periphery of distribution and close to the geographical center, mainly on valley floors (dark red areas in Figs. 1C and 1E; Table 1). This scattered distribution of central and marginal environmental conditions in different parts of the distributional range of this species are probably favored by the succession of high and low elevation areas. In fact, in the study region, tectonic, glacial and fluvial processes have shaped the deeply incised valleys that are characteristic of the southern Alps (Sanchez et al., 2010). As a result, the lack of correlation between geographical and environmental distances is not surprising (Fig. 2).

In addition, the distribution range of L. pomponium is largely west-east orientated because the north-south extension is constrained by the sharp altitudinal gradient of the mountain chain. As a result, the east and west peripheral populations are further away from the center than south and north peripheral populations (Figs. 1B–1F). For this reason, harsh conditions in the latter are closer to the geographical center than harsh conditions in the former peripheral populations (dark red areas in Figs. 1C and 1E). This asymmetry may explain why correlations are higher between environmental gradients and distance from the edge of distribution than between environmental gradients and the distance from the geographical center (Figs. 2B, 2D and 2F). In L. pomponium, the weak association between geographical and environmental distances is thus closely affected by the complex topography of the region that causes a discrete and strongly asymmetrical distribution of environmental conditions across the distributional range. Correlations between geographical and ecological gradients thus may not occur in areas of high localized spatial heterogeneity. Indeed, the occurrence of central environmental conditions at the geographical margins (white areas in Figs. 1C and 1E) suggests that in L. pomponium range limits may be due to other factors, such as a gradient in habitat availability, rather than a gradient in habitat quality. This supports the idea that local environmental conditions (such as in microrefugia) explain the occurrence of a patchy distribution (Patsiou et al., 2014).

Genetic structure and the CPH

In general, intra-population genetic diversity is expected to be associated with inter-population differentiation because stochastic genetic drift due to change in population size causes both a decline in genetic diversity within populations and divergence among populations (Eckert, Samis & Lougheed, 2008). In fact, peripheral populations having a low genetic diversity usually have also a high genetic differentiation (Pironon et al., 2017). However, our results do not provide support for the CPH prediction of a reduction in within-population genetic diversity towards the periphery but do provide evidence for an increase in inter-population genetic differentiation in peripheral populations relative to those in the center.

For L. pomponium, although Rho measure of genetic differentiation illustrate high genetic differentiation near the geographical periphery, suggesting a decrease in gene flow among peripheral populations, we found no evidence for a decline in genetic diversity in peripheral populations (Fig. 3; Tables 3 and 4). The difference between Rho and both F_ST and Jost’s D is probably because of Rho is independent from the selfing rate, while the other measures increase when self-fertilization occurs within populations (Ronfort et al., 1998). This may explain why some small-sized populations near the geographical center but far from environmental center (i.e., P04, P05 and P08 in Table 3) have high F_ST and Jost’D values and low Rho values. In fact, the high number of flowers per inflorescence detected in these small populations may favor geitonogamous pollination among flowers on the same plant (Macrì et al., 2021). In small sized populations of facultative-autogamous species, like L. pomponium (Casazza et al., 2018), geitonogamy may assure the production of a small, but regular number of seeds (Roberts et al., 2014), but at the cost in terms of reduced outcrossing (Lloyd, 1992; Harder & Barrett, 1995).

The high genetic differentiation among peripheral populations suggested by Rho index may be the result of a decrease in the availability of suitable habitat and therefore a more (naturally) fragmented distribution of populations rather than the result of a decline in habitat quality and population size (Pironon et al., 2017). In fact, the reduction in density of suitable habitat may reduce population connectivity both by pollen and seeds (Hargreaves & Eckert, 2014; Young, Boyle & Brown, 1996), favoring isolation of populations and resulting in higher inter-population differentiation. Moreover, in L. pomponium population connectivity may be further reduced by a decrease in suitable habitat due to (i) shrubland or forest expansion in mid-elevation of high elevation sites (Carlson et al., 2014) and/or (ii) urbanization in peripheral coastal lowland areas (Noble & Diadema, 2011). Given that gene flow in perennial plants mainly occurs by pollen transfer (Levin & Kerster, 1974; Ennos, 1994; Tarayre et al., 1997), pollen limitation detected in the southernmost populations of L. pomponium may further drive the genetic differentiation among these populations (Macrì et al., 2021).

In general, diversity is expected to decline towards the geographical periphery, regardless of environmental conditions (Eckert, Samis & Lougheed, 2008; Pironon et al., 2017). However, in L. pomponium genetic diversity does not decrease with the distance from the geographical center but does decline with distance from climatic centroid (Table 4). The lack of association between genetic diversity and distance from the geographical periphery is indicative of historical range stability, and of the existence of patchwork of many local refugia, the so-called refugia within refugia hypothesis. The occurrence of widespread local refugia characterizes the Mediterranean region, explaining the persistence of species during the Pleistocene glaciation cycles (Médail & Diadema, 2009; Thompson, 2020) and their ability to rapidly recolonize new locations at higher altitude or latitude at the end of glaciation (Villellas et al., 2014; Ferreira et al., 2015).

In species survived in few glacial refugia, patterns of genetic variation have been strongly shaped by past climate-driven range dynamics, and consequently a decrease of within-population genetic diversity is expected towards recently colonized areas at high altitude or latitude because of the stochasticity of colonization processes (Hampe & Petit, 2005; Thompson, 2020). However, contrary to this expectation, in L. pomponium we detected the highly genetically diverse populations at high latitude and elevation and the less genetically diverse populations at low altitude or latitude edge. This finding supports the hypothesis of a past range expansion and a recent contraction in relation with the complex phylogeographical history of the Maritime and Ligurian Alps (Médail & Diadema, 2009). A further analysis aimed precisely at reconstructing past range dynamics would be most interesting here. In contrast to what we have observed in L. pomponium, high genetic diversity near the geographical center has been recorded in other Mediterranean plant species in which a central–marginal pattern may be associated with a history of expansion from a center of diversity (Pouget et al., 2013). Elsewhere, the pattern we observed in L. pomponium has been previously detected in 40 species in which environmental distance more closely predicts genetic diversity than geographical distance (Lira-Noriega & Manthey, 2014).

The reduction of neutral genetic diversity in environmentally marginal conditions may be due to a cascade effect of environmental impacts on population dynamics (Lira-Noriega & Manthey, 2014; Pironon et al., 2017). The quality of the local habitat (i.e., distance from the environmental center) may affect demographic processes (Jaquiéry et al., 2008), ultimately resulting in a reduction in genetic variation due to selection and/or genetic drift. Nevertheless, in L. pomponium some small populations have relatively high genetic diversity (e.g., P03, P14 and P16 in Table 3). The perennial nature of L. pomponium (more than 10 years in L. martagon; Lundqvist, 1991) may reduce the effect of genetic drift (Ellstrand & Elam, 1993; Nybom & Bartish, 2000). In addition, natural selection may directly play an important role in shaping neutral genetic diversity (Gillespie, 2000; Leffler et al., 2012; Corbett-Detig, Hartl & Sackton, 2015). In large populations that incur strong environmental filtering, natural selection may directly affect levels of neutral variation and favor variants that are linked to beneficial mutations (Gillespie, 2000) and/or purge those linked to deleterious mutations (Kaplan, Hudson & Langley, 1989), the so-called genetic draft (Gillespie, 2000). In L. pomponium genetic differentiation may thus be affected by habitat availability that reduces population connectivity while genetic diversity may be affected by habitat quality that affects population dynamics.