Introduction

The endemic and disjunct distribution of closely related species in the Mediterranean flora is often thought to have arisen from range fragmentation (Favarger and Contandriopoulos 1961; Thompson et al. 2005; Mayol et al. 2012), mostly due to the spatial isolation of tectonic microplates that occurred during the Tertiary (Cardona and Contandriopoulos 1979; Verlaque et al. 1991; Médail and Quézel 1997). Present-day distribution patterns and speciation dynamics have also been heavily affected by successive palaeogeographical events (climate oscillations and eustatic variations), which fostered dispersal, recolonisation and secondary contacts (e.g. hybridisation) (Médail and Diadema 2009). This results in a more heterogeneous evolutionary history than the one expected from simple range fragmentation and allopatric speciation (Taberlet et al. 1998; Thompson 2005; Weiss and Ferrand 2007; Médail and Diadema 2009). The taxonomic units generated by these processes are not always well separated; the Mediterranean flora is rich in recently diverged and poorly defined taxa with scattered patterns (e.g. Celtis L., De Castro and Maugeri 2006; Anchusa L., Coppi et al. 2008; Nigella L., Bittkau and Comes 2009; Senecio L., Molins et al. 2009; Genista L., De Castro et al. 2015) rich in congeneric endemic/widespread pairs (Lavergne et al. 2004; Matesanz et al. 2009). Conservation and decision making for these kinds of group are therefore complex, but essential to guarantee that diversity be maintained together with the main features of Mediterranean vegetation (Quézel 1985; Fréville et al. 1998; Gitzendanner and Soltis 2000). This observation can be easily adapted to other biodiversity hotspots (Myers et al. 2000).

One of the recently diverged groups of plants whose current distribution has been influenced by range fragmentation is Asperula sect. Cynanchicae (DC.) Boiss. (Rubiaceae) (Gargiulo et al. 2015). Among the species of this section, the present study is focused on A. crassifolia L., a narrow endemic to southern Italy (Campania region). This species is one of the several endemic members of the section Cynanchicae and it was previously included in the informal “ser. Palaeomediterraneae”; the series grouped schizoendemic species (i.e. derived from range fragmentation and speciation of a more widespread ancestor; Favarger and Contandriopoulos 1961) distributed in different sites (mostly islands) of the eastern, central and western Mediterranean (Gutermann and Ehrendorfer 2000). However, the monophyly of this group has not been confirmed because of the concurrence of several biological processes affecting the phylogenetic boundaries within sect. Cynanchicae (Gargiulo et al. 2015).

As for many other allied Asperula endemics, the occurrence in stable habitats (steep rocky slopes) is suggestive of poor competitiveness and avoidance of human-affected sites (Polunin 1980; Lavergne et al. 2004; Thompson et al. 2005). Although a direct count for the species has not been carried out, A. crassifolia might be considered as an “extremely narrow endemic” (ENE, López-Pujol et al. 2013; Fernández-Mazuecos et al. 2014) with a very restricted number of populations and small census sizes (< 500). The inclusion of A. crassifolia in the IUCN Red List as a vulnerable species is recommended (Santangelo, pers. comm.)

In the present study, we assess genetic diversity and divergence of A. crassifolia populations, by employing plastid markers and nuclear microsatellites (nrSSRs), to investigate: (1) the relationships among the haplotypes detected in A. crassifolia and in other allied species of A. sect. Cynanchicae; (2) the population history, by assessing genetic variation within this species (and among its populations). Moreover, (3) we discuss the conservation implications of population diversity in A. crassifolia.

Materials and methods

Study species

