Genetic diversity and differentiation in narrow versus widespread taxa of Helianthemum (Cistaceae) in a hotspot: The role of geographic range, habitat, and reproductive traits

Abstract Unraveling the relationships between ecological, functional traits and genetic diversity of narrow endemic plants provide opportunities for understanding how evolutionary processes operate over local spatial scales and ultimately how diversity is created and maintained. To explore these aspects in Sierra Nevada, the core of the Mediterranean Betic‐Rifean hotspot, we have analyzed nuclear DNA microsatellite diversity and a set of biological and environmental factors (physicochemical soil parameters, floral traits, and community composition) in two strictly endemic taxa from dolomite outcrops of Sierra Nevada (Helianthemum pannosum and H. apenninum subsp. estevei) and two congeneric widespread taxa (H. cinereum subsp. rotundifolium and H. apenninum subsp. apenninum) that further belong to two different lineages (subgenera) of Helianthemum. We obtained rather unexpected results contrasting with the theory: (a) The narrow endemic taxa showed higher values of genetic diversity as well as higher average values of pollen production per flower and pollen‐to‐ovule ratio than their widespread relatives; and (b) the two taxa of subg. Helianthemum, with larger corollas, approach herkogamy and higher pollen production than the two taxa of subg. Plectolobum, displayed lower genetic diversity and higher values of inbreeding. Altogether, these results disclose how genetic diversity may be affected simultaneously by a large number of intrinsic and extrinsic factors, especially in Pleistocene glacial refugia in mountains where the spatial context harbors a great ecological heterogeneity. On the other hand, differences in mating system and the significant effect of the substrate profile, both being highly diverse in the genus Helianthemum, in the genetic variability illustrate about the importance of these two factors in the diversification and species differentiation of this paradigmatic genus in the Mediterranean and open the field to formulate and test new hypotheses of local adaptation, trait evolution, and habitat diversification.


| INTRODUC TI ON
It is well known that biological diversity is concentrated in particular regions around the world, so-called biodiversity hotspots.
Despite scanty information being available about the building up of biodiversity in such regions, attention is recently being paid to the evolutionary history of species assemblages and lineages by considering phylogenetic relationships for whole floristic datasets (e.g., Molina-Venegas, Aparicio, Pina, Valdés, & Arroyo, 2013), communities (Anacker & Harrison, 2012), or lineages radiated in particular hotspots (Valente et al., 2010). These approaches have permitted analysis of the effect of environmental factors in promoting differentiation and diversity and thus the building up of such hotspots (see Anacker, Whittall, Goldberg, & Harrison, 2011 for the California hotspot).
The Mediterranean region is one of the largest biodiversity hotspots in the world (Myers, Mittermeier, Mittermeier, da Fonseca, & Kent, 2000), further composed by a number of subhotspots. One of them lies in the Betic-Rifean area (Médail & Quézel, 1997) and harbors a range of altitudes, lithologies, and climatic conditions that mostly reproduces those present in the whole Mediterranean Basin (Molina-Venegas et al., 2013). In addition to many endemic plants locally evolved, it also harbors many refuged taxa that avoided extinction during drier and colder periods in the Miocene and the Pleistocene, respectively (Médail & Diadema, 2009;Rodríguez-Sánchez & Arroyo, 2011). Within this area, Sierra Nevada is the core of the Betic-Rifean hotspot and shows one of the highest floristic diversities in the Iberian Peninsula and around the Mediterranean Basin (Médail & Diadema, 2009;Molina-Venegas et al., 2013). Very recently, its woody flora has been subjected to an evolutionary account of biodiversity using barcoding techniques and the construction of megaphylogenies, and it has been possible to determine the effects of elevation and substrate mosaicism on the phylogenetic structure of its species assemblages (Molina-Venegas, Aparicio, Lavergne, & Arroyo, 2015;Simón-Porcar et al., 2018). Whereas these studies have allowed detecting the prevailing evolutionary patterns, further insight requires scaling down to particular lineages diversified in the region, and even to species and population level of taxa subjected to systematic discussion, to disentangle how diversity in this hotspot is created and maintained, and the drivers involved. In this regard, unraveling the relationships between ecological, functional traits and genetic diversity of narrow endemic plants provide opportunities for understanding how evolutionary processes operate over these local spatial scales (Arafeh et al., 2002;Simón-Porcar et al., 2018).
Large reviews of genetic data have provided valuable insights into the patterns of genetic variation and their correlates, showing their dependence on many factors related to geographic range, effective population size, life form, mating and breeding system, environmental changes and biogeographical events, demographic processes and history of populations, polyploidization, hybridization, and natural selection (Ellegren & Galtier, 2016). Nevertheless, it is also necessary to underline the existence of a strong phylogenetic signal in the levels of genetic diversity between related species due to the conservativeness of functional traits in evolution, particularly in pollen and seed dispersal mechanisms (Duminil et al., 2007). Overall, it could be expected that these factors act synergistically in biodiversity hotspots as drivers of genetic diversity and diversification.
This study focuses on four taxa of Helianthemum from Sierra Nevada (southern Spain) that are two pairs of relatives with disparate distribution area (two local endemics vs. two widespread) and soil preferences (two dolomite specialists vs. two soil generalist), which moreover represent two different lineages (subgenera) within Helianthemum (Aparicio et al., 2017). The aim of this study was to unravel the effect of environmental, reproductive, and phylogenetic factors on the patterns of genetic diversity and differentiation in our precise case study. Specifically, for each taxa, we have gathered data about soil characteristics, community composition, and floral traits (as subrogates of the breeding system) to assess the effect of these factors on their genetic diversity and spatial genetic distribution of nuclear DNA microsatellite variation. Beyond deriving possible implications for the implementation of conservation action plans for the two stenochorous taxa in the Sierra Nevada National Park involved in this study, the assessment of the microevolutionary forces that drive the species divergence and differentiation in Helianthemum in a context of recent radiation can further shed light on why many species in this genus are prone to endemism.

