Podocarpus in the palaeogeographically complex island of Hispaniola: A stepping-stone colonization and conservation recommendations

Aim: Hispaniola is the second largest island in the Caribbean and a hot spot of bio-diversity. The island was formed by the fusion of a northern and southern palaeo-islands during the mid-Miocene (15 Ma). The historical split of Hispaniola together with repeated marine incursions during the


| INTRODUC TI ON
The Caribbean archipelago is the largest, most species-rich insular system in the Neotropics (Maunder et al., 2008), and one of the world's 34 biodiversity hot pots sensu Mittermeier et al. (2004).In situ speciation is one of the leading hypotheses proposed to explain its high biodiversity (Roncal et al., 2020) which is supported by the high number of endemics (up to 8000 seed plant species; Acevedo-Rodríguez & Strong, 2008).Endemic clades thus provide opportunities for understanding the evolutionary history of insular biotas.Previous studies have aimed at reconstructing the colonization history and adaptive radiations of Caribbean lineages with a bias for vertebrates (e.g.Losos et al., 2006;Nieto-Blázquez et al., 2021;Reynolds et al., 2016;Turvey et al., 2016).The heterogeneous environment of large and geologically complex islands has been associated with historical and modern barriers to gene flow that can drive speciation (Turvey et al., 2016).The evolutionary history of intraisland endemic plant lineages, that may or may not surpass these barriers, will thus increase our understanding of how effective these barriers are for intra-island diversification.Moreover, the distribution of genetic diversity in large islands will provide insight on gene flow in space and will inform local conservation strategies for insular taxa (Brown et al., 2014).
The island of Hispaniola is the result of the juncture of two palaeo-islands that collided tectonically during the middle Miocene, ca. 15 Ma (Graham, 2003;Mann et al., 1991).These northern and southern palaeo-islands remained separated by the Neiba Valley-Cul de Sac Plain (Figure 1), an arid and deep rift that was inundated by marine incursions repeatedly over the late Pleistocene (Maurrasse et al., 1980;McLaughlin et al., 1991).The southern palaeo-island is further divided into two physiographic provinces (Massif de la Hotte to the west and Massif de la Selle-Sierra de Bahoruco to the east) separated by the Jacmel-Fauche' depression (known as Bond's line), which was also inundated during the Plio-Pleistocene through a sea channel (Maurrasse et al., 1980).These geological events left three distinct biogeographic regions disconnected during the incursion periods generating allopatric populations and genetic divergence in a few endemic birds and mammals (Brace et al., 2012;Sly et al., 2010Sly et al., , 2011;;Townsend et al., 2007;Turvey et al., 2016), but the genetic or evolutionary effects on other lineages are less understood including vascular plants and thus preventing the generalization of this model.Hispaniola's current topography consists of mountain chains (cordilleras) that run in parallel from NW to SE (Latta et al., 2006) separated by lowland xeric valleys (Heubeck & Mann, 1991;Townsend et al., 2007).The largest mountain system is the Cordillera Central, which was uplifted during the middle to late Eocene and is dominated dispersal barriers, Hispaniola, historical demography, in situ conservation, IUCN endangered, montane forest, Podocarpus, population structure, SNPs, stepping-stone colonization F I G U R E 1 Topographical map of Hispaniola indicating the main mountain ranges and sampling collection sites for Podocarpus on Hispaniola.CC, Cordillera Central; CS, Cordillera Septentrional; SB, Sierra Bahoruco; SN, Sierra Neiba.Black dotted lines indicate the location of the Neiba Valley-Cul de Sac Plain that was subject to marine incursions during the Pleistocene separating the south and north palaeo-islands, and the Jacmel-Fauche depression (Bond's line), which divides the southern palaeo-island into east and west.Inset shows location of Hispaniola within the Caribbean and tropical America.Red dots indicate collection sites of P.buchii and blue dots collections of P.hispaniolensis.Squares indicate the five populations used in the historical demographic analysis.Black squares and circles correspond to P.buchii and P.hispaniolensis specimen localities, respectively, representing the entire distribution range for each species (Mill, 2015) by igneous and volcanic materials.The Sierra Neiba, a karstic range in the south of Hispaniola (or southern range of the northern palaeoisland), and the Sierra Bahoruco, a limestone mountain in the southcentral part of Hispaniola, are of similar age to the Cordillera Central.
The Cordillera Septentrional, located in the north of Hispaniola (or northern range of the northern palaeo-island), is of later, Oligocene-Miocene origin and presents mostly sedimentary rocks (Figure 1; Cano-Carmona et al., 2010;Cano-Ortiz et al., 2016).These mountain chains have been invoked as the explanation for the high levels of genetic differentiation among populations of a dry lowland lizard species (Gifford et al., 2004).Other current ecological factors that have caused lineage divergence of Hispaniolan populations are geographic macrohabitats (in lizards, Glor et al., 2003) and unoccupied feeding niches (e.g. in butterflies, Matos-Maraví et al., 2014).
The conifer genus Podocarpus is represented by two endemic sister species on Hispaniola, P.buchii Urb. and P.hispaniolensis de Laub., that diverged 15 mya (Nieto-Blázquez et al., 2021), and coinciding with the junction of the two palaeo-islands.Podocarpus colonized the Greater Antilles from the Andes 45 to 31 mya, and a vicariant origin for the ancestor of Hispaniolan Podocarpus attributed to the tectonic separation of eastern Cuba and western Hispaniola during the middle Miocene has been proposed (Nieto-Blázquez et al., 2021).The sister species differ in morphology, elevational range and distribution.
Podocarpus buchii can be distinguished from P.hispaniolensis by the absence of a prominent deep groove on the adaxial leaf surface, the smaller terminal buds usually less than 3.5 mm (not >5 mm), the short female peduncles less than 4 mm long (not >6 mm) and the smaller male cones less than 10 mm (not >12 mm, Mill, 2015).Podocarpus buchii tends to be found at higher elevations (1100-2500 m) than P.hispaniolensis (800-1200 m).Podocarpus hispaniolensis is found in the northernmost mountain chain (Cordillera Septentrional) of the Dominican Republic and historically in the Massif du Nord in Haiti (Mill, 2015).Podocarpus buchii is found in the south, particularly in the Sierra de Bahoruco and Neiba of the Dominican Republic, and also in the Massif de la Hotte and Chaine de la Selle in Haiti (Mill, 2015, Figure 1).Both species occur in the Cordillera Central, and although the two species do not typically occur in sympatry, collections that include both species have been made in relatively small areas.Seeds are probably dispersed by frugivorous birds and small mammals attracted by the colour and/or swollen cone bract (Enright & Jaffré, 2011;Mill, 2003).Several bird families feed on the fleshy female cones of the congeneric P.parlatorei Pilg.transporting its seeds (Blendinger, 2017).
Both Hispaniolan species are IUCN endangered because of habitat loss due to logging, forest clearance for agricultural and commercial plantations, fire and mining, and selective logging of its highly valuable timber (Gardner, 2013).The most severe habitat degradation and loss has occurred in Haiti.The trend of population decline for P.buchii has been estimated at 50% over three generations, resulting in currently severely fragmented populations (Gardner, 2013).
Podocarpus hispaniolensis has an estimated area of occupancy of less than 5000 km 2 with about five locations severely fragmented.Both species are recorded in a few protected areas where threats are still present.There is no action recovery plan, and no local species management including ex situ conservation for either species (Gardner, 2013).
Vulnerable and endangered species usually consist of small and relatively isolated populations that maybe impacted by the detrimental effects of reduced genetic diversity due to drift, limited gene flow and inbreeding that can threaten their long-term persistence (Frankham et al., 2002).Although the impact of genetic factors on species' well-being has been controversial (Spielman et al., 2004;Teixeira & Huber, 2021), populations occurring in small numbers should be the focus of conservation actions because they may be prone to suffer from environmental and demographic stochasticity (Lande, 1993), the latter usually driven by anthropogenic factors.
Nonetheless, with the advent of the genomic era, we have a unique opportunity to analyse in more detail populations and species under threat in order to contribute to their better conservation and management (Supple & Shapiro, 2018).
Our study aims to advance our knowledge of whether multiple Hispaniolan lineages exhibit similar evolutionary history and population structure in response to historical allopatric barriers (e.g.palaeo-islands, marine incursions), or whether the significance of these historical barriers in regulating gene flow is now masked by current landscape factors (e.g.topography, lowland xeric valleys).
Our research (1) estimates genetic diversity of populations, (2) identifies genetic structure within and between species, (3) infers the phylogenetic relationship among populations, (4) elucidates the demographic history of the genus on the island and (5) recommends strategies of long-term genomic conservation of Hispaniolan Podocarpus.We tested four alternative hypotheses for analysing the historical demography: (1) an ancestral population from Cordillera Central gave rise to both Podocarpus species from where one species migrated south and the other north, (2) allopatric speciation in the two palaeo-islands followed by secondary contact in Cordillera Central, (3) a stepping-stone colonization from south to north and (4) a stepping-stone colonization from north to south.