Asperula crassifolia (Online Resource 1), recently typified by Peruzzi et al. (2013), is a caespitose perennial, with four leaves per node (leafy stipules are here counted as leaves) at the base of the stem, in pairs from middle of the stem upwards. Inflorescences are pyramidal with hermaphroditic flowers; corollas are hypocrateriform, yellowish or reddish, densely patent-hairy as the whole plant, with tube 2–3 times as long as lobes (Guadagno 1913; Ehrendorfer and Krendl 1976). Fruits are dry and indehiscent (two mericarps). Studies on the reproductive and spreading strategies of A. crassifolia are lacking, but self-incompatibility through mechanisms of dimorphic heterostyly is reported for the perennial members of the tribe Rubieae (Bahadur 1968; Anderson 1973; Bremer and Eriksson 1992; Tao and Ehrendorfer 2011). Asperula crassifolia is morphologically similar to its widespread and sympatric congener A. aristata subsp. aristata (e.g. Tenore 1811; Parlatore 1848), especially in the case of individuals with reduced hairiness (in the past indicated as A. commutata C.Presl; Roemer and Schultes 1818; Béguinot 1907; Fiori 1923). However, we comply with the traditional taxonomic treatment, which describes A. crassifolia as a separate species, especially on the basis of the following unique morphological traits: yellowish/reddish, hairy corollas and fleshy leaves (Ehrendorfer and Krendl 1976; Ehrendorfer 1982; Del Guacchio et al. 2016). A presumed close relationship with A. aristata subsp. calabra (Fiori) Del Guacchio & Caputo has been recently disproved (Del Guacchio et al. 2016). Limited morphological divergence, however, is common in A. sect. Cynanchicae (e.g. Ehrendorfer and Krendl 1976; Del Guacchio and Caputo 2005, 2013).

Asperula crassifolia grows on rocky calcareous substrates of the southern Italian island of Capri, in the neighbouring Sirenusae archipelago, and in a single locality of the Sorrentine Peninsula (Nerano). Another locality in inland Campania, i.e. Mount Fellino (668 m a.s.l., Moraldo and La Valva 1989), is not confirmed (Santangelo, pers. comm.) and appears doubtful for geobotanical and ecological reasons.

Sampling

The rarity of A. crassifolia prevented the collection of a conspicuous number of individuals. Non-destructive sampling was conducted during the flowering season of the species and one voucher specimen was collected per population (NAP); sampling sites are shown in Fig. 1. The island of Capri was explored in each of the sites reported in the literature (Boccone 1697; Tenore 1811; Guadagno 1922; Ricciardi 1996); at Nerano, the population was discovered in the 1990s (Caputo et al. 1989–1990) and it is restricted to a rocky slope. At Sirenusae, the only accessible island of the archipelago, Gallo Lungo, was explored.

Fig. 1
figure 1

Sampling sites with haplotype distribution for Asperula crassifolia populations. Asterisks indicate sampling sites at Capri (from left to right: Monte Solaro, Via Krupp and Pizzolungo-Via Tragara)

Full size image

Forty-eight individuals were sampled, representative of the populations of A. crassifolia (Table 1). The sampling strategy adopted was aimed at avoiding clonal individuals. Additionally, 10 historical herbarium specimens preserved in NAP and B herbaria (acronyms according to Thiers 2016) and collected in 1808–1933 from the island of Capri were included in the analysis (Online Resource 2). In order to assess the relationships between A. crassifolia haplotypes and those detected in closely related species, we included in the analysis the following members of A. sect. Cynanchicae: A. aristata subsp. aristata, A. aristata subsp. calabra, A. aristata subsp. oreophila (Briq.) Hayek (= A. aristata subsp. scabra Nyman), A. borbasiana (Korica) Korica, A. cynanchica subsp. cynanchica L., A. cynanchica subsp. neglecta (Guss.) Arcang., A. deficiens Viv., A. garganica Huter, Porta & Rigo ex Ehrend. & Krendl, A. gussonei Boiss., A. lutea Sm., A. naufraga Ehrend. & Gutermann, A. paui subsp. paui Font Quer, A. peloritana C.Brullo, Brullo, Giusso & Scuderi, A. pumila Moris, A. staliana subsp. diomedea Korica, Lausi and Ehrend., A. visianii Korica (see Online Resource 2). This selection was carried out according to the results of the phylogenetic investigation on A. sect. Cynanchicae, choosing some of the taxa most closely related to A. crassifolia and the geographically contiguous subspecies of the two most widespread species of the section, i.e. A. aristata and A. cynanchica (Gargiulo et al. 2015).