| Study area
Sierra Nevada is located at the core of the Betic-Rifean region in SE Spain (37.07°N, 3.18°W) occupying an area of ca. 2,100 km 2 , the altitude ranging from 250 m a.s.l. to the highest peak in the Iberian Peninsula, the Mulhacén, at 3,482 m a.s.l. The climate is Mediterranean, characterized by cold winters and hot summers with a marked summer drought (July to August). The annual average temperature decreases in altitude from 12 to 16°C below 1,500 m to 0°C above 3,000 m a.s.l., and the annual average precipitation is about 600 mm, ranging from less than 250 mm in the lowest parts of the mountain range to more than 700 mm in summit areas. The number of vascular-plant taxa recorded in Sierra Nevada is 2,353 (Lorite, 2016), about 12% restricted to the Betic mountains (Blanca et al., 2002) and c. 80 local endemics (Lorite, Navarro, & Valle, 2007). This mountain range is currently considered one of the most important biodiversity hotspots in the Mediterranean region (Médail & Diadema, 2009)

| Study taxa and sampling
The four study taxa are chamaephytic shrubs belonging to the genus Helianthemum (Cistaceae; Figure 1) (Table 1). Nevertheless, it is necessary to stress that these WTs can never be found cohabiting with the former congeneric ETs.
Full or partial self-incompatibility has been reported in the few perennial species of Helianthemum whose breeding system have been studied (Tébar, Gil, & Llorens, 1997;Rodríguez-Pérez, 2005 who considered H. cinereum to be self-compatible) and, although the pollination biology of most species remains unknown, flowers are visited by generalist insects, mostly bees and beetles (Proctor, 1956;Rodríguez-Pérez, 2005;Agulló et al., 2015;personal observations).
In Sierra Nevada, we selected two areas differing in soil type  (Table 1), and to avoid biases in the measures of genetic variation and structure we sampled all the four taxa in a similar area (i.e., 132 and 197 km 2 ) and interindividual distance (cf. Cole, 2003). Whenever possible, depending on orographic circumstances, we aimed to keep the mean interindividual distance of the collected the same day coinciding with the beginning of blooming period (Blanca et al., 2002;López-González, 1993). We also col-