| Sample collection and DNA extraction
We collected silica-dried leaf samples from 20 P.buchii and 16 P.hispaniolensis, across 11 collection sites in the Dominican Republic across all four mountain chains (Table 1, Figure 1).The low sample size reflects their rarity in the field, as recognized by their status as IUCN Endangered and the alarming estimate that forest cover in the Dominican Republic is down to 14% (Russell, 1991).Very few individuals were thus available for sampling.Duplicate herbarium specimens from each individual were deposited at the National Herbarium of the Jardín Botánico Nacional Dr. Rafael Ma.Moscoso (JBSD) and the Royal Botanic Garden of Edinburgh (E, Table S1).We selected adult individuals with healthy-looking leaves for both silica-dried leaf samples and herbarium specimens.We selected as out-group single specimens from 11 of the 13 Podocarpus species in the tropical American clade to which Hispaniola's Podocarpus belongs (P.angustifolius, P.aristulatus, P.coriaceus, P.ekmanii, P.guatemalensis, P.matudae, P.oleifolius, P.purdieanus, P.sellowii, P.trinitensis, P.urbanii; Quiroga et al., 2016).We obtained these out-group samples from the Royal Botanic Garden Edinburgh and the Montgomery Botanical Center (Table S1).
We isolated DNA from 35 to 40 mg of dry leaf tissue with the DNeasy Plant MiniKit (Qiagen).We modified the manufacturer's standard protocol to improve DNA recovery.We increased AP1 buffer (for cell lysis) from 600 to 750 μl, the time of incubation for cell lysis up to 60 min and P3 buffer (for precipitation of polysaccharides, detergent, and proteins) from 195 to 225 μl.We diluted DNA extractions in EB buffer to 20 ng/μl.