Table 1 Genetic diversity parameters computed on the microsatellite dataset for Asperula crassifolia populations
Full size table

Total genomic DNA was extracted after lysis in liquid nitrogen with CTAB method (Doyle and Doyle 1987); for recalcitrant material, the extraction was followed by purification using GENECLEAN® II Kit (MP Biomedicals).

Plastid markers

Several plastid markers (trnC-petN, psbA-trnH and trnQ-rps16 intergenic spacers, matK gene and rps16 intron) were screened for variability, by analysing the most divergent species of the dataset (following the phylogenetic outline of Gargiulo et al. 2015); eventually, rps16 intron and trnC-petN intergenic spacer were selected. For rps16 amplifications, we used the primers described by Oxelman et al. (1997) and an internal primer designed ad hoc for Asperula; primer pairs are listed in Online Resource 3. Polymerase chain amplifications were carried out in a final volume of 25 μL, 1 U DreamTaq DNA Polymerase (Fermentas, Thermo Scientific), 1X DreamTaq Buffer, 0.5–2.5 mM MgCl2, 0.2 mM dNTPs, 0.5 μM primers, DNA at various concentrations, water to the final volume and the following thermal programme: denaturation at 94 °C for 3 min, 30 cycles at 94 °C for 30 s, annealing at 55 °C for 30 s, extension for 45 s at 72 °C, followed by 5 min of final extension. PCR products were purified using PEG 8000 precipitation (http://www.mcdb.lsa.umich.edu/labs/olsen/files/PCR.pdf, last accessed January 2018). Sequencing was conducted using BigDye Terminator Cycle Sequencing Kit v3.1 (Applied Biosystems, Life Technologies), with the following conditions: 96 °C for 10 s, 50 °C for 5 s and 60 °C for 4 min, for 25 cycles. Sequence purification was conducted according to Di Maio and De Castro (2013). Sequences were electrophoresed on an ABI 3130 Genetic Analyzer (Applied Biosystems, Life Technologies). The resulting sequences have been submitted to GenBank (see Online Resource 2 for the accession numbers).

Plastid sequences were aligned using Muscle (Edgar 2004), as implemented in Mesquite v3.01 (Maddison and Maddison 2014); rps16 alignment required further manual editing, because of the presence of insertions and deletions. Two distinct matrices were generated: one with the sequences of the entire A. crassifolia sampling (A. crassifolia dataset) and the other including the detected A. crassifolia haplotypes plus the sequences of the other Asperula accessions (total dataset). Accessions in the total dataset were assigned to different haplotypes using DnaSP v5.10 (Librado and Rozas 2009). Relationships among A. crassifolia haplotypes were presented as a statistical parsimony network, generated using the program TCS v1.21 (Clement et al. 2000). Gaps were treated as fifth state, by converting each indel to a one-nucleotide state, with a 94%-connection limit, as it gave more resolution than the default value (95%). Asperula aristata subsp. aristata from Sorrentine Peninsula (i.e. DEI) was used as an outgroup.

The relationships among all the Asperula accessions were investigated in MrBayes v3.2.6 (Huelsenbeck and Ronquist 2001; Ronquist and Huelsenbeck 2003) using the online facility CIPRES Science Gateway (Miller et al. 2010). Asperula lutea was chosen as outgroup. The model GTR + I was selected as best-fitting according to jModeltest v3.7 (Posada and Crandall 1998), using the Akaike information criterion (Akaike 1974); indels were treated as missing data. Markov chain Monte Carlo sampling was conducted for 1 × 107 number of generations, with a sampling frequency of 1000 generations and a burn-in of 25%.

Plastid DNA haplotype diversity for each population (hS) and for the total range (hT), and two coefficients of gene differentiation (GST, NST) were estimated according to Pons and Petit (1996), as implemented in PermutCpSSR v2, with 1000 permutations.

Nuclear microsatellite genotyping

Nine loci were selected among the microsatellite markers developed for A. crassifolia: GA_1, GA_11, GA_17, GA_5A, GA_55A, GA_30C, GA_51C, GTT_1, GTT_ 22 (Gargiulo and De Castro 2015) and amplified in sampled populations; herbarium specimens were not genotyped as they could not be considered as discrete populations. Amplifications were carried out in a final volume of 25 μL, with 50 ng genomic DNA, 4 pmol reverse primer, 4 pmol labelled-M13 fluorescent dye (6-FAM, PET, VIC, NED), 1 pmol M13-tailed forward primer, 1X Taq Buffer B (Solis BioDyne), 0.8 mM dNTP mix (Promega), 2 mM MgCl2, 1U FIREPol DNA Polymerase (Solis BioDyne), on a 2700 thermal cycler (Applied Biosystems). Thermal programme is described in Gargiulo and De Castro (2015). One millilitre of diluted PCR product was added to 12 μL of formamide/GeneTrace 500 LIZ size standard (Carolina Biosystems) and run on the ABI 3130 Genetic Analyzer (Applied Biosystems). Fragments sizing was conducted in PeakScanner v1.0 (Applied Biosystems).

Presence of null alleles and genotyping errors (large allele dropout or stuttering) was assessed with MicroChecker v2.2.3 (van Oosterhout et al. 2004) using 1000 randomisations; anomalies were rectified using the Brookfield 2 algorithm (Brookfield 1996). We assessed the occurrence of clonal genotypes in the R (R Development Core Team 2008) package Poppr v2.8 (Kamvar et al. 2014). In order to assess the effect of null alleles and to correct the bias induced by their presence on the FST estimation, we used the Excluding Null Alleles (ENA) method as implemented in the program FreeNA (Chapuis and Estoup 2007), using 10,000 bootstrap replications. The software estimates the null allele frequencies and the global and pairwise FST values (Weir 1996). We conducted this assessment on the uncorrected dataset and after Brookfield correction, to verify the congruence of all the pairwise FST outcomes.

Departures from Hardy–Weinberg equilibrium (HWE), linkage disequilibrium (LD) and FIS values were evaluated with the exact tests implemented in GenePop v4.2.2 (Raymond and Rousset 1995). Significance level was adjusted using sequential Bonferroni correction for multiple comparisons (Rice 1989). GenAlEx v6.5 (Peakall and Smouse 2006) was used to calculate genetic diversity parameters: i.e. observed and (unbiased) expected heterozygosities (respectively HO, HE), percentage of polymorphic loci and total number of alleles. Allelic richness was calculated by using the rarefaction approach (Petit et al. 1998), as implemented in FSTAT v2.9.3.2 (Goudet 2002).

Analysis of molecular variance (AMOVA), as implemented in Arlequin v3.5 (Excoffier and Lischer 2010) Locus by Locus Amova using the number of different alleles—FST), was carried out to verify partitioning of total genetic variation at the infra- and interpopulation levels, considering population subdivision at Capri based on the locality of sampling (Fig. 1). Statistical significance was tested with 10,000 permutations.

