Influence of natural habitat fragmentation on the genetic structure of Canarian populations of Juniperus turbinata

Oceanic archipelagos provide an important platform from which to evaluate the effects of isolation and fragmentation on the genetic structure of species. As a result of oceanic isolation, such species usually show lower levels of genetic diversity and higher genetic differentiation than their mainland congeners. However, this is not necessarily the case for long distance dispersal species, whose genetic structure is not strictly defined by population isolation. We assessed the level and distribution of genetic diversity among Canarian populations of Juniperus turbinata in order to evaluate the influence of population isolation on its genetic structure. Using Amplified Fragment Length Polymorphism markers, we analyzed molecular diversity among 175 individuals from five populations occurring across the Canary Island and three Moroccan populations. Principal Coordinate Analysis, neighbor joining clustering, AMOVA and Bayesian-based analysis were applied to examine population structure. Despite the documented habitat loss and decline in Canarian populations, Amplified Fragment Length Polymorphism markers revealed levels of intra-population genetic diversity that were similar to those from mainland populations, and low levels of genetic differentiation. Bayesian analysis of population structure showed three main clusters, one comprising El Hierro population and a few individuals from several islands, a second cluster that grouped the remaining Canarian populations together, and a third cluster grouping Moroccan populations. Our results suggest that the main force driving the genetic structure of Canarian populations of J. turbinata is its capacity for long distance dispersal.


Introduction
Juniperus phoenicea group (family Cupressaceae) is an assortment of small monoecious or rarely dioecious Mediterranean trees whose distribution extends encompasses the Mediterranean basin and Macaronesian regions, from Madeira and Canary Islands in the west to Jordan and Saudi Arabia in the east, and is most prevalent throughout the Iberian Peninsula and North Africa (Quézel and Pesson 1980;Greuter et al. 1984;do Amaral Franco 1986;Farjon 2005).The species comprises two subspecies, J. phoenicea subsp.phoenicea and J. phoenicea subsp.turbinata (Guss.)Nyman, which were defined based on taxonomic (Farjon 2005), morphological (Mazur et al. 2010) and molecular studies (Adams et al. 2002;Boratyński et al. 2009;Dzialuk et al. 2011).Alternative interpretations define two independent species, J. phoenicea s. str.and J. turbinata Guss., with J. turbinata being discussed by several authors (see Farjon 2005), and its recognition as a separate being supported by morphological and biochemical studies (Lebreton and Pérez de Paz 2001;Mazur et al. 2016) and recent molecular studies (Adams and Schwarzbach 2013;Adams 2014;Adams et al. 2013Adams et al. , 2014)).As such, J. phoenicea ≡ J. phoenicea subsp.phoenicea is generally considered as being endemic to the Iberian Peninsula, Northwest Italy and the south of France, while J. turbinata ≡ J. phoenicea subsp.turbinata is widespread throughout the Mediterranean basin and Macaronesia.Furthermore, while the Macaronesian populations of J. turbinata were formerly described as J. canariensis Guyot & Mathou, they were later combined as a subspecies, J. turbinata subsp.canariensis (Guyot & Mathou) Rivas Mart., Wildpret & P. Pérez by Rivas-Martínez et al. (1993).However, its taxonomic status has recently been challenged based on several molecular markers (RAPDs, Adams et al. 2006;cpDNA, Adams et al. 2010;nrDNA and cpDNA, Adams et al. 2013a), leading to the conclusion that J. turbinata subsp.canariensis should not be supported at the subspecific level.
Regardless of whether Macaronesian populations represent a well-defined subspecies, they are scarce and threatened, and the islands where they occur are largely considered priority natural habitats by the European Directive 92/43, due to the particular richness of endemic species (Otto et al. 2012).At present, Canarian populations of J. turbinata are located in the five westernmost islands of the archipelago, and only occupy some 10% of the species potential distribution range due to widespread human disturbance (Fernández-Palacios et al. 2008;Del Arco et al. 2010).The taxon is reported as having been formerly widespread throughout the Canary Islands (including the eastern islands of Lanzarote and Fuerteventura), until their colonization, approximately six centuries ago.Populations of Tenerife, La Palma and Gran Canaria are usually small and frequently formed by isolated individuals, whereas the largest populations are located on La Gomera and El Hierro islands.Del Arco et al. (2010) calculated a distribution area of 3943 ha for J. turbinata across the Canary Islands, with 70% occurring in La Gomera.
Although there are several studies that point out the dynamics of Canarian populations of J. turbinata (Fernández-Palacios et al. 2008;Otto et al. 2012), there is a lack of knowledge on the genetic structure of these populations.As such, research into the amount and distribution of genetic diversity is required to a better understanding of associated ecological processes such as island colonization, gene flow among populations or historical bottlenecks, and to enhance and inform existing conservation measures aimed at protecting the populations (Fernández-Palacios et al 2008;Otto et al. 2010).
It is well known that oceanic islands provide ample opportunity to test hypotheses on the ecological factors affecting species' genetic structure.These may include colonization patterns, species origin/diversity, island biogeography and population differentiation (Francisco-Ortega et al. 2000;Emerson 2002;Stuessy et al. 2014).On oceanic islands, taxa are usually somewhat differentiated from those of mainland populations and they tend to display lower levels of genetic diversity because each island can be considered as an isolated subpopulation, and the potential area of population distribution is restricted, so population sizes are usually limited, making populations more prone to genetic risks, such as genetic drift, loss of genetic diversity or even extinction by stochastic factors or human disturbance (Ellstrand and Elam 1993;Franckham 1997;Whittaker et al. 2008).In addition, genetic differentiation is higher between oceanic islands than within them and populations mostly exhibit isolation-by-distance patterns of differentiation (Franckham 1997;Whittaker and González-Palacios 2007).However, some studies have also shown similar or even higher levels of genetic variation in widespread island taxa in comparison with mainland relatives (Francisco-Ortega et al. 2000;Su et al. 2010).In species capable of long-distance dispersal genetic differentiation among islands may be reduced, and intra-population genetic diversity might be higher than otherwise expected (García Verdugo et al. 2014;García-Verdugo et al. 2015).
With a view to contributing to more global study on J. phoenicea, for the current study, we report on J. turbinata using AFLP markers.Our main aim was to determine the genetic structure of the Canarian population revealed by neutral nuclear markers, and to evaluate whether its structure is more generally affected by geographical isolation or by gene flow.As such, subpopulations of J. turbinata were analyzed: a) to assess genetic diversity across the Canarian populations, and b) to infer the influence of isolation in the genetic structure of J. turbinata.