| Genotyping-by-sequencing, SNP discovery and genetic diversity
The Institut de Biologie Intégrative et de Systèmes (IBIS) of the University of Laval in Canada conducted the genotyping-bysequencing (GBS).Genomic libraries were prepared for the 47 DNA samples with two restriction enzymes, SbfI (high fidelity) and MspI (New England BioLabs Inc.).Unique barcodes 10-12 bp long were added to each sample to facilitate subsequent demultiplexing.Library preparation and sequencing followed the protocol of Abed et al. (2019).We obtained single-end sequencing reads of variable length (up to 200 bp) from two chips of an Ion Proton system.We inspected data quality with FastQC 0.11.5 (Babraham Bioinformatics), which reports Phred quality scores and % GC content.
We conducted SNP discovery with Stacks 1.47 (Catchen et al., 2013).We demultiplexed and filtered sequencing reads using the process_radtags function.We trimmed reads to 64 bp length and eliminated uncalled reads and reads with low-quality scores (phred score of 10).Since we did not have a reference genome, we conducted a "de novo" assembly of loci with the module "ustacks." In this step, reads were aligned into matching blocks or stacks per sample.Then, the stacks were compared, a set of putative loci produced and SNPs detected at each locus.We followed the recommendations of Paris et al. (2017) on the selection of parameter values through the pipeline, keeping those that rendered the highest number of loci.We used -m (minimum depth coverage to create a stack) = 5; -M (maximum distance allowed to create a stack) = 3; -p (parallel execution of several threads) = 15; and the -gapped option which allowed gapped alignments between stacks.Subsequently, the module "cstacks" built a catalog from the loci produced in "ustacks" by merging stacks with at most 3 (-n = 3) mismatches between loci and allowed parallel execution of several threads (-p = 15).The module "sstacks" matches the loci produced by "ustacks" with the catalog produced by "cstacks," and Note: Collection site numbers as in Figure 1, with acronyms in parenthesis indicating the designated population for historical demographic analysis in DIYABC.# indiv = number of individuals collected per site.Individuals from collection sites marked with an * occur in protected areas (Figure S5).