| DNA isolation and microsatellite genotyping
We ground 30 mg of dehydrated leaf tissue with a Retsch MM200 shaker mill and isolated total genomic DNA from each specimen with

| Soil and environmental sampling
We also collected 125 soil samples at 10-20 cm depth close to each of the 140 sampled plants (except for plants occurring less than 5 m apart). These soil samples were left to dry for 24 hr at room temperature and then passed through a 2-mm sieve before determining the following physicochemical parameters: texture (percentage of sand, silt and clay), pH (at 25°C 1:5), electrical conductivity (μS/cm, at 25°C 1:5), percentage of carbonates and organic matter, content of macronutrients (mg/kg of N-Kjeldahl and P-Olsen), micronutrients (mg/ kg of Cu, Fe, Mn, Zn, and Bo), and assimilable cations (meq/100 g of Ca, Mg, K, and Na). Soil analyses were performed by "Laboratorio Agrama SL" (Seville, Spain).
We also recorded in the field the community neighborhood and composition of each individual sampled plant by recording (a) the number of conspecifics in a 5-m radius buffer, (b) the distance to the nearest conspecific plant, and (c) the percentage of shrub coverage also in a 5-m radius buffer following a visual scale with four categories: 1 = 1%-25%, 2 = 26%-50%, 3 = 51%-75%, and 4 = 76%-100%.
Note. Sample size (N), number of polymorphic loci (Nloci), number of alleles (Na), effective number of alleles (Ne), unbiased expected heterozygosity (He), and inbreeding coefficient (Fis). Differences for genetic diversity parameters between widespread and endemic taxa within each subgenus assessed through two-sample nonparametric studentized permutation tests. Significant values highlighted in bold.

| Floral traits
For each taxon, we averaged the number of stamens and ovules per ovary after dissecting one floral bud from 5 to 21 individual plants under a dissecting microscope totaling 59 flower buds (Table 3). For each sample, we also estimated the number of pollen grains per anther by squashing one anther randomly selected on a slide containing Isoton II (Beckman Coulter, Fullerton, CA) and counting the number of pollen grains contained in aliquots of 500 μl using a particle counter (Coulter Multisizer 3, Beckman Coulter). Pollen production per flower was estimated as the number of pollen grains per anther (mean value of three aliquots) multiplied by the number of anthers in the flower; then, we computed the pollen-to-ovule (P/O) ratios accordingly (Cruden, 1977). We also dissected 1-3 recently- negative values when the stigma is below the anthers (i.e., reverse herkogamy).

| Genetic diversity and inference of historical events on genetic structure
We estimated the following standard parameters of genetic diversity: observed number of alleles (Na), effective number of alleles (Ne), unbiased Nei's gene diversity (He), and inbreeding coefficient (Fis) with GenAlEx 6.5 (Peakall & Smouse, 2012  To assess the number of independent genetic and evolutionary entities, we inferred the genetic structure by two methodologically contrasting approaches: the Bayesian clustering procedure of Structure (Hubisz, Falush, Stephens, & Pritchard, 2009;Pritchard, Stephens, & Donnelly, 2000) and the discriminant analysis of principal components (DAPC: Jombart, Devillard, & Balloux, 2010). Both methods allow to identify the optimal number of clusters (K) in the data set and to assign simultaneously the sampled individuals to each of the inferred clusters. Nonetheless, each of these methods has limitations that can affect the validity of their results. The Bayesian clustering methods infer subtle population structure by minimizing linkage disequilibrium and departures from Hardy-Weinberg equilibrium within each inferred cluster (Pritchard et al., 2000), which are often difficult to verify and can restrict their applicability. In fact, the assumption of panmixia does not necessarily hold in these taxa, which can exhibit a mixed-mating system in natural populations (see Discussion). Alternatively, Principal Components Analysis (PCA)-based methods construct low-dimensional projections of the data that maximally retain the variance-covariance structure among the sample genotypes. However, while these low-dimensional projections allow for straightforward visualization of the underlying population structure, it is not always straightforward to derive and interpret estimates for global ancestry of sample individuals from their projection coordinates (Novembre & Stephens, 2008). We restricted these analyses to the set of five microsatellite loci common to all four taxa (see Results) in order to avoid potential effects of the nonrandom distribution of missing data. Since the main aim of our study is to detect possible admixture between individuals of the four studied taxa, the number of potential K clusters assessed ranged from 1 (assuming a single panmictic entity) to 4 (assuming every taxon formed its own cluster). The Bayesian clustering was con- Pophelper (Francis, 2017). The DAPC analysis was conducted with the R package adegenet v.2.0.1 (Jombart, 2008). The optimal number of genetic clusters (K) was estimated with the function find.cluster and identified as the one with the lowest BIC value.
We used the program Bottleneck 1.2.02 (Cornuet & Luikart, 1996) to detect recent demographical changes in the effective population size of the study taxa. We employed Wilcoxon's tests to detect heterozygosity excess in comparison with simulated values under mutation-drift equilibrium with two mutational models of microsatellite variation: the infinite allele model (IAM) and the two-phase mutational model (TPM), the latter allowing for 10% of single-step changes (Cristescu, Sherwin, Handasyde, Cahill, & Cooper, 2010).
The variation rate was set to 12, as recommended by Piry, Luikart, and Cornuet (1999) for microsatellite markers.
To assess significant relationships among environmental variables and genetic diversity for each taxa, we constructed a genetic distance (among individuals) matrix with GenAlEx 6.5 that was subsequently subjected to multiple regression on distance matrices