Plant Sampling
As J. turbinata occurs mainly as dispersed individuals in the Canary Islands, the strategy was to sample at locations in which individuals were grouped forming populations.So, a population from each of the five Canarian islands where J. turbinata occurs was sampled.Three populations from mainland Morocco were also sampled to compare values of genetic diversity with those obtained for the island populations.So, 175 individuals belonging to eight populations were sampled.The locations of sampled populations are shown in Fig. 1 and their characteristics are listed in Table 1.To avoid damaging individuals, the sampling strategy consisted of collecting apexes of branches measuring about 10 cm in length.Furthermore, we only sampled mature individuals that were sufficiently spaced.All plant material was stored in labeled ziplock bags with silica gel at room temperature to dry it until DNA extraction (Chase and Hills 1991).

DNA analyses
DNA was extracted from silica gel dried leaves using the 2 × cetyltrimethylammonium bromide (CTAB) method (Doyle and Doyle 1987).Total DNA extracts were quantified using a Nanodrop 2000 spectrophotometer.Amplified fragment length polymorphism (AFLP) analysis was carried out following Vos et al. (1995) with minor modifications as follows.The analysis was performed with fluorescence-labeled primers (VIC, NED, PET, Applied Biosystems, Madrid, Spain) instead of radioactive labeled primers.Fragments were selectively amplified with the primer pairs EcoR1-ATG/MseI-CGT, EcoRI-ACA/MseI-CGT, EcoRI-ACG/MseI-CAC.Multiplex products were run for 4 h on an ABI 377 sequencer to separate fragments together with an internal size standard (GeneScan 600 LIZ, ABI).Sizing and peak identification were performed using Genemapper 4.0 software (Applied Biosystems).In the finally assembled binary matrix the presence of a band was scored as 1 while the absence of a band was scored as 0.
As AFLP markers are dominant, we assumed that null bands were homologous and that populations were in Hardy-Weinberg equilibrium (Lynch and Milligan 1994) in order to compute diversity indices [percentage of polymorphic markers (PLP) and Nei's (1978) unbiased expected heterozygosity (He)] and genetic distance among populations (Fst).These parameters were inferred with AFLP-surv 1.0 (Vekemans 2002).Significance of Fst values was determined using 1000 bootstrapped data sets.
A Bayesian model-based analysis was performed to infer population structure with Structure version 2.2 (Pritchard et al. 2000;Falush et al. 2007).The F model, based on an admixture ancestry model with correlated allele frequencies, was imposed to estimate the posterior probabilities [LnP(D)] of K groups (Pritchard and Wen 2004) and the individual percentages of membership assigned to them according to their molecular multilocus profiles (Falush et al. 2003(Falush et al. , 2007)).Probabilities for a range of K were examined starting from one to the number of sampled populations plus one (K = 1-9), using a burn-in period and run length of the Markov chain Monte Carlo We used CLUMPP 1.1.2(Jakobsson and Rosenberg 2007) to reach a consensus on the results of the independent runs for the optimal K.For the consensus, we used the Greedy option with random input order and 10 000 repeats.The consensus was visualized in DISTRUCT 1.1 (Rosenberg 2004).Alternatively, another bayesian analysis of the genetic structure was conducted with BAPS 5.4 (Corander and Marttinen 2006), using stochastic optimization instead of Markov chain Monte Carlo to find the optimal partition.We performed a mixture analysis of individuals without using the geographic origin of the samples as a prior informative (clustering of individuals).BAPS simulations were run with the maximal number of groups (K) set from 1 to 10.Each run was replicated 10 times, and results were averaged according to the resultant likelihood scores.The output of the mixture analysis was used as the input for population admixture analysis (Corander and Marttinen 2006), with the default settings to detect admixture between clusters.
In addition, two Principal Coordinate Analyses (PCoA) were performed to illustrate overall similarity among individuals using Multivariate Statistical Package software (MVSP version 3.12d http://www.kovcomp.com/mvsp).The first one included all the sampled individuals, while the second one was performed using only Canarian populations, in order to detect any substructure of genetic diversity.PCoA was inferred from the pairwise Nei's genetic distance (Nei 1978) between all pairs of AFLP phenotypes.
Pairwise Fst-values between all populations were also calculated and tested for significance by resampling with 1000 random permutations.Mantel tests (Mantel 1967) were carried out to assess correlation among genetic and geographic distances.These analyses were performed using all populations, and alternatively using only Canarian populations.Thereafter, the fine-scale genetic structure among the Canarian populations was investigated using spatial autocorrelation analysis in GenAlEx 6.5 (Peakall and Smouse 2006).The spatial autocorrelation coefficient (r) was computed using the geographic distance and the genetic distance between individuals.Tests for statistical significance were conducted by random shuffling (999 times) of individual geographic locations to define the upper and lower bounds of the 95% confidence interval for each distance class, and estimating 95% confidence intervals around mean r values by bootstrapping pair-wise comparisons within each distance class (1000 repeats).Analyses were performed using the variable distance classes option based on the geographical distances estimated in each population.

Results
The three selective AFLP primer combinations amplified a total of 589 fragments.In total 94.1% of the fragments were polymorphic and every of 175 individuals had unique AFLP profiles.No private markers for island were detected, but we found five exclusive bands for Canary Islands, which appeared in more than 80% of sampled individuals, and two exclusive bands for the Moroccan populations.The proportion of polymorphic loci (PLP) for a single population ranged from 65.0% (MOR2) to 73.7% (PALM).Expected heterozygosity values showed that PALM was the most diverse population (He = 0.218), whereas HIER showed the lowest within-population diversity (0.192).The average gene diversity within populations (Hs) was 0.209, whereas the total diversity (Ht) value was 0.243 (Table 1).The fixation index was highly significant (Fst = 0.1398, P < 0.000), suggesting low to moderate genetic differentiation among populations.
The Bayesian analysis of genetic structure of J. turbinata populations conducted with Structure found the highest LnP(D) and ΔK values for K = 3. Moroccan populations were assigned to one cluster and Canarian populations were assigned to other two clusters; the first of them grouped individuals of GOME, TENE, and PALM with high proportion of membership and the last cluster grouped HIER, whereas GCAN individuals showed proportion of membership intermediate between these two clusters.The genetic structure revealed by BAPS was highly congruent with STRUCTURE results identifying three clusters (P = 0.999).However, only 5.14% of individuals were admixed (Fig. 1).The first cluster grouped the individuals from Morocco.The second one grouped GOME, TENE, PALM and seven individuals from GCAN, whereas the third cluster included HIER and 18 individuals from GCAN.The individuals admixed were four from GCAN, three from PALM and two from TENE.PCoA analyses gave similar results (Fig. 2).PCoA plot performed including Moroccan population revealed two clusters, one comprising the samples of Morocco and a second cluster that grouped Canarian populations (Fig. 2a).For this analysis the three first axes accounted for 22.76% of the variation (15.51%, 4.33% and 2.92%, respectively).PCoA analysis for the Canarian populations revealed three clusters (three first axes explained 8.67%, 4.05% and 2.80% of the variation respectively).The first cluster was constituted by individuals from TENE, GOME and GCAN.The second grouped individuals form populations of TENE, GOME, GCAN and PALM.The third was formed by HIER population, including some individuals from GCAN (Fig. 2b).
Pairwise Fst among populations showed high differentiation between Moroccan and Canarian populations.In addition HIER population was moderately differentiated from TENE and GOME populations (Table 2).The remaining populations revealed low pairwise genetic differentiation.The Mantel test in which were included Moroccan populations showed a great isolation by distance pattern (R = 0.909, P = 0.001).When these populations were eliminated from the analysis the correlation among genetic and geographic distances was not significant (R = 0.301, P = 0.085).Spatial autocorrelation analysis showed positive isolation by distance at distances of 90 km approximately (Fig. 3).

Genetic diversity
It is well documented that the genetic diversity of small island species is usually lower than those of their mainland congeners (Frankham 1997;Whittaker et al. 2008), which is a result of founder effects or bottlenecks.In our study, genetic diversity estimates obtained with AFLP do not follow predicted patterns.Although Canarian populations are differentiated from Moroccan populations, PLP and He suggest that genetic diversity is not impoverished.This evidence is not unusual, as García-Verdugo et al. (2015) suggest that in more than 50% of studies in which genetic diversity of island plants is addressed, island populations show similar or higher genetic diversity than their mainland relatives, with possible causes being the ability of some plant species to recurrently overcome sea barriers, or the role of islands as climatic refugia (García-Vergudo et al. 2015 and references therein).With respect to other studies on the genetic diversity of J. turbinata populations is acknowledged that genetic diversity values obtained with molecular markers are difficult to compare between studies, because the statistical methods used to process the data often vary, particularly with dominant and co-dominant markers.However, parameters such as intra-population genetic markers may serve to establish comparisons between species (Nybom 2004).Regarding the He values found in our study, studies on J. phoenicea group by Boratyński et al. (2009) and Dzialuk et al. (2011) produced similar results to ours, while Meloni et al. (2006) obtained higher Hs values.The intra-population genetic diversity observed in our study falls within the range of the values reported by Nybom (2004) for outcrossing, long-lived endemic species.
Our results therefore imply that, despite the low number of individuals sampled, as well as the documented historical decrease of populations throughout the Canarian archipelago (Fernández-Palacios et al. 2008), J. turbinata has moderate levels of genetic diversity (evidenced by AFLP) with no apparent evidence pointing to genetic bottlenecks.Long-lived woody perennials are expected to be resilient to changes in genetic diversity due to ample gene flow and long generation times.Furthermore, fleshy fruited species typically show widespread distributions coupled with extensive gene flow among islands (García-Verdugo et al. 2014).Therefore, Canarian populations of J. turbinata may have been able to avoid the effects of genetic drift through long-distance dispersal, which has previously been described for several Macaronesian species (Ferreira et al. 2001;Vargas 2007;García-Verdugo et al. 2010;García-Verdugo et al. 2014;Rumeu et al. 2014).

The genetic structure of populations
The Fst pairwise estimates obtained in this study showed high levels of genetic differentiation between MOR and Canarian populations and low to moderate population differentiation among Canarian populations (paiwise Fst < 0.13).Furthermore, the Mantel test and spatial autocorrelation analysis failed to predict an isolation-by-distance pattern when only Canarian populations were analyzed.This pattern of differentiation suggests an ancient colonization of the Canary Islands by mainland populations, with a posterior interruption of gene flow with African populations.This is plausible if we taking into account that the estimated age of J. phoenicea s.l. is more than 30 Mya (Mao et al. 2010).According to the order of emergence of islands, it is reasonable to think that colonization occurred from eastern to western islands.However, with our available data is difficult predict if it was a single or multiple colonization events, or the role of each island in these hypothetical events.The levels of genetic differentiation between mainland and island populations suggest an anagenetic type of evolution for the Canarian populations, in which genetic variation would have been accumulating for generations through mutation and recombination until the overall level of genetic variation may approximate to that of source population (Stuessy et al. 2006).
However, the distribution of genetic variation suggests substantial gene flow between almost all the islands, which seems to be less effective between El Hierro and the others.Surprisingly, the PCoA, the admixed individuals depicted by BAPS and the pairwise Fst suggest that the Gran Canaria population is most similar to the El Hierro population, which is the most geographically remote of the islands.This finding would suggest more complex patterns of gene flow among islands have taken place than the predicted stepping-stone model of differentiation from east to west, such as that observed for Gallotia spp.(González et al. 1996) or Olea europaea L. (Hess et al. 2000).Historical geologic events might have played an important role in the partitioning of genetic diversity by promoting local extinctions and re-colonizations, as it has been observed for several groups of animals and plants (Brown and Pestano 1998;Nogales et al. 1998;Juan et al. 2000;Fernández-Mazuecos and Vargas 2011).Unfortunately, due to the scarcity of well defined populations in Canary Islands, we sampled a population per island, instead individuals from different locations in each island.Therefore, observed genetic differences could be due to this sampling strategy.Adding georeferenced individuals from different locations in further studies might elucidate patterns of gene flow among islands, colonization events, etc.
Furthermore, the documented arrival of settlers to Canary Islands might have been a cause of intra-population genetic diversity decrease and differentiation increase, through the transformation of favourable habitats in grazing (Fernández Palacios et al. 2011).However, as demonstrated for several Canarian endemics (García-Verdugo et al. 2010;Ferreira et al. 2011;Rumeu et al. 2011;Sosa et al. 2013), the recurrent gene flow between islands have revealed moderate levels of genetic diversity and absence of genetic structure across islands.The ability of J. phoenicea s.l. to disperse both pollen (by wind) and seeds (by birds) to distant locations has already been documented in mainland populations (Jordano 1993).The known dispersers of J. turbinata in Canary Islands are ravens, blackbirds and lizards (Nogales et al. 1999;Fernández-Palacios et al. 2008;Otto et al. 2010).Long distance dispersal (LDD) mechanisms in oceanic islands involves several types of dispersal; anemochory, zoochory, by birds capable to fly over long distances, as well as diplochory or secondary seed dispersal, in which the primary disperser, is hunted by a predator capable to travel between islands, and sea dispersal or hydrochory.(Nogales et al. 2012;Vargas et al. 2014).In this case dispersal by birds and diplochory might be effective methods of LDD (Nogales et al. 1999(Nogales et al. , 2002(Nogales et al. , 2007;;Padilla et al. 2012).Furthermore, although J. turbinata do not exhibit specialized floating structures, dry galbulus with fertile seeds are capable to float, so hydrochory could also be considered as an effective LDD mechanism.These abilities could have led to an important interexchange of pollen and seeds between islands that may mask any ancient relationships.
As yet, while we do not have any data on the distribution of cpDNA haplotypes for our study populations, such information may enable us to infer the underlying processes driving colonization more accurately, as previously demonstrated for several other Canarian endemics such as J. cedrus Webb.& Berth.(Rumeu et al. 2014), or O. europaea subsp.guanchica P. Vargas, J. Hess, Muñoz Garm. & Kadereit (García Verdugo et al. 2009).Further studies are therefore recommended to shed more light on the specific processes underlying island colonization.

Adknowledgements
Authors wish to thank to David López García and Juan Bautista Vera for field assistance in collecting samples.This work has been supported by a Ministry of Economy and Competitiveness CGL2011-30099 project grant to Pedro Sánchez Gómez.

Fig. 1 .
Fig. 1.(a) Location map and Bayesian analysis of the genetic structure of 175 individuals from eight sampled populations of Juniperus turbinata.The great rectangle that includes Canary Islands represent an enlargement of the same area in the small rectangle Pie charts represent the mean proportion of membership of each population to the predefined K = 3 clusters with BAPS.(b) Plots show the membership of individuals to the predefined K = 3 clusters obtanined with b1, STRUCTURE and b2, BAPS.Number and codes of populations are listed inTable 1.

Fig. 2 .
Fig. 2. Principal coordinate analysis of the 175 individuals from eight sampled populations of J. turbinata based on pair-wise Nei's (1978) genetic distances.(a) Scatter plot showing eight populations based on two first components of the PCoA.(b) Scatter plot showing 115 individuals from five Canarian populations.Number and codes of populations are listed inTable 1.

Fig. 3 .
Fig. 3. Correlogram from spatial autocorrelation analysis using the correlation coefficient r, and variable distance classes.95% confidence error bars for r were estimated by bootstrapping over pairs of samples; dashed lines (U, L) represent upper and lower bounds of a 95% CI for r generated under the null hypothesis of random geographic distribution.

Table 2 .
Pairwise estimate Fst values among eight populations of J. turbinata.All values are significant at 0.001 level.Number and codes of populations are listed in Table 1.