Population structure was assessed by using the software MavericK v1.0 (Verity and Nichols 2016), using a procedure known as Thermodynamic Integration (TI), which is considered as more accurate in determining the real value of genetic clusters (K) (Verity and Nichols 2016). Both admixture and without-admixture models were tested, with K ranging from 1 to 8, 15 × 104 MCMC iterations, a burn-in of 5000 samples, repeating the MCMC procedure 5 times. TI parameters were: 50 rungs, 5000 burn-in and 55 × 103 thermodynamic samples. Bar plots were obtained by employing Distruct v1.1 (Rosenberg 2004). Results obtained in MavericK were compared to the ones obtained in Structure v2.3 (Pritchard et al. 2000; see “Results” section).

Demographic history of the populations, with particular emphasis on recent (i.e. last 2Ne–4Ne generations; Luikart et al. 1998a) reduction in effective population size was explored with Bottleneck v1.2 (Piry et al. 1999). Tests were performed with 10,000 replications, under the step-wise mutation model (SMM) and the two-phase model (TPM), which are believed to be appropriate for microsatellite evolution (Luikart and England 1999; Balloux and Lugon-Moulin 2002; Williamson-Natesan 2005). TPM model included a proportion of SMM set to 70% and the variance to 30% (multiple step mutation). Sign test and the one-tailed Wilcoxon signed rank were employed to verify significant deviations in observed heterozygosity, as predicted in the cases of recent bottleneck. Bottleneck method is based on the assumption that allele frequency is the result of equilibrium between mutation and genetic drift (Piry et al. 1999); consequently, we also tested the departure (Mode-shift) from the L-shaped distribution of allelic frequencies, which is typical of an equilibrium situation (Luikart et al. 1998b). One deme at Capri (Pizzolungo-Tragara) was excluded from the Bottleneck computations, as the number of individuals (n = 4) was not sufficient to perform the analysis.