TA B L E 1 Sample collection information for Podocarpus on Hispaniola
Abbreviation: DR, Dominican Republic.
again parallel execution of several threads (-p = 15) was allowed.
We used the "populations" module of Stacks to obtain populationlevel summary statistics, genetic diversity indexes (i.e.observed and expected heterozygosity, nucleotide polymorphism diversity), F ST values and a SNP dataset.We chose the following parameters: -p (minimum number of populations containing each locus) = 1; -r (minimum percentage of individuals present in a population to process a locus for that population) = 0.4; and -min_maf (minimum allele frequency) = 0.1.We ran the "populations" module of Stacks twice.In the first run, we assigned samples to the 11 collection sites plus the out-group, and in the second, we assigned samples to five populations (see Section 3).The first SNP dataset was used to run a phylogenetic reconstruction, while the second to run the genetic structure, and demographic history analyses described below.
We performed a Mantel test to detect isolation by distance (IBD).
We used the R package "ade4" (Chessel et al., 2004) to measure the correlation between a matrix of pairwise F ST values among the 11 collection sites and their geographic distances.We calculated pairwise geographic distances in metres with the Geographic Distance Matrix Generator (Ersts, 2011).The mantel test uses the Pearson coefficient and a default alpha value of 0.05.

| Genetic diversity and structure of Podocarpus on Hispaniola
To identify loci under selection, we performed an outlier analysis using the Bayesian-likelihood method implemented in BayeScan 2.1 (Foll & Gaggiotti, 2008) using the following parameters: number of outputted iterations was 5 × 10 3 , thinning interval size of 10, 20 pilot runs of length 5 × 10 3 each and burn-in length of 50 × 10 3 .We applied a false discovery rate of 0.05.Evanno et al., 2005).Finally, we used distruct OSX10 (Rosenberg, 2004) to graphically display the STRUCTURE results.

| Phylogenetic relationships
To determine the phylogenetic relationships among collection sites, we ran maximum-likelihood (ML) and Bayesian inference (BI) analyses using the CIPRES Science Getaway (Miller et al., 2010), including all individuals per collection site.We ran the ML analysis on RAxML-HPC2 v8.2.10 (Stamatakis, 2006) where we used the nucleotide substitution model GTR+Γ and a rapid bootstrap (BS) algorithm with 500 replicates.We ran the BI analysis with the GTR+G substitution model on MrBayes v3.2.6 (Ronquist et al., 2012) for 50 million MCMC generations sampling every 10,000 generations, and a burnin of 25% generations discarded.We selected 11 tropical American Podocarpus species as the out-group.To further elucidate the pattern of Podocarpus population divergence, we also ran ML and BI phylogenetic analyses for the five populations used in the historical demographic analysis using the same parameters as described above.

| Demographic history of Podocarpus on Hispaniola
To infer the demographic history of Podocarpus on Hispaniola, we conducted an approximate Bayesian computation analysis in DIYABC v 2.1.0(Cornuet et al., 2014).We defined five populations based on

| Genotyping-by-sequencing and SNP discovery
The  1).Pairwise F ST values among the 11 collection sites ranged from <0.001 to 0.160 (Table S3).The Mantel test showed a significant correlation between the F ST matrix for the 11 collection sites and their geographic location (R = 0.442, p-value = 0.003), which suggests the presence of isolation by geographic distance.
The STRUCTURE analyses with and without sample assignment into populations gave very similar results, and thus, we present results with population assignments.The full dataset with both species combined contained 22,657 SNPs, of which 22,609 were neutral, and 48 were under selection.The genetic structure analysis for P.buchii used 15,524 SNPs and for P.hispaniolensis used 5687 SNPs.For P.buchii, the highest supported K was 4 followed by K of 2. Both K = 2 and K = 4 showed individuals from Sierra Bahoruco as a distinct genetic cluster.Relative to K = 2, K = 4 distinguished individuals in Cordillera Central from those in Sierra de Neiba (Figure S1).Two individuals from Sierra Bahoruco and the four from Cordillera Central showed admixture for K = 2, while individuals from Sierra Neiba did not show any admixture.
For K = 4, four individuals from Sierra Bahoruco, two from Sierra Neiba and one from Cordillera Central showed admixture (Figure S1).For P.hispaniolensis, the highest supported number of clusters was K = 2 followed by K = 5.Septentrional showed admixture (Figure S1).
When both species and all loci were analysed, K = 7 received the highest support followed by K = 5.For both K values, P.buchii individuals from Sierra Bahoruco, Sierra Neiba and Cordillera Central belonged to different genetic clusters, but P.hispaniolensis individuals from Cordillera Central and Cordillera Septentrional belonged to one genetic cluster (Figure 3a).Podocarpus hispaniolensis shared genetic similarity with all P.buchii genetic clusters.
The individuals that showed the most admixture corresponded to P.buchii from Sierra Bahoruco, while all but one P.hispaniolensis individuals did not show any admixture (Figure 3a).When only neutral loci were analysed, K = 5 received the highest support, while a K = 4 was preferred when loci under selection were analysed (Figure 3b).STRUCTURE analyses with all and neutral loci separately showed a similar geographic genetic structure by mountain chain and species identity.However, this geographic genetic structure was not clearly observed with the loci under selection (Figure 3c).Genetic structure results for the two Podocarpus species separately and combined supported the assignment of five populations for the historical demographic analysis as defined in the Section 2 and Table 1.

| Phylogenetic relationships
ML and BI analyses for the 11 collection sites recovered a congruent topology and showed that P.hispaniolensis was nested within P.buchii  Both analyses also recovered Site 1 (P.buchii from Sierra Bahoruco Oriental) as sister to the P.hispaniolensis clade with high support (PP = 1.0 and BS = 99%).The relationships of P.buchii sites closer to the root showed some discrepancies between the two phylogenetic reconstruction analyses, but recovered both P.buchii Cordillera Central sites (7 and 8) as sister to each other with high support

| Demographic history of Podocarpus on Hispaniola
We discarded loci that had missing data in all individuals from a given population, as required by DIYABC, resulting in 244 SNPs for this analysis.The approximate Bayesian computation identified Models 6 and 8, the stepping-stone hypothesis with bottlenecks for each mountain colonization, as the best-fit under the direct and logistic regression approaches (Table S4 and Figure S2).Models that hypothesized an ancestral population from Cordillera Central that gave rise to two Podocarpus species or an allopatric origin for the two species in the two palaeo-islands reaching secondary contact in the Cordillera Central received no support (Figure 2, Table S4).The posterior predictive global error was 0.008 under the direct approach and 0.004 under the logistic regression approach.The prior predictive global error was 0.244 under the direct and logistic regression approaches.Direct approach type I errors were 0.06 and 0.2 for Models 6 and 8, respectively.Logistic type I errors were 0.036 and 0.1 for Models 6 and 8, respectively.Datasets from the posterior predictive distribution were close enough to the observed dataset for both Models 6 and 8 (Figures S3 and S4).Some summary statistics for both models had a proportion of simulated datasets below the observed one (Tables S5 and S6).The estimation of mean population sample sizes ranged from 386 to 1270 for Model 6 and 373 to 1840 for Model 8. Mean coalescent times were 15,500 years at t4, 1837 years at t3, 1010 years at t2 and 843 years at t1 for Model 6 (Table S7) and 3674 years at t4, 3325 years at t3, 3475 years at t2, and 1925 years at t1 for Model 8 (Table S8).Parameter estimates for Models 6 and 8 were accurate and unbiased, with low mean relative bias and RRMISE values (Tables S9 and S10, respectively).

| DISCUSS ION
In the paleogeographically complex Caribbean island of Hispaniola,

| Patterns of genetic diversity and structure of Podocarpus on Hispaniola
The higher genetic diversity in P.buchii is likely due to its earlier diverging position with respect to the derived species P.hispaniolensis as shown in the phylogenetic tree of collection sites.A similar geographic distribution of genetic diversity was reported for Hispaniolan hutias, where southern populations had a higher nucleotide diversity (π) compared with northern populations (Brace et al., 2012).Yet, other taxa, such as the palm Pseudophoenix vinifera (Mart.)Becc., did not show significant differences in genetic diversity (H e ) between northern and southern populations (Rodríguez-Peña et al., 2014).
Geographic distance is a major cause of genetic isolation (isolation by distance), which seems to occur in Podocarpus populations from Hispaniola.Even in small islands, such as Lord Howe Island in Australia, isolation by distance plays a role in lineage divergence in plants (Papadopulos et al., 2014).However, genetic distances are not always linearly correlated with geographic distances and other factors such as habitat configuration and maximum migration distance might influence isolation among populations (van Strien et al., 2015).
This might explain why our most geographically distant collection sites (Sites 1 and 2 vs. 10 and 11) did not show the greatest genetic differentiation in our study.Overall genetic differentiation among collection sites (pairwise F ST values) was low.
Contrasting spatial patterns of genetic structure have been observed on Hispaniolan fauna.For example, in mammals of the genus Plagiodontia (Brace et al., 2012), and Solenodon (Turvey et al., 2016), and in birds (Sly et al., 2011), the distribution of genetically distinct allopatric lineages coincides with the geological fragmentation of the south and north palaeo-islands.Similarly, Matos-Maraví et al. (2014) showed that two vicariant events prompted the divergence between Calisto butterflies, one related to the uplift of the Cordillera Central, and a second one related to the repeated marine incursions in Neiba Valley-Cul de Sac Plain.In contrast, it is likely that for Podocarpus on Hispaniola, the dry valleys that lay in between the mountain chains acted as dispersal barriers that prompted divergence.This geographical structure is concordant with that of the lizard Ameiva chrysolaema (Gifford et al., 2004), where the complex topography of

| Demographic history of Podocarpus on Hispaniola
The approximate Bayesian computation analysis supports a steppingstone colonization model with bottleneck events with each new mountain range migration.However, the direction of colonization (south to north or vice versa) could not be elucidated unequivocally.
A north to south colonization direction within Hispaniola was shown for chat-tanagers (Calyptophilus frugivorus) from the north palaeoisland into Sierra Bahoruco in the south palaeo-island (Townsend et al., 2007).
The estimated population coalescence times for Podocarpus are young compared with those of the Plagiodontia mammals, where the south and north populations diverged 0.594 Ma, during the time when marine incursions in the Neiba Valley-Cul de Sac Plain were still occurring (Brace et al., 2012).Given the inferred coalescence times, Podocarpus divergence history on Hispaniola occurred after the island achieved its current configuration: the north and south palaeoislands were already connected, the main mountain ranges were in place, and all marine incursions no longer occurred.Therefore, the marine incursions and geological reconfigurations did not leave a signature on the population divergence history of Podocarpus.We acknowledge that the small sample size in the present study could have rendered the rather young coalescence times, and thus, they should be corroborated with a larger Podocarpus sampling including seedlings when adults are unavailable.
The historical demography, phylogenetic analysis and species range distribution suggest that an ecological speciation event occurred in Cordillera Central probably triggered by different selection pressures under different elevations.This hypothesis, however, requires further testing by identifying loci under divergent selection in an elevation gradient.Glor et al. (2003) showed that lineage divergence and genetic structure of Anolis lizards on Hispaniola were more likely the result of macrohabitat specialization than the result of geological events.Other examples of ecologically driven divergence on Hispaniola occurred between two sister bird genera, Xenoligea a montane taxon and Microligea a generalist taxon (Sly et al., 2011), and in a clade of Calisto butterflies, where the uplift of Cordillera Central potentially provided new ecological habitats that increased diversification rates (Matos-Maraví et al., 2014).This type of population isolation by environment (IBE) has been a major force of ecological speciation in insular plant genera (Huang et al., 2017), and it has been shown that an altitudinal gradient can drive speciation by ecological specialization in islands of the Mediterranean region (e.g.Corsica and Crete, Steinbauer et al., 2013).Moreover, in other Podocarpus species it has been postulated that past climatic events may have been the cause of restrictions to gene flow between populations, which is now detected by the presence of unique haplotidic variants in nearby populations (Mellick et al., 2012;Quiroga et al., 2012) or distant populations (Dantas et al., 2015), even in divergent lineages since the Pleistocene (Migliore et al., 2020).Although we do not rule out the hypothesis that climatic events may have played a role in the divergence of these lineages, new studies with a historical/climatic focus should be conducted to analyse this hypothesis.

| Management units for conservation
Podocarpus on Hispaniola as other members of the genus, wrongly commonly named pines around the world, usually occur in mixed forests as small populations and isolated individuals that may result in different types of rarity and genetic polymorphism in conifers (Premoli et al., 2001).To the best of our knowledge, this is the first study using genotyping-by-sequencing of the whole genome within the genus Podocarpus, and SNP polymorphisms of endangered P.buchii and P.hispaniolensis yielded elevated within-population genomic diversity.This is despite the low number of individuals analysed per population that in some cases reflected the existence of solitary individuals in the field such as in collection Sites 2 and 9.
Therefore, guidelines to protect them should focus at the individual and/or grove level.Moreover, not only Hispaniola Podocarpus had elevated heterozygosity but also all analysed populations had hundreds or in some cases thousands of private, that is unique, alleles.
In contrast, reduced within-population genetic diversity was found in the only conifer found in north-eastern Brazil Podocarpus sellowii using SSR and ISSR markers (Dantas et al., 2015).Thus, highly variable genomic markers are useful in describing polymorphism levels even of populations occurring in small numbers.In addition, the genetic structure clusters (see Section 3) and their geographic distribution: SB-b = P.buchii from Sierra Bahoruco, SN-b = P.buchii from Sierra Neiba, CC-b = P.buchii from Cordillera Central, CC-h = P. hispaniolensis from Cordillera Central and CS-h = P.hispaniolensis from Cordillera Septentrional (Figure1).Population sample sizes ranged from 4 to 10 (Table2).Population sizes are designated as N (e.g.NSB-b corresponds to P.buchii from Sierra Bahoruco).Population sizes during bottlenecks from founder events are designated as Nf, and population size of an unsampled ancestor is designated as NA.We designed eight competing demographic models to represent the four different evolutionary hypotheses (Figure2).Model 1 represents an ancestral population from Cordillera Central that gave rise to both Podocarpus species from where subsequently one species migrated south and the other north.Model 2 has identical divergence patterns to Model 1 but includes bottlenecks for each migration event.Model 3 proposes that the two species originated allopatrically in the two palaeoislands, reaching secondary contact in the Cordillera Central.Model 4 has identical topology to Model 3 but includes bottlenecks.Model 5 proposes a stepping-stone migration pattern from the south to the north.Model 6 has identical topology to Model 5 but with bottlenecks.Model 7 proposes a stepping-stone migration pattern from the north to the south.Model 8 has identical topology to Model 7 but with bottlenecks (Figure2).Bottlenecks in Models 2, 4, 6 and 8 were included after each colonization step to model founder effects.Additional information for the DIYABC analysis can be found in the Supplementary Material.

TA B L E 2
Genetic diversity statistics for the 11 Podocarpus collection sites on Hispaniola Island and K = 5 showed a clear separation between individuals from collection Sites 10 and 11 within Cordillera Septentrional, and those from Cordillera Central (Figure S1).For K = 2, individuals from Cordillera Central did not show any admixture but individuals from Cordillera Septentrional did, with one individual showing ca.90% association with the Cordillera Central cluster.For K = 5, two individuals from Cordillera Central and five from Cordillera Figure 4a,b).Both analyses recovered with high support (PP = 1.0 and BS = 100%) a clade that included all P.hispaniolensis sites.F I G U R E 2 Schematic diagram showing the eight historical demographic models compared in DIYABC for Podocarpus on Hispaniola.Models 2, 4, 6 and 8 have the same divergence pattern as Models 1, 3, 5 and 7 but with bottlenecks indicated with an asterisk.Models 1 and 2: Speciation in Cordillera Central and subsequent migration to the north and south.Models 3 and 4: Allopatric speciation in the two palaeo-islands followed by secondary contact in Cordillera Central.Models 5 and 6: Stepping-stone colonization from south to north.Models 7 and 8 Phylogenetic resolution and support were better for P.hispaniolensis sites than for P.buchii sites.The two southernmost P.hispaniolensis sites in Cordillera Central (5 and 6) were sister to a clade formed by the three northernmost P.hispaniolensis sites (9, 10 and 11).

F
Bar plots showing the genetic structure of both Podocarpus species on Hispaniola.Analyses using (a) all loci, (b) neutral loci and (c) loci under selection.Numbers at the top of the barplots indicate collection sites a) Maximum-likelihood (RAxML) phylogenetic tree of the 11 Podocarpus collection sites on Hispaniola Island.Numbers at nodes are bootstrap support values.(b) Bayesian (MrBayes) phylogenetic tree of the 11 collection sites.Numbers at nodes are posterior probability branch support values.Numbers in parenthesis indicate collection sites as in Table 1 1.0 and BS = 95%).The ML and BI analyses for the five populations provided poor resolution and support due to extremely short branches (trees not shown).
species and population divergence have been influenced by historical and current dispersal barriers.The genetic diversity of Podocarpus on Hispaniola was geographically structured by mountain ranges.This geographic structure, however, could not be attributed to Pliocene-Pleistocene marine incursions or the fusion of palaeo-islands because individuals did not form two genetic clusters or phylogenetic clades corresponding to the north and south palaeo-islands and because of the very recent estimated population divergence times.The genetic structure instead was likely influenced by contemporary topographic barriers of the dry valleys separating the different mountain ranges where Podocarpus grow.The demographic history of Podocarpus on Hispaniola was likely the product of a single colonization with subsequent dispersals to cloud forests in a stepping-stone manner with bottlenecks associated with each new mountain range migration.The phylogenetic reconstruction suggests a speciation scenario where P.hispaniolensis emerged from P.buchii.
Hispaniola formed three different evolutionary lineages.The genetic admixture observed in Podocarpus suggests the presence of gene flow across the different mountain ranges or remnant ancestral similarity given the life cycle of Podocarpus and the fact that gymnosperms tend to have lower evolution rates (De La Torre et al., 2017).Other studies on Hispaniola have shown gene flow among bird populations occupying different mountain ranges,suggesting over-land dispersal between the south and north palaeo-islands(Sly et al., 2011).Despite the presence of dry valleys that might act as barriers for Podocarpus, gene flow might have homogenized some populations; for example, the genetic cluster in blue (Figure3a) is present in all populations.The different sets of loci used in the combined versus individual genetic structure analyses could explain the lack of admixture for P.hispaniolensis in the former but not the latter.
each pure Podocarpus taxon was clearly identified and amongpopulation divergence, although relatively low, was clustered by location which can be used to design conservation actions.Similar geographically controlled genetic structure was documented in Podocarpus parlatorei from the subtropical Yungas of northern Argentina and southern Bolivia(Quiroga et al., 2012).Thus, conservation strategies for endangered Podocarpus should follow a multi-population approach where distinct population clusters of each species should be locally protected throughout Hispaniola.Southern populations of P.buchii at Sierra Bahoruco (south palaeoisland) and Sierra Neiba (south of north palaeo-island) although geographically close are genetically different and thus need to be conserved as distinct management units.In this area, only collection Site 3 at Loma del Quince is located into a National Park Sierra de Neiba (FigureS5).Similarly, nearly sympatric high-elevation P.buchii and low-elevation P.hispaniolensis from Cordillera Central, although they had lower number of unique alleles compared to allopatric southernmost P.buchii and northernmost P.hispaniolensis, were clearly identifiable by STRUCTURE and phylogenetic analyses.Therefore, Cordillera Central seems to reflect how the distinct species have differentiated by ecological divergence.Yet,

of genetic diversity and structure of Podocarpus on Hispaniola
Ion Proton sequencing generated 32.36 GB of raw data containing ca. 102.6 million reads.After demultiplexing, quality filtering and discarding reads with ambiguous barcodes, a total of ca.93.7 million reads (91.3% of initial number of reads) were kept for further analyses.The number of reads per sample ranged from ca. 258 K to ca. 10.3 M, with an average number of ca.2.7 M reads per sample.Read quality inspection in FastQC reported a Phred score average of 24, and GC content of 57%.The catalog built from Stacks pro-Population genetic statistics of observed (H o ) and expected heterozygosity (H e ) and nucleotide polymorphism diversity (π) showed that P.buchii had higher genetic diversity than P.hispaniolensis (Table2).This result was consistent with the higher number of genetic clusters found in P.buchii compared with P.hispaniolensis.When individuals were grouped into five populations, SN-b had the highest H e followed by SB-b, while the northernmost population (CS-h) had the lowest H e .CC-b and SB-b showed the highest H o and π, and CS-h showed the lowest π (TableS2).When individuals were grouped by collection site, H e and π were highest for the two Sierra Neiba sites (3 and 4) with P.buchii, while the lowest values were observed on collection Site 6 for P.hispaniolensis in Cordillera Central.Site 2 (P.buchii in SB) had highest H o , while Site 5 (P.hispaniolensis in CC) had lowestH o (Table