| Floral traits
We assessed differences among taxa in the petal length, stigmaanther separation, and the P/O ratio through general linear mixed models (GLMMs) with the package lme4 (Bates, Maechler, Bolker, & Walker, 2015) in R. We set individuals as random factor to account for measures accomplished on different flowers of the same individual, log-transformed P/O ratio before model adjustment, and validated models by checking the absence of significant patterns in the residuals of the fitted models. We accomplished significant differences among pairs of taxa through sequential Bonferroni correction.

| Genetic diversity and inferred historical events
Genetic diversity values for the individual microsatellite loci amplified in the four Helianthemum taxa of this study are in Supporting Information
Helianthemum apenninum subsp. apenninum showed a significant and positive autocorrelation in the pairwise genetic coefficient at the first distance class (0-1 km), followed by a decline with significant but negative spatial autocorrelation coefficients above 20 km.
Helianthemum pannosum showed a significant and positive autocorrelation in the pairwise genetic coefficient at the first distance class (0-1 km) and at 2-4 km. No significant spatial signal was detected in any distance class for the other two taxa. Results of multiple regressions on distance matrices are shown in Table 2. Genetic distance showed marginally significant correlation with soil parameters (all soil variables standardized) in H. apenninum subsp. estevei (pvalue = 0.054). The neighborhood-composition distance matrix and the geographic distance matrix were uncorrelated with genetic distance in all cases.

21
58.2 ± 2.2 (21) 6.1 ± 0.3 (21) 16,760 ± 1,636 (21) 2,840 ± 379 (21) 4.31 ± 0.94 (21) 3.60 ± 0.79 (21) −1.40 ± 0.14 (21) Hp 6.33 ± 2.00 (10) 10 31.3 ± 4.6 (5) 5.5 ± 0.2 (5) 17,569 ± 7,209 (5) 2,957 ± 1,322 (5) 5.04 ± 1.59 (10) 4.27 ± 1.35 (10) −0.77 ± 0.18 (10) Note. Sample size (number of flowers) in brackets. . The boxplot displays the smallest and largest values as well as the first quartile, the median, and the third quartile. Outliers are indicated with dots. Different letters indicate statistical significant differences among taxa after a sequential Bonferroni correction on regression coefficients from GLMM models mating system rather than geographic distribution is driving the differences in genetic variability between the two pairs of relatives in this study and the existence of a (albeit marginally) significant effect of the substrate profile in the genetic variability of the dolomite specialist H. apenninum subsp. estevei may be further indicating the important role of these factors on the diversification and species differentiation in Helianthemum, also showing that proneness to endemism in the genus is associated with ecological specialization, given the apparent distinct distribution and habitat types of many endemics (López-González, 1992). Altogether, from a conservation biology perspective our results underline that management and conservation action plans relying only on simple proxies for genetic diversity such as rarity, geographic range, successful sexual reproduction, or outbreeding may be masking the actual genetic patterns.

| Genetic diversity and historical events
In this study, we have found that genetic diversity of the ETs (H. pannosum and H. apenninum subsp. estevei) was at least equal or higher than in the WTs (H. cinereum and H. apenninum subsp. apenninum; Table 1). It is usually assumed that rare species tend to have low levels of genetic variability due to the small population size (e.g., Arafeh et al., 2002;Edwards, Lindsay, Bailey, & Lance, 2014;Turchetto et al., 2016). However, small population size coupled with substantial genetic diversity is not unusual in plants and has been reported for several rare species (Barrett & Kohn, 1991;Ellstrand & Elam, 1993;Lutz, Schneller, & Holderegger, 2000). Therefore, range extension on its own is not a reliable predictor of genetic diversity (Premoli, Souto, Allnutt, & Newton, 2001) and we must invoke alternative explanations to explain the unexpected high levels of genetic diversity found in the ETs.
Historical changes in population size and long-term isolation have often been proposed as possible explanations for such a lack of correlation between population size and genetic diversity (e.g., Landergott, Holderegger, Kozlowski, & Schneller, 2001). For example, the high genetic diversity found in small populations of some alpine species (also in the study area of Sierra Nevada) such as Gentiana alpina, Kernera saxatilis, and Papaver alpinum has been linked to shifts in their distribution ranges and vicariance processes occurred during the Pleistocene (Kropf, Comes, & Kadereit, 2006).
Alternatively, the high levels of genetic diversity that we have found in the ETs could be due to the existence of gene flow between the related taxa (Smith & Voung-Pham, 1996) and to recent diversification or phylogenetic constraints. The genetic distinctiveness of H. pannosum and H. cinereum detected by Structure and DAPC (see Figure 3) led us to discard contemporary hybridization events between the two species of subg. Plectolobum. However, the fact that individuals from the two subspecies of H. apenninum were considered as a single cluster by Structure or integrated into two not welldefined clusters by DAPC could indicate a recent divergence with persistence of gene flow between both subspecies. This is expected because the diversification of subg. Helianthemum did not likely occur until the Pleistocene (Aparicio et al., 2017) and the reproductive barriers between their members may be still labile, particularly at infraspecific levels (Nieto-Feliner, 2014). Although the taxonomic splitting of species subjected to systematic discussion such as, for example, H. apenninun, H. nummularium, H. marifolium, H. cinereum, or H. oelandicum into species and subspecies has been widely criticized because diagnostic characters are not stable but correlated with ecological conditions in most of the cases (Soubani et al., 2015;Volkova et al., 2016;Widén, 2015Widén, , 2018 Marrero-Gómez, Bañares, & Sosa, 2015;Santana-López, 2015). All these species have been considered relicts of ancient larger populations and wider distribution, but further studies specifically addressing these issues are required.

| Genetic diversity and environmental data
Helianthemum apenninum subsp. apenninum and H. pannosum showed a significant and positive spatial autocorrelation in the first distance class (0-1 km; Figure 4), an isolation-by-distance pattern usually interpreted to be consequence of limited dispersal, in agreement with the prevalent gravity dispersal mechanism reported for some species of Helianthemum (González-Pérez et al., 2013;Tébar et al., 1997). However, the absence of the isolation-by-distance result can indicate the existence of active mechanisms of adaptation to the stressful conditions since dolomitic soils drain much more efficiently promoting strong xericity (Mota et al., 2008;Salmerón-Sánchez et al., 2014). Despite the usual view regarding endemic taxa to be local specialists with reduced genetic variation (e.g., Thompson, 2005), numerous studies have shown that adaptation to stressing soil conditions can promote diversification within lineages eventually leading to increased levels of diversity at different spatial scales (Molina-Venegas et al., 2015;Rajakaruna et al., 2003), which is congruent with the higher levels of genetic diversity retrieved for the dolomite specialist H. apenninum subsp. estevei.
The prevalence of endemic species is outstanding in some Mediterranean taxonomic groups (Thompson, 2005), as is also no-

| Genetic diversity and floral traits
Perianth size, stigma-anther separation, P/O ratio, or life form are biological traits usually considered to be indicative of mating and breeding systems (Barrett, 2003;Cruden, 1977), also associated with range size (Lavergne, Thompson, Garnier, & Debussche, 2004).
For instance, narrowly distributed species may exhibit reproductive traits prone to reduce outcrossing compared to their widespread relatives (e.g., fewer and smaller flowers, less stigma-anther separation, and lower pollen-ovule ratios) as a strategy for persistence after local adaptation (Lavergne et al., 2004 showing that variation in floral traits is common in Helianthemum, even within populations (Aragón & Escudero, 2008;Martín-Hernanz et al., unpublished), and the prevalence of mixed-mating systems in these species (Aragón & Escudero, 2008;Rodríguez-Pérez, 2005).
Moreover, since H. pannosum and H. apenninum subsp. estevei also displayed a higher number of alleles (Na), effective number of alleles (Ne), and unbiased expected heterozygosity (He), it seems that mixed-mating systems are promoting survival and contributing to the maintenance, or even to the increase, of the genetic diversity of these two stenochorous taxa. This striking result has also been documented in the Cretan endemic Cyclamen creticum, whose genetic variability is not lower than in its widespread relative C. repandum, being the levels of genetic diversity more influenced by the mating system rather than by the amplitude of the geographic distribution (Affre & Thompson, 1997).

| Conservation implications
We have found that in Helianthemum, a recently evolved lineage where many species are currently just differentiating (Aparicio et al., 2017), genetic diversity is actually being driven by many intrinsic and extrinsic factors that act synergistically in promoting diversification in Sierra Nevada. This finding stresses the essential role of robust and dated phylogenies embracing the species under study to interpret patterns and to develop action plans for locally evolved of refuged species in biodiversity hotspots. Nevertheless, the core question in conservation genetics about how does genetic variability compare among rare and their widespread related species still persists (Gitzendanner & Soltis, 2000;Turchetto et al., 2016).
The two strictly dolomitic endemics in Sierra Nevada here studied (H. pannosum and H. apenninum subsp. estevei) have been evaluated as Vulnerable by the Spanish Red List of Endangered Plants (Moreno, 2008) due to their restricted geographic distribution range, the low number of populations, and the impact by grazing animals, albeit the rarity of these plants is probably due to natural causes related to the ecological specificity and the scarcity and the natural fragmentation of their habitats (Blanca et al., 2002). In our study, we have found that facultative xenogamy and mechanisms of adaptation to the stressful conditions imposed by dolomite soils can stand for the persistence on the study taxa (see Table 1 and Table 3); however, if we take into account the global warming scenario these populations, at the verge of alpine vegetation, can be at risk of rapid extinction because the Mediterranean summits are even at higher risk compared to other European mountains (Parmesan, 2006;Pauli et al., 2012). Therefore, H. pannosum and H. apenninum subsp. estevei require ongoing monitoring of their conservation status and it would be advisable the development of appropriate conservation programs to integrate both in situ and ex situ conservation techniques (Volis & Blecher, 2010). One key ex situ strategy for biodiversity conservation is the implementation of germplasm collections which have to be designed aiming to capture as much as genetic diversity as possible, especially regarding the rarest alleles (Caujapé-Castells & Pedrola- Monfort, 2004). The kind of precise genetic information at spatial scale that we provide in this study could be easily included in the implementation of germplasm collection programs by the environmental managers to maximize their cost-effectiveness and to ensure the conservation of endangered species.

ACK N OWLED G M ENTS
We thank the authorities of the Sierra Nevada National Park and the Junta de Andalucía (Andalusian Regional Government) for granting permission to collect samples even of species under protection.
We also thank the Centro de Investigación Tecnología e Innovación

DATA ACCE SS I B I LIT Y
Microsatellite data can be found at Dryad Digital Repository: doi.