Introduction

Genomic tools are now available for many organisms from all over the tree of life and enable ecologists and evolutionary biologists to perform fundamental tasks such as inferring recent and historical demographic trends of populations in response to environmental changes (Aitken et al. 2008; Neale and Kremer 2011). However, the large size and the complexity of many tree genomes render the development of species-specific genomic resources for this group particularly challenging. For example, many gymnosperms have huge genomes with a high proportion (up to 60%) of repetitive content, and polyploidy is common among angiosperms (Ellegren 2014; Soltis et al. 2015). Although the advent of increasingly fast and affordable next-generation sequencing (NGS) technologies has partly solved these obstacles, genomic research on forest trees still remains largely focused on a few species of economic importance, whereas most ecologically relevant tree species still lack basic genomic resources. In addition, a strong geographic bias towards temperate tree species constrains our ability to decipher the long-lasting evolutionary impact of climate change across forested lands worldwide, a task of utmost importance to forecast the response of forests to emerging mega-disturbances (Millar and Stephenson 2015).

Forests and woodlands have been deforested for centuries in favor of agriculture and farming worldwide, but rural abandonment has also prompted the expansion of native forest species from isolated remnant patches in the past decades, particularly in the Northern Hemisphere (Trumbore et al. 2015). Tracking the demographic expansion of formerly isolated forests and tailing the advance of invasive species are two major challenges for ecologists interested in maintaining the biodiversity of managed landscapes and the ecosystem services they provide. Accomplishing those challenges requires obtaining reliable demographic estimates to design scientifically informed management plans. Estimates such as the effective population size (Ne) can be inferred from genetic data of extensively genotyped populations, as already proven in endangered animal populations with fluctuating demographic trends (McCoy et al. 2013). Yet, unlike forestry and agronomic studies, ecological studies typically rely on a limited number (~ 10–20) of simple sequence repeat (SSR) markers that often prove insufficient for gaining reliable demographic inferences from genetic data. These studies could enormously benefit from the incorporation of single-nucleotide polymorphisms (SNPs) because of the following: (i) these provide less biased estimates of gene flow (Singh et al. 2013) and perform better in inferring the timing of recent and past demographic bottlenecks (McCoy et al. 2013); (ii) SNP calling is less prone to ambiguous visual scoring compared to SSR calling, although it is also strongly dependent on the parameters set in the bioinformatic pipelines; and (iii) SNPs allow not only demographic inferences but can potentially also characterize functional genetic variation (Eckert et al. 2010). Therefore, the development of SNP markers and the combination of different molecular marker types will improve our understanding of the evolutionary trends of tree species inhabiting managed landscapes, including species with large and complex genomes.

Several Juniperus species have increased their population sizes and geographical distribution during the past decades driven by land-use changes, namely rural abandonment, and increasingly arid conditions (Van Auken 2008). For example, fragmented populations of Juniperus phoenicea have rapidly expanded across a matrix of surrounding Mediterranean shrublands (García et al. 2014); former remnant patches of J. thurifera now have developed into extensive savanna-like forests in former agricultural landscapes of the central Iberian Peninsula (Escribano-Ávila et al. 2012; García-Cervigon et al. 2016), and J. monosperma has become the dominant species of the piñon-juniper ecotone in northern New Mexico following an extreme drought event in the 1950s (Allen and Breshears 1998). At a very different temporal scale, phylogeographical studies have retraced the postglacial expansions of juniper taxa that represent today the major tree species of many central Asian highlands and mountain ranges (Opgenoorth et al. 2010). All these junipers act as foundation species structuring their community and creating stable local conditions for other plant species (Whitham et al. 2006), and therefore, they are fundamental to forecast the response of plant communities to recent climatic trends. The genus Juniperus is moreover central to evolutionary studies examining past tree distributional ranges driven by long-term geological and climate changes (Mao et al. 2010, 2012). Its complex evolutionary history still puzzles evolutionary biologists whose studies on deciphering the evolutionary history of this complex group would benefit from a greater availability of transferable genetic markers among closely related species (Li et al. 2012).