Results

Plastid data

The A. crassifolia dataset consisted of 521 positions for the rps16 intron and 729 positions for trnC-petN, whereas for the total dataset, rps16 alignment resulted in 574 positions and trnC-petN in 685 positions (as a result of removing the initial portion because it was incomplete for some of the accessions). Total alignment is shown in Online Resource 4. The rps16 region turned out to be the most variable, consisting of different indel polymorphisms, each one considered as a different state in the TCS analysis (see Online Resource 5 for the polymorphic positions in A. crassifolia). On the basis of these polymorphisms, 17 plastid haplotypes were recognised for the Asperula accessions included in the analysis, four of them detected in A. crassifolia. TCS network is shown in Fig. 2. Haplotype A was the only one occurring at Nerano and Sirenusae; at Capri, it was present in a single accession, the herbarium specimen PAS2 collected in 1868 (see Online Resource 2). Noteworthily, haplotype B was the most represented in the historical herbarium specimens (8 out of 10 specimens), whereas it is less represented in the present-day collection at Capri (3 out of 19 specimens). Haplotype C was the most abundant, whereas haplotype D was private and present in a single individual at Capri (Online Resource 2).

Fig. 2
figure 2

TCS network, showing the relationship among the haplotypes of Asperula crassifolia. Asperula aristata subsp. aristata DEI was chosen as outgroup

Full size image

The Bayesian analysis on the total dataset showed that Haplotype C is clustered with A. cynanchica accessions (A. cynanchica subsp. cynanchica AN1, A. cynanchica subsp. neglecta) and the Adriatic taxa A. borbasiana and A. staliana subsp. diomedea in a strongly supported group (posterior probability—PP = 0.99), whereas the other haplotypes are shared with several other taxa (Figs. 3, 4).

Fig. 3
figure 3

Bayesian phylogram showing the relationships among haplotypes in Asperula crassifolia and other members of Asperula, based on plastid regions. Note that gaps were treated as missing data

Full size image
Fig. 4
figure 4

Map showing the distribution of the Asperula accessions analysed in the present study. Taxa sharing the same haplotype as Asperula crassifolia are marked with boxes indicating the haplotype name. ARI2A. aristata 2H; ARI8A. aristata 8H; BORA. borbasiana; CALA. calabra; CRAA. crassifolia; CYNA. cynanchica subsp. cynanchica AN1; DIOA. staliana subsp. diomedea; DEIA. aristata subsp. aristata DEI; DEFA. deficiens; GARA. garganica; GUSA. gussonei; LUTA. lutea; NAUA. naufraga; NEGA.cynanchica subsp. neglecta; ORE1A. aristata subsp. oreophila 131; ORETA. aristata subsp. oreophila TRM11; PAUA. paui; PUMA. pumila; PELA. peloritana; VISA. visianii. For explanations about the accessions, see Online Resource 2

Full size image

Genetic diversity indices are shown in Online Resource 6; hT and hS were respectively 0.70 and 0.20. The coefficients of differentiation measured over the three populations were GST = 0.705 and NST = 0.720.

Nuclear microsatellite data

Among the nine microsatellite loci, no evidence of large allele dropout or stuttering was detected; the presence of null alleles was suggested by the homozygote excess for the loci GA_17, GA_5A, GA_55A, GA_30C, differentially displayed in each population (Online Resource 7). Multilocus genotypes were variable in 47 out of 48 individuals, as two individuals at Nerano displayed the same genotype (although with a large amount of missing data). Null alleles did not affect significantly pairwise FST values (Online Resource 8).

Deviation from HWE was found for the locus GA_1 at Capri and Nerano, and for GA_17 and GA_5A at Nerano and Sirenusae (see Online Resource 9 for the Bonferroni correction); there was no evidence for linkage disequilibrium. FIS values significantly departed from random mating for almost all loci (Online Resource 10); with a mixture of negative and positive FIS values at Capri and predominantly positive values at Nerano and Sirenusae. Mean values of the genetic diversity indices for each population and total are summarised in Table 1 (details for each locus are shown in Online Resource 7). Globally, Nerano and Sirenusae populations exhibit less genetic diversity than Capri population in terms of heterozygosity. Only at Capri, HO was comparable to HE (respectively 0.468 and 0.459), whereas for Nerano and Sirenusae populations, HO was lower than expected (see Table 1). Allelic richness computed on the minimum sample size of one diploid individual was highest in Capri’s locations (values averaged over loci > 1.40, Table 1; detailed values for each locus are shown in Online Resource 11). AMOVA analysis showed moderate genetic differentiation among populations (21%), with most of the variation partitioned within individuals (67%; Online Resource 12).

Analysis of genetic structure as implemented in MavericK revealed a higher posterior probability associated with the without-admixture model, with the most likely K being equal to 2 (Fig. 5). Admixture levels were higher in Nerano, where individuals displayed a mixed proportion of the two genetic clusters. The most homogeneous populations in terms of genetic structure are Sirenusae and Capri-Via Krupp (Fig. 5). A comparison with the results obtained in the software Structure is shown in Online Resource 13.

Fig. 5
figure 5

Results of the analysis of genetic differentiation in Asperula crassifolia, as implemented in MavericK v1.0. a Posterior probabilities associated with the without-admixture model and the admixture model. b Posterior probabilities associated with the number of genetic clusters, K. c Bar plot for K = 2

Full size image

Bottleneck analysis did not indicate excess of heterozygosity associated with recent bottlenecks at Capri (considered as a single population) and at Nerano, either under SMM or TPM models. However, there was evidence of a shift in the normal L-shaped distribution at each locality at Capri and at Sirenusae; no evidence of bottleneck was found according to the other two tests (Online Resource 14).

Discussion

Genetic diversity and structure of the Asperula crassifolia populations

Both plastid and nuclear markers revealed a limited genetic diversity in Nerano and Sirenusae populations. One plastid haplotype was detected, and heterozygosity values computed on the nSSRs were lower than the values reported for narrow/endemic populations by Nybom (2004, i.e. HO = 0.52/0.32 and HE = 0.56/0.42, with N < 16). Moreover, inbreeding coefficients per locus were mainly positive. If we exclude the presence of clonal individuals, either the occurrence of null alleles, non-random mating or unrecognised internal structure producing an apparent heterozygote deficit (Wahlund effect), may be the causes of such values of FIS. Although null alleles are likely for GA_17 and GA_5A, as the homozygote excess seems to be locus-specific, the population-specific trend at Nerano and Sirenusae suggests that limited gene flow occurs among the populations. This is also supported by pairwise FST values (Online Resource 8).

At Capri, instead, three different haplotypes were detected (apart from the haplotype detected in a single herbarium specimen) and variation at nuclear microsatellites was expressed by higher heterozygosity values (Table 1) and comparable allelic richness, when alleles were normalised to the minimum sample size (Table 1; Online Resource 11). However, allelic richness may be underestimated because of the constraints imposed by the rarefaction method. The analysis of genetic differentiation revealed that the most homogeneous populations in terms of genetic clusters were Capri-Krupp and Sirenusae (Fig. 5), where almost all individuals have proportions of ancestry close to 100%. At Nerano, some individuals displayed a shared ancestry with individuals at Sirenusae; however, most of them were members of the main cluster detected in Capri, in line with the FST values (0.09 for Capri/Nerano vs. 0.21 Nerano/Sirenusae; Online Resource 8). These results were also in line with the AMOVA results, which defined a moderate partition of genetic variation among populations (ca. 20%; Online Resource 12).

The Bottleneck analysis showed contrasting results, in supporting a recent reduction in population size at Sirenusae and in the demes at Capri through the allele frequency distribution analysis, but not through the other tests (Online Resource 14). Moreover, the small population sizes at Solaro and Krupp and the fact that interbreeding among demes might occur in the island of Capri may have biased the outcome of the analysis, by determining spurious signals of bottlenecks (Wakeley 1999; Chikhi et al. 2010). At Sirenusae, on the other hand, a recent bottleneck may be supported by the observation of a particularly abundant population at the beginning of the last century (Guadagno 1913). We believe that Sirenusae is more threatened with genetic drift, as the observed levels of genetic variation are lower (Table 1, Online Resource 11). The evidence obtained from the nSSRs and, especially, from the plastid haplotypes suggests that genetic diversity at Nerano and Sirenusae is lower but not distinct from genetic diversity in Capri. However, our results do not clarify whether the lower genetic variation at Nerano and Sirenusae is either a consequence of a bottleneck or of a founder event from Capri.

Phylogeographic history and origin of Asperula crassifolia

Asperula crassifolia shares haplotypes with other members of the sect. Cynanchicae. In particular, Haplotype C is shared by the single island endemic A. borbasiana (Krk, Croatia) and Haplotype B is shared by various species of Asperula, including two individuals of A. aristata s.l. from southern Italy. The remaining two haplotypes are exclusive to A. crassifolia (Fig. 4; Online Resource 5). Haplotype sharing among different species may be explained as a consequence of (1) persistence of ancestral haplotypes (i.e. incomplete lineage sorting; Pamilo and Nei 1988; Maddison 1997; Avise 2000) or (2) hybridisation (and introgression) with the only sympatric taxon of the section, i.e. A. aristata subsp. aristata.

It is noteworthy that the two accessions of A. aristata subsp. aristata sampled in the Sorrentine Peninsula bear two distinct haplotypes, whereas the Haplotype B of A. crassifolia is identical to the haplotypes found in two individuals from other locations of southern Italy (subsp. aristata 2H and subsp. oreophila TRM11) (Fig. 4). However, our recent search for A. aristata at Capri, in the sites indicated under various synonyms by Tenore (1811), Guadagno (1922) and, even more recently, by Ricciardi (1996) has been unsuccessful. Consequently, it was not possible to confirm hybridisation between A. aristata and A. crassifolia. In other parts of A. aristata’s range, however, hybridisation and/or secondary contacts are suggested (Gutermann and Ehrendorfer 2000; Gargiulo et al. 2015), and the occurrence of individuals with intermediate features are reported (e.g. Tenore 1831, sub A. tomentosa). Nevertheless, we presume that ptDNA variation in A. crassifolia at Capri is also attributable to the persistence of ancestral haplotypes due to insufficient divergence time (Avise 2000). Hence, genetic variation in A. crassifolia is arguably part of the genetic variation of the ancestor of different species within sect. Cynanchicae. Population genetics theory predicts that, in small populations, polymorphisms will be lost more rapidly by stochastic events (i.e. genetic drift) and thus ancestral polymorphisms are more likely to occur in situation of recent and ongoing diversification (Ellstrand and Elam 1993). Recent diversification is hypothesised for other Mediterranean schizoendemic groups, such as Senecio sect. Senecio (Comes and Abbott 2001), Digitalis minor L. (Sales et al. 2001) and Crepis triasii (Mayol et al. 2012), where ancestral polymorphism persists by virtue of recent diversification.

Our data suggest incomplete lineage sorting related to recent and ongoing diversification in the genus Asperula, as previously reported; however, our results do not allow us to determine when A. crassifolia speciated. Capri Island and Sorrentine Peninsula were intensely shaped by eustatic and tectonic variations until the upper Pleistocene. During the middle Pleistocene (700,000–120,000 ya), a tectonic depression determined the insularity of Capri. Afterwards, an alternation of sea rises and regressions shaped the island profile and made repeatedly possible a connection via land bridges of Capri, Sirenusae, and Sorrentine Peninsula, up to 12,000 ya (Barattolo et al. 1992; Guida and Vallario 2003). Consequently, it is possible that the splitting of A. crassifolia occurred at any point during a period of isolation. Higher genetic diversity at Capri may suggest that the origin of the species occurred in the island. However, we cannot exclude that speciation occurred in Sorrentine Peninsula, where individuals and genetic variation were subsequently lost, according to the current diversity observed at Nerano and Sirenusae. Further population genetic research is necessary to clarify variation and relationships in the A. aristata/A. crassifolia species pair. However, we support the current recognition of A. crassifolia as a separate species, on the basis of morphological traits (Ehrendorfer and Krendl 1976) and of the exclusive haplotypes found in the present study.

Conclusions and conservation remarks

The Quaternary period is responsible for the genetic diversification and speciation of many endemic species (Hewitt 2000, 2004; Willis and Niklas 2004), including A. crassifolia (Gargiulo et al. 2015). Rarity in most of these species is mainly an innate condition, not descending from genetic erosion, but from local adaptation and species divergence in a peculiar ecological context (Rabinowitz et al. 1986; Gaston 1994). However, rarity triggers further rarefaction (O’Grady et al. 2004), especially in a human-affected and changing landscape. Despite the problems with taxonomic circumscription due to the “speciation in action” (Gargiulo et al. 2015), these taxa deserve protection as they are extremely valuable not only as part of the endemic/rare flora, but also as a present record of the intricate pattern of phenomena which lead to speciation in highly dynamic ecosystems. In A. crassifolia, we have observed small populations with some possible size reductions compared to past observations, and low genetic diversity in Nerano and Sirenusae. Although the number of individuals growing in inaccessible sites has not been estimated, it is quite improbable that population sizes will be sufficient to counterbalance the effects of genetic drift, inbreeding depression and stochastic oscillations of the habitat. These processes may be exacerbated by small population sizes. From a conservation perspective, human activity may contribute to further rarefaction of the species. Risk factors for the A. crassifolia are especially connected to the management systems against the natural erosion of the rocky slopes where the species grows (i.e. containing nets), and to the maintenance of tourist paths, especially at Capri and Nerano. In the Sirenusae archipelago, on the contrary, the habitat is relatively protected as not easily accessible; however, the population at Sirenusae was previously reported as more abundant than in Capri (Guadagno 1913, 1922). Given our findings, we recommend that the diversity in Capri be prioritised for conservation protection. Little is known about the actual number of individuals, which is presumably very low (R. Gargiulo, pers. obs.), but some effort should be made in order to detect individuals growing in inaccessible coastal slopes, both in Capri and in the Sirenusae archipelago. As the species is considered as Vulnerable according to IUCN (2012) [VU B1 ab (I, II, III, IV); B2 ab (I, II, III, IV), Santangelo, pers. comm.], ex situ preservation of the germplasm might be beneficial for eventual reintroductions, together with research aimed at investigating germination conditions and pollination mechanisms and constraints.