Previous studies have reported the development of a handful of SSR markers for closely related juniper species that comprises the section Sabina (Mao et al. 2010), such as J. sabina (Geng et al. 2017), J. thurifera (Teixeira et al. 2014), and J. phoenicea (Molecular Ecology Resources Primer Development Consortium et al. 2013) but they resulted insufficient to perform demographic inference and parentage analyses. Here, we aim to contrast the performance of two types of molecular markers on reliably inferring population demographic parameters in long-lived tree species that inhabit managed landscapes. To attain that goal, we designed 19 SSR and 147 SNP markers for Phoenician juniper (Juniperus phoenicea ssp. turbinata), a woodland-forming conifer species that is widely distributed across the Mediterranean Basin and the Macaronesian archipelagos (Adams 2011). Phoenician juniper woodlands have been severely damaged for centuries by rural activities but are currently expanding in protected areas (García et al. 2014) and on abandoned agricultural lands (Bello-Rodríguez et al. 2016). We derive a series of population genetic diversity estimates for a well-known expanding population (García et al. 2014) using the two marker sets either independently or in combination to compare their respective performance. Then, we test the transferability of the new SNP markers to the closely related tetraploid species J. thurifera. Finally, we search for functionally equivalent genes based on the extensive genome sequences available for two gymnosperm model species, Pinus pinaster and Pinus taeda in order to facilitate prospective landscape genomic studies with juniper species.

Material and methods

Plant material and DNA extraction

Juniperus phoenicea ssp. turbinata grows primarily on stabilized coastal dunes and rocky seashores (Adams 2011). The species is monoecious, wind-pollinated and produces berry-like cones (arcestides) that are dispersed by frugivorous vertebrates. Populations greatly differ in their conservation state ranging from well-preserved stands (mainly in protected or inaccessible areas) to chronically fragmented and jeopardized sites (Bello-Rodríguez et al. 2016). We used plant material from a 1.2-ha permanent study plot previously established in an expanding juniper population that occupies coastal stabilized dunes in the Reserva Biológica de Doñana (lat. 37.017085, long. − 6.554601; Huelva province, Spain) (García et al. 2014). For the present study, we randomly chose 83 individuals from this long-term study plot to be genotyped. We also included 11 samples of Juniperus thurifera collected in the Parque Natural del Alto Tajo (lat. 41.063114, long. 2.198333; Guadalajara province, Spain) to test for the cross-species transferability of the newly developed SNP markers.

Nuclear DNA was extracted from 5 to 10 mg of dry leaf tissue per individual. This material was placed into a 2 ml screwed-top tube with two steel beads (4 mm in diameter). Tubes were frozen in liquid nitrogen before grinding at 30 Hz for 2 min and 30 s with a Mixer Mill MM300 (Retsch, Germany). When tissue failed to grind into powder, we applied a second grinding cycle. For DNA isolation, we used a genomic DNA plant tissue kit NucleoSpin Plant II (Macherey-Nagel, Duren) and followed the manufacturer’s instructions.

SSR marker development and genotyping

Marker development

An equimolar mixture of genomic DNA from 15 individuals sampled from the same focal population with a final concentration of 16.5 ng/μL measured with a fluorometer (Invitrogen Qubit 4, Life Technologies, Singapore) and an average A260/A280 ratio ca. 1.8 was sent to Ecogenics (www.ecogenics.ch) to proceed with high-throughput microsatellite isolation following Malausa et al. (2011). Size-selected fragments from genomic DNA were enriched for SSR content by using magnetic streptavidin beads and biotin-labeled CT and GT repeat oligonucleotides. The SSR-enriched library was analyzed on an Illumina MiSeq platform using the Nano 2 × 250 v2 format. After the assembly using ABySS software (Simpson et al. 2009), SSRs were detected using MISA Perl script (Thiel et al. 2003) that yielded 9313 contigs or singlets containing a microsatellite insert with a tetra- or a trinucleotide of at least 6 repeat units or a dinucleotide of at least 10 repeat units. Suitable primer design was possible in 4451 SSR candidates. Primers were designed with Primer3 using standard settings (Rozen and Skaletsky 1999). Initially, we selected 96 loci with lengths between 7 and 21 repeats following van Asch et al. (2010) and with amplicon sizes ranging from 100 to 400 bp to facilitate multiplex setting. Among them, 21 loci amplified successfully and 19 loci showed polymorphism with a subset of 15 individuals randomly sampled in our study area. Previous studies have shown a similar percentage of success in developing SSRs for non-model tree species (De Bellis et al. 2016; Schoebel et al. 2013).

Genotyping

We genotyped all 83 J. phoenicea ssp. turbinata individuals at the 21 selected SSR loci. For this purpose, we optimized multiplex PCR conditions for the set of primers following a three-primer approach (Schuelke 2000). Multiplexed PCRs were performed in a 20-μL final volume containing 1× buffer [67 mM Tris–HCl pH 8.8, 16 mM (NH4)2SO4, 0.01% Tween-20], 3.5 mM MgCl2, 0.01% BSA (Roche Diagnostics, Germany), 0.25 mM dNTP, 0.5 U Taq DNA polymerase (Bioline, UK), 2 μL of the primer premix, and 5 μL (50 ng) of genomic DNA. The primer premix contained 0.5 μM of the primer M13 and a final concentration of the forward and reverse primers according to Table SM1. We used Multiplex Manager 1.2 (Holleley and Geerts 2009) to compute the temperature of annealing (Table SM1) and the complementary threshold (the maximum number of AT or CG matches for any two primers within a multiplex reaction), which was set to seven. Amplified fragments were analyzed on an ABI 3130xl Genetic Analyzer and sized using GeneMapper 4.0 (Applied Biosystems, USA) and LIZ 500 size standard.

SNP development and genotyping

Double-digest restriction-site-associated DNA library preparation and sequencing

Nine of the 83 individuals were used for double-digest restriction-site-associated DNA (ddRAD) library preparation (Peterson et al. 2012) and one sample was duplicate. ddRAD library preparation protocol followed the methods described by Pukk et al. (2015) with minor modifications, including purification by solid-phase reversible immobilization (SPRI) bead solution (CleanNA, Netherlands) after each step. In short, 650 ng of DNA was digested for 3 h at 37 °C with two rare-cutting restriction enzymes, 10 U of AseI (restriction site 5′ ATTAAT 3′) and PstI (restriction site 5′ CTGCAG 3′) (New England Biolabs, USA). The ligation consisted of 0.04 μM of P1-AseI and A-PstI adapters, which were added to 8 μL of restriction reaction together with 0.5 mM of ATP, 1× of T4 DNA Ligase Buffer, and 800 U of T4 DNA Ligase (New England Biolabs). A barcode (Ion Xpress 1-9) was added to the A-PstI adapter to identify each sample a posteriori. The 20-μL ligation reaction was carried out at 22 °C for 2 h and heat-inactivated for 11 min at 65 °C before cooling (1 °C per minute). Libraries were then quantified with the Ion Library TaqMan Quantitation Kit (Thermo Fisher Scientific) before equimolar library pooling. The pool was loaded on an automated size-selection system (Pippin Prep—Sage Science, USA) with a 2% agarose cartridge to extract DNA fragments of 290–310 bp. The sized pool (30 μL) was amplified in a 100-μL reaction containing 1× Q5 High Fidelity PCR Master mix and 0.6 μM of Ion Torrent primers A and P1 (New England Biolabs). PCR consisted of 98 °C for 30 s followed by 12 cycles of 98 °C for 10 s, 58 °C for 30 s, and 65 °C for 30 s. The quality and quantity of the pool were measured on a Bioanalyzer 2100 using High Sensitivity DNA kit (Agilent Technologies, USA) and the final library was diluted to 10 pM before sequencing on an Ion Torrent PROTON (Thermo Fisher Scientific).

Quality control and SNP discovery

Raw sequences were demultiplexed based on their barcodes and quality filtered using the default settings of the Ion Torrent BaseCaller (> Q16 with a window size of 30 bases). Filtered reads were used for de novo identification of putative SNPs using three programs implemented in STACKS (Catchen et al. 2013). First, we applied process_radtags to trim all reads to 200 bp. Then, we used denovo_map to build the catalog of loci and call SNPs with the following parameters: minimum number of identical reads m = 15, number of mismatches allowed between loci when processing a single individual M = 2, number of mismatches allowed between loci when building the catalog n = 3. Finally, populations generated the FASTA and VCF files used for downstream analysis. Putative SNPs were filtered out using different criteria: (i) SNP called for all nine samples; (ii) unique SNP within each stack; (iii) an identical genotype for the two technical replicates; and (iv) three genotypic classes identified and SNP not present in the first or last 20 bases of the sequence.

SNP selection for genotyping, SNP calling, and genotyping of individuals

A total of 187 candidate SNPs were submitted for assay design using the MassARRAY® Assay Designer version 4.0.0.2 (Agena Biosciences). Four multiplexes of 156 SNPs (three 40-plex and one 36-plex) were designed for the SNP genotyping, which was performed using the iPLEX Gold chemistry following Gabriel et al. (2009) on a MassArray System (Agena Biosciences, USA). Data analysis was completed using MassARRAY Typer Analyzer 4.0.26.75 (Agena Biosciences). We filtered out all monomorphic SNPs, loci with weak or ambiguous signal (i.e., displaying more than three clusters of genotypes or unclear cluster delimitation) and loci with > 6% missing data.

Functional annotation of SNP markers

An orthology analysis was performed to search for functional annotation of the 187 sequences of the candidate SNPs using Blastall 2.2.15 (Zhang et al. 2000). The BlastN algorithm was used to align the flanking sequence of the SNPs with the transcriptome of Pinus pinaster (Canales et al. 2014, SustainPine v3.0, http://www.scbi.uma.es/sustainpinedb/). For Pinus taeda, the BlastX algorithm was used with the ConGenIE database (http://congenie.org/). The E value cut-off was set at 10−10. We considered only the best blast hit for biological interpretation.

Evaluation of SSR and SNP markers to perform population genetic analyses

We first examined the quality of the 21 SSR markers based on the analysis of the multilocus genotypes of 15 individuals run on an ABI3730 (Applied Biosystems, California) (Table SM1). Specifically, we evaluated the genetic correlation among SSR markers by estimating the Hardy-Weinberg equilibrium (HWE, Table SM2) and by testing for genotypic linkage disequilibrium as implemented in GENEPOP (Rousset 2008) (Table SM3). Additionally, PopGenReport version 3.0 served to infer the frequency of null alleles based on the Brookfield estimator (Brookfield 1996) (Table SM1). Finally, we gauged the polymorphism of each marker by assessing (Table SM1): (i) the number of alleles (A); (ii) the effective number of alleles (Ae), an estimate of the number of equally frequent alleles in an ideal population, Ae is of interest for comparison of allelic diversity across loci with diverse allele frequency distributions; (iii) unbiased expected and observed heterozygosity (uHe and Ho, respectively); (iv) the polymorphic information content (PIC), a measure of the informativeness of a genetic marker (Botstein et al. 1980); (v) the probability of identity (PI, i.e., the average probability that two unrelated individuals drawn from the same randomly mating population will have the same multilocus genotype by chance); and (vi) the probability of exclusion (PE, the probability of excluding a putative parent pair). We used PICcalc (Nagy et al. 2012) to obtain the PIC per locus and we implemented GenAlEx 6.502 (Peakall and Smouse 2012) to obtain all other estimates.

Secondly, we performed population genetic analyses on the 83 genotyped J. phoenicea ssp. turbinata individuals for each marker type independently and for both types combined. We recorded the overall number of different alleles (allelic richness, A) and estimated the unbiased expected and observed heterozygosity (uHe and Ho) as implemented in GenAlEx. We also calculated unbiased multilocus estimates of population inbreeding coefficient (FIS) as implemented in INEST 2.2 (Chybicki et al. 2011) that has proved to be robust to the presence of null alleles (Campagne et al. 2012). As a measure of the usefulness of each set of markers for assignment studies, we measured average polymorphic information content (PIC), the overall parentage exclusion probability, the informativeness for inferring relationships (IR), and relatedness (Ir) as implemented in Coancestry (Wang 2011). The suitability of each set of markers to provide reliable demographic inferences was measured by estimating the contemporary effective population size (Ne) based on the linkage disequilibrium method as implemented in NeEstimator v2 (Do et al. 2014). We estimated Ne for the complete set of individuals successfully genotyped. We also evaluated the ability of each set of markers in estimating Ne of small-sized populations. To that end, we sampled genotypes randomly from our data set (N = 83) to create 100 populations of sizes N = 10, 25, and 50, respectively. Then, we estimated Ne as above and recorded the average Ne and average confidence intervals over 100 populations for each size class. Finally, we tested the ability of each set of markers in detecting a recent bottleneck based on a sign test as implemented in Bottleneck (Cornuet and Luikart 1996).

Results

SSR and SNP development and selection

Two of the 21 new SSR markers (Junpho_017898 and Junpho_068482) were monomorphic and hence were discarded from all further analyses (Table SM1). Five of the remaining 19 loci did not meet HWE (Table SM2) and one pair of loci showed statistical evidence of linkage disequilibrium (LD) after applying the B–Y correction (Narum 2006) for multiple tests (Table SM3). A total of 102 million filtered reads (median size = 229 bp) were generated for the nine samples and 49,457 putative SNPs were identified. The SNP calling error rate, based on the replicated sample, was 5.5% for the 5661 loci for which all nine samples were available (19.6% for all SNPs for which at least the two technical replicates were called). The multi-criteria filtering generated a total of 187 SNPs retained for multiplex design. We discarded 40 SNPs because they did not properly amplify. We finally used a set of 147 SNPs for the population genetics analyses. Our final data set for these analyses contained ca. 2% (range 0–6%) of missing data per locus.

Estimates of genetic diversity

The overall number of SNP alleles more than tripled the number of SSR alleles (Table 1). Values of mean observed heterozygosity (Ho) across loci were similar for both sets of markers, whereas unbiased expected heterozygosity (uHe) was higher for SSRs than for SNPs. The SSR-based estimate of the mean population inbreeding coefficient (Fis) exceeded the estimate obtained with SNPs. PIC values showed that SSR markers were on average more informative than SNP markers and, as a result, SSRs yielded higher values of IR, Ir, and parentage exclusion probability (Table 1). In turn, SSRs failed to detect a recent bottleneck that was successfully detected with SNPs either alone or in combination with SSRs. In addition, Ne estimates based on SNPs or the combination of SNPs and SSRs were more reliable both for the full data set (N = 83) as well as for the smaller populations subsampled from the original data set (N = 10, 25, and 50) (Table 1). For small population sizes (N = 10 and N = 25), SSRs frequently failed to provide a reliable Ne estimate. On the contrary, the SNPs, alone or in combination with SSRs, always successfully estimated Ne when we subsampled as low as N = 25 individuals from the original data set.

Table 1 Genetic diversity estimates for N = 83 juniper (Juniperus phoenicea ssp. turbinata) individuals based on 19 simple sequence repeat (SSR) markers (left), 147 SNP markers (center), and both types of markers (right). For each set of markers, we report the number of polymorphic markers attained, the total number of different alleles, the range of the number of alleles per locus, the observed heterozygosity (Ho), the unbiased expected heterozygosity (He), the posterior mean estimate of the population inbreeding that takes into account the presence of null alleles (Fis), the mean polymorphic information content (PIC), the mean informativeness for inferring relationships (IR) and relatedness (Ir) according to Wang (2002), and the probability associated to a sign test to identify a recent bottleneck. We also report estimates of the contemporary effective population size (Ne) and their associated confidence intervals based on the linkage disequilibrium method for three different sample sizes: N = 83 (the full data set); N = 10 (subset of 10 individuals randomly sampled 100 times from the full data set); N = 25 (subset of 25 individuals randomly sampled 100 times from the full data set); N = 50 (subset of 50 individuals randomly sampled 100 times from the full data set). Estimated Ne values for N = 10, N = 25, and N = 50 are averaged for 100 simulated populations. Percentages indicate the number of times that Ne was successfully estimated over 100 simulated populations. Note that for both sets of markers, the estimated Ne sometimes exceeded the simulated population size (N) because the samples are integrated within a larger remnant forest patch composed of few hundreds of trees (N ≈ 750 individuals)
Full size table

Transferability and functional annotation of SNP markers

A total of 112 out of the 156 tested SNP loci amplified in the sister species J. thurifera. We detected only five polymorphic loci, possibly due to some extent to the relatively small number of individuals tested. Only five (P. pinaster) and four (P. taeda) loci, respectively, could be annotated with strong E values (6e−11 to 6e−30) and identities (75 to 88.7%) (Table 2). Two loci (1334_130 and 23712_40) emerged in both transcriptomes, whereas the remaining 182 loci (96.3%) did not match any sequence.

Table 2 Identification of putative functional genes after applying an orthology analysis that compares the sequences of the new set of SNP markers developed for Juniperus phoenicea ssp. turbinata with the sequences obtained based on the transcriptome of Pinus pinaster (http://www.scbi.uma.es/sustainpinedb) and the ConGenIE database for Pinus taeda (http://congenie.org/)
Full size table

Discussion

The Aichi Biodiversity Targets include the preservation of the genetic diversity as one strategic goal, yet most non-model, but ecologically relevant, tree species lack genetic and genomic resources that allow researchers forecasting the fate of remnant populations in a changing world (Neale and Kremer 2011). As a result, the majority of recovery plans for endangered plant species overlook genetic factors as drivers of population extinction (Pierson et al. 2016), in spite of the ample existing evidence that underpins the role of genetic erosion and inbreeding in accelerating demographic decline (Frankham 2005). Therefore, by providing new genomic resources for a foundation tree species (J. phoeniceas spp. turbinata), our study can contribute to advancing both basic and applied ecological and evolutionary research (Mao et al. 2010).

Many juniper species are foundation species of semi-arid and arid ecosystems across the Northern Hemisphere, both in warm low-elevation and cold high-mountain habitats (Mao et al. 2010). To date, only 11 SSRs have been described for juniper species (Molecular Ecology Resources Primer Development Consortium et al. 2013), but these were not polymorphic enough to perform demographic inference or parentage analyses. The new set of SSR and SNP markers proved fully suitable to estimate genetic diversity, perform assignment or parentage analyses, and obtain reliable demographic inference. The issue of whether a small number of SSRs would do better in performing population genetic analysis than a large set of genome-wide distributed SNPs remains contentious. Our results show that mean values and variance of SSR-based estimates of genetic diversity (uHe) exceeded those based on SNPs, as found in the previous studies (Fischer et al. 2017; Glover et al. 2010; Hamblin et al. 2007). This result confirms that estimates of genetic diversity typically show a wide variation among SSR markers (in terms of uHe), which make SNPs more reliable when it comes to infer genome-wide genetic diversity, particularly when at least a few thousands of random SNPs are available (Fischer et al. 2017). The average PIC values show that both sets of markers, taken alone, still have a relatively moderate information content. However, the higher values of IR, Ir, and parentage exclusion probability obtained for SSRs (Table 1) indicate that they would be the marker of choice for assignment techniques such as parentage, maternity, and paternity analyses. On the other hand, the SNP marker set clearly outperformed the SSRs in obtaining reliable demographic inferences, both in gaining robust estimates of Ne and in identifying a recent bottleneck event. This finding concurs with a previous study showing that SSRs typically fail to identify recent genetic bottlenecks (Peery et al. 2012). SNPs yielded reliable Ne estimates even for small-sized samples: for N = 25, SNPs provided robust Ne values in 100% of the simulated populations, whereas SSRs only yielded reliable Ne estimates in 47%. This result suggests that, given a comparable budget to develop each set of markers (Hodel et al. 2016), SNPs are the marker type of choice for performing demographic inferences, particularly in small populations of non-model species. Note that this study is based on a remnant age-structured forest patch that has undergone successive cycles of population contraction and expansion. Therefore, genetic drift, founder effects, and non-random mating patterns are expected to impede allele frequencies to meet the Hardy-Weinberg proportions (HWP), which is the case for five out of the 19 SSRs. Similarly, these natural processes increase the probability that different markers show linkage disequilibrium (LD), as we found for one of 171 pairwise comparisons. As pointed out by Waples (2015), once different sources of genotyping errors have been ruled out and the frequency of null alleles remains low (< 10%), the presence of markers that do not meet the HWP and LD in natural populations could be the result of biological processes rather than genotyping errors. Although our results regarding the suitability of SSRs and SNPs for different types of evolutionary analyses might be influenced by the particular characteristics of the sample population, our results concur with previous studies concluding that the relative utility of each set of markers depends on the main goals of the study (Emanuelli et al. 2013; Singh et al. 2013; Yang et al. 2011), the availability of previous genetic resources, and the number of individuals sampled (Hodel et al. 2016; Nazareno et al. 2017).

In conclusion, if the main goal of the study is to describe the distribution of genetic variation across the landscape (i.e., landscape genetics) or ascertaining parentage and pedigree relationships (for example in mating system studies), a handful of SSR markers represent a cost-efficient tool that can provide reasonably good results (Hodel et al. 2016). However, if the main interest focuses on obtaining reliable demographic inferences, then SNPs would be the marker of choice, particularly when the study entails small-sized populations such as in conservation and recovery plans Waters et al. (2013). Studies would of course benefit from applying both types of markers whenever possible (Graudal et al. 2014; Olsson et al. 2016). The transferability of the new set of SNPs to J. thurifera suggests that SNP markers could potentially be applied to address ecological and evolutionary studies entailing other closely related species of this key genus (Mao et al. 2010). Yet, we found a low level of polymorphisms in J. thurifera, which encourages the use of increasingly cost-effective de novo development of species-specific SNPs and sequence data to address evolutionary studies across taxa in non-model species.

The functional annotation of SNPs links population genetics and population genomics as it allows to dissect the contribution of neutral processes (such as gene flow) and selective processes (such as local adaptation) in determining current patterns of genetic variation across heterogeneous gradients (Eckert et al. 2010). Non-model species typically lack a reference genome, but an orthology analysis based on gymnosperm model species allowed us to determine whether any of our SNPs is located at a functionally relevant position in the genome. Given the large genome size of most gymnosperm species, finding SNPs associated to genes that code for ecologically relevant proteins is challenging and, as a result, we only identified genes with generic functions, such as nucleic acid–binding proteins. Moreover, the Sequenom sequences we used were very short (80–120 bp), a factor that increases the likelihood of failure when performing functional annotation. SNP development from transcriptomes would certainly increase the chances to find polymorphism linked to genes that mediate plant responses to environmental factors, such as prolonged droughts (Neale and Ingvarsson 2008). In a context of environmental change, identifying locally adapted genotypes is of utmost importance to pinpoint source populations in breeding programs or to design assisted migration plans. By doing so, managers could potentially move best-adapted genotypes to locations where environmental conditions are expected to shift in the near future, for example by becoming more arid (Aitken et al. 2008; Jordan et al. 2017). Lastly, the possibility of transferring genomic resources among closely related species expands the application of molecular markers to address long-lasting evolutionary studies to identify hybridization events, depict diversification patterns, or compare genetic diversity among closely related species to infer the ecological and evolutionary factors that have shaped current patterns of genetic diversity among forest species with complex genomes (Petit and Hampe 2006).