Pinus subsection Cembroides Engelm. is a clade of North American pinyon pines occurring in arid or semi-arid environments from the southwestern United States to south central Mexico (Critchfield and Little 1966). The species of this subsection are small to medium-sized trees or shrubs characterized by secondary leaves with deciduous fascicle sheaths. Except for P. rzedowskii Madrigal, Caball.M., all taxa have enlarged seeds with a thickened sclerotesta (Madrigal and Caballero 1969) and are functionally wingless. The morphological characters that have been used in species identification in this group include the number of needles per fascicle, needle length and width, the distribution of stomata on the adaxial and abaxial leaf surfaces, and cone morphology (Malusa 1992; Farjon and Styles 1997). Pinus subsection Cembroides is classified in section Parrya Mayr together with two other North American subsections, Balfourianae Engelm. and Nelsoniae Burgh (Gernandt et al. 2005). Pinus subsection Nelsoniae is monotypic, with the pinyon pine, P. nelsonii Shaw distinguished by persistent fascicle sheaths and connate needles (Little and Critchfield 1969; Gernandt et al. 2001). Section Parrya is classified together with section Quinquefoliae Duha-mel in subgenus Strobus Lemmon (Gernandt et al. 2005).

Seeds of P. cembroides subsp. cembroides Zucc. and P. edulis Engelm. are important sources of food in Mexico and the United States (Lanner 1981; Farjon and Styles 1997). The high nutritive value in seeds encourages interactions among pinyon pines, rodents, and corvid birds. The needles are sometimes used for medical treatments (Lanner 1981). The International Union for Conservation of Nature (IUCN 2017) lists five pinyon pine taxa as vulnerable or endangered and the Mexican government lists nine as protected (SEMARNAT 2010; Table 1).

Pinus subsection Cembroides has been the focus of several phylogenetic studies (Malusa 1992; Farjon and Styles 1997; Gernandt et al. 2001, 2003; Parks et al. 2012; Flores-Rentería et al. 2013) and the number of recognized species and infraspecific taxa differs in recent works (Malusa 1992; Eckenwalder 2009; Farjon and Filer 2013; Table 2). Phylogenetic results in the subsection have varied due to the use of different morphological, ecological, and molecular characters and differences in sampling (Malusa 1992; Farjon and Styles 1997; Gernandt et al. 2001, 2003, 2005; Syring et al. 2005). Based on a cladistic analysis of morphological and ecological characters, Malusa (1992) divided subsection Cembroides into a group of eight species with small seed cones and a second group of four species with large cones. In a restriction site study of noncoding plastid DNA, Pérez de la Rosa et al. (1995) recovered P. nelsonii as separate from subsection Cembroides. A cladistic analysis of morphological characters in Neotropical species of Pinus subgenus Strobus by Farjon and Styles (1997) recovered the small-coned pinyons as monophyletic, and the four large-cones species as paraphyletic to subsection Strobus. Two of the large-cone species, P. nelsonii and P. pinceana Gordon, formed a clade, leading the authors to classify them together in subsection Nelsoniae. In contrast, Price et al. (1998) classified P. nelsonii and P. pinceana in subsection Cembroides and P. rzedowskii in the monotypic subsection Rzedowskianae Carvajal.

More recently, phylogenetic relationships of pinyon pines have been inferred from sequences of plastid and nuclear ribosomal DNA (Gernandt et al. 2001, 2003, 2005; Parks et al. 2012; Ortiz-Medrano et al. 2016). Gernandt et al. (2001) reported phylogenetic analyses of the ITS region of nrDNA. The authors found that divergent copies of the ITS region in individuals of the same species do not group together. Nevertheless, results from the ITS region study and two subsequent phylogenetic studies using plastid DNA sequences corroborated the separation of P. nelsonii from the other pinyon pines. Plastid DNA studies also have suggested that P. johannis M.F.Robert and P. discolor D.K.Bailey & Hawksw., are not infraspecific taxa of P. cembroides but instead close relatives of P. culminicola Andresen & Beaman (Gernandt et al. 2003, 2005; Parks et al. 2012; Ortiz-Medrano et al. 2016).

Table 1.

Taxa classified in Pinus subsection Cembroides with geographic distribution and conservation risk category. * Subject to special protection, and ** Listed as endangered by the Mexican government (SEMARNAT 2010). † Vulnerable and †† Endangered by the IUCN (2017).

Some, but not all the results from phylogenetic studies have been followed in subsequent taxonomic treatments (e.g. Eckenwalder 2009; Debreczy and Rácz 2011; Farjon and Filer 2013). Farjon and Filer (2013) recognized that P. nelsonii and P. pinceana are not closely related, but that P. nelsonii belongs to a more distant group, and they recognized that P. pinceana is closely related to P. maximartinezii Rzed. Farjon and Filer (2013) treated P. johannis as a subspecies of P. cembroides although this has been contradicted by phylogenetic analysis of plastid DNA. Another plastid DNA study left in question whether the single-leaf pinyon pines are monophyletic, recovering P. californiarum D.K.Bailey as sister to P. edulis rather than to P. monophylla Torr. & Frém. (Gernandt et al. 2007). Pinus californiarum has been treated as a synonym of P. monophylla for sharing a single needle, or separated as P. californiarum or P. monophylla var. californiarum (D.K.Bailey) Silba based principally on differences in the number of resin canals, number of stomatal lines, and diameter of the needle (Silba 1990; Farjon and Styles 1997; Price et al. 1998).

Table 2.

Two classifications of the species recognized in Pinus subsection Cembroides.

Introgressive hybridization and gene flow have been reported in species of Pinus subsection Cembroides (e.g. Mirov 1967; Lanner 1974a, 1974b; Lanner and Phillips 1992; Malusa 1992). Some populations of P. edulis (predominantly two needles per fascicle) and P. monophylla (predominantly single-needled) occur in sympatry (Fig. 1) in the eastern Great Basin where trees of P. edulis with both single needles and two needles per fascicle have been observed (Lanner 1974a). In a study of natural hybridization in pinyon pines in northwestern Arizona, Lanner and Phillips (1992) analyzed variation in morphological characters over different years. Based on differences in the frequency of needle and resin canal numbers at 22 sites, they concluded that bidirectional introgression was occurring between P. edulis and P. monophylla, and proposed that overlapping phenology, wind-dispersed pollen, and weak premating barriers were responsible.

Lanner (1974b) proposed that P. quadrifolia Parl. ex Sudw. originated from hybridization between P. monophylla (treated here as P. californiarium; predominantly single-needled) and a previously undescribed species with five needles per fascicle, P. juarezensis Lanner (treated here as a synonym of P. quadrifolia). According to Lanner (1974b), this would explain extreme needle number variation in P. quadrifolia, a characteristic that is frequently observed in pine artificial hybrids when parents differ in this character (Keng and Little 1961). The geographical distributions of the taxa overlap broadly (Fig. 1), with sympatric populations common in Baja California (e.g. in the Sierra Juárez), suggesting that putative hybrids of several types co-exist in a hybrid swarm (Lanner 1974b). Farjon and Styles (1997) reported that pollen dispersal occurs in April and May in P. monophylla and in March and April in P. quadrifolia, which would allow for interspecific gene flow. Nonetheless, the proposal to recognize P. juarezensis as one of the parental species has not been widely accepted. Pinus juarezensis was considered a synonym of P. quadrifolia by subsequent authors (Farjon and Styles 1997; Gernandt et al. 2003; Eckenwalder 2009). Plastid DNA from three samples of P. quadrifolia (two of which were from the type locality of P. juarezensis) grouped together and formed a sister group to P. monophylla from California, indicating that P. quadrifolia has not captured pollen of P. californiarum (Gernandt et al. 2003).

Fig. 1.

Map showing the distribution of Pinus section Parrya based on Critchfield and Little (1966). Species records are based on herbarium collections and data from field studies. The individual sample sites are indicated by triangles (Appendix 1). Colors represent North American subsections that belong to section Parrya.

In his description of P. remota, Little (1968) suggested the possibility that in the past its populations had been in contact with those of P. edulis, permitting introgressive hybridization. The geographic distribution of P. remota and P. edulis may have been more extensive in the past (Late Quaternary), bringing the species into contact, maybe in the Chihuahuan Desert (Van Devender 1990). Additionally, there is a possible overlap in phenology between P. remota, which disperses its pollen between March and April, and P. edulis, which disperses pollen in a short period in the spring (Lanner 1970; Farjon and Styles 1997). In the cladistic analysis of morphological characters by Malusa (1992), P. remota, which has three needles per fascicle, formed a clade with the single-needle pinyon pines (P. californiarum, P. fallax (Businský) Little, and P. monophylla). The morphological character P. remota shares with the single-leaf pinyons is an elevated number of resin canals. An increase in the number of resin canals may have arisen in the species not by independent evolution but by hybridization, likely between P. remota and P. fallax (Malusa 1992). Plastid DNA of P. remota is unique, ruling out the possibility that it was acquired recently through hybridization (Gernandt et al. 2003).

Plastid DNA can provide important insights into phylogenetic relationships, but exclusive dependence on it as a phylogenetic marker can be misleading because plastid capture through introgression has been documented in diverse plant lineages (e.g. Delgado et al. 2007; Gernandt et al. 2018; Morales-Briones et al. 2018). Reticulation through hybridization and incomplete lineage sorting (ILS) also play important roles in gene tree discordance (Rieseberg and Soltis 1991; Maddison 1997) by producing very similar patterns of shared genetic diversity and obscuring phylogenetic relationships, which greatly limits our understanding of the processes of diversification in many lineages (Maddison and Knowles 2006). On one hand, weak premating barriers to gene flow can permit the introduction of alleles from other species (Mirov 1967). On the other, ILS can reduce the differentiation between species. Therefore, conifer species with long generation times (Petit and Hampe 2006) often share genetic variation (e.g. Wang et al. 2011; DeGiorgio et al. 2014). New DNA sequencing technologies and the development of coalescent-based frameworks allow us to better study hybridization and ILS. A promising approach is the use of low-copy nuclear genes to increase phylogenetic resolution and explore the problems of discordance between gene trees and species trees. For example, Syring et al. (2005) found that phylogenetic inference based on low-copy nuclear genes in Pinus subgenus Strobus increased resolution and robustness in most of the subsectional clades recovered in their study. Gernandt et al. (2018) demonstrated the utility of low-copy nuclear genes to explore ILS and reticulation at the species level in Pinus subsection Australes. For this reason, our aims were to infer the phylogenetic relationships for species of Pinus subsection Cembroides and explore the relative importance of incomplete lineage sorting and reticulation as causes of phylogenetic discordance. We used targeted sequence capture (Gnirke et al. 2009), also known as Hyb-Seq (Weitemier et al. 2014), to characterize nuclear DNA sequences from multiple individuals per species and perform concatenated and multispecies coalescent analyses. This study represents the most complete taxonomic sampling to date for Pinus subsection Cembroides and its close relatives.

Materials and Methods

Taxon Sampling—We sampled 60 individuals, with most taxa in subsections Cembroides and Nelsoniae represented by multiple populations (Appendix 1). Vouchers were deposited in the Herbario Nacional de México (MEXU), Instituto de Biología, Universidad Nacional Autónoma of México and the Oregon State University Herbarium (OSC). The individuals represent all taxa of Pinus subsection Cembroides recognized by Gernandt et al. (2005). For outgroups we included three individuals representing two of three species of subsection Balfourianae and three of P. nelsonii, the only member of subsection Nelsoniae (Fig. 1). These two subsections were recovered as the sister group of subsection Cembroides in previous phylogenetic analyses and together with subsection Cembroides are classified in section Parrya (Gernandt et al. 2005). We also included representative species from section Quinquefoliae: subsections Gerardianae Loudon (2), Krempfianae Little and Critchfield (1), and Strobus as more distant outgroups (Appendix 1).

DNA Extraction and Quantification—For extraction of genomic DNA, we followed the modified CTAB method of Doyle and Doyle (1987) for diploid leaf tissue and used a Wizard genomic DNA purification kit (Promega, Madison, Wisconsin) for haploid seed megagametophyte tissue. We used a Nanodrop spectrophotometer 2000/2000c (Thermo Scientific, Waltham, Massachusetts) to measure the absorbance maxima ratio. Samples with 800 ng or more of DNA and an A260/A280 between 2.0 and 2.2 were selected for sequencing. We used a Qubit fluorometer v. 3.0 and dsDNA HS assay kit (Life Technologies, Carlsbad, California) to measure DNA concentration.

Probe Design—Details on probe design, library preparation, sequencing, and gene assembly were described by Gernandt et al. (2018). Briefly, a total of 1045 putative single copy nuclear genes were screened for probe design from Pinus species (P. taeda L., P. pinaster Aiton, and P. sylvestris L.; Neves et al. 2013; Willyard et al. 2007). The exon sequences were submitted to Arbor Biosciences (Ann Arbor, Michigan, USA), and 120 bp RNA bait sequences were used to perform a BLAST search on the P. taeda draft genome v. 1.0 (Neale et al. 2014; Wegrzyn et al. 2014).

Illumina Library Preparation and Target Enrichment—Genomic libraries were prepared with between 100 and 500 ng of DNA per sample. The DNA was fragmented into ca. 250 bp with a bioruptor and barcode adapters were ligated for sequencing on the Illumina using a TruSeq library prep kit (Illumina, San Diego, California). Libraries were pooled into 24 samples in equimolar ratios and enrichment was carried out with MYbaits biotinylated RNA baits following the manufacturer's protocol v. 2.3.1 (Arbor Biosciences, Ann Arbor, Michigan). The samples were combined in equal concentrations (48×) and sequenced using an Illumina Hi-Seq 2500 with the 100 bp module with paired reads.

Data Selection—We used two different data sets as input for phylogenetic analyses. The principal data set consisted of genes assembled with HybPiper v. 1.2 (Johnson et al. 2016). Data are available in the Dryad Digital Repository (Montes et al. 2019). A second data set consisted of single nucleotide polymorphisms (SNPs) identified with SAMtools (see below). For both, Illumina reads were filtered in Trimmomatic v. 0.36 (Bolger et al. 2014). For nuclear genes, a total of 60 paired R1 and R2 files in fastq format were filtered in Trimmomatic, removing bases at read ends with qualities < Q20 using a 4 bp sliding window, and removing reads with a length < 30 bp following trimming (Weitemier et al. 2014). Only reads with both pairs surviving were assembled into individual alignments with HybPiper (Johnson et al. 2016). HybPiper used BWA v. 0.7.1 (Li and Durbin 2009) to align reads to the reference nuclear gene sequences. SAMtools v. 0.1.19 (Li et al. 2009) was used to sort the reads into separate directories for each gene. The pipeline subsequently used SPAdes v. 3.10.1 (Bankevich et al. 2012) for de novo assembly of each gene individually using the retrieved reads. A total of 969 genes from 996 targets were assembled successfully for at least one sample.

The resulting gene assemblies were imported into Geneious v. R9 (Kearse et al. 2012). Individual nuclear gene alignments were performed with MAFFT v. 7.0 (Katoh et al. 2002). We used the following criteria proposed by Gernandt et al. (2018) to filter the multiple sequence alignments: 1) missing one or more samples, 2) fewer than 50% of sites, 3) pairwise similarity less than 93%, and 4) putative paralogs. The first three criteria were applied in Geneious, whereas the paralogs were detected when assembling the only two haploid references, P. cembroides USA and P. bungeana CA, with HybPiper (Johnson et al. 2016). Paralogs script identified contigs with lengths ≥ 85% of the reference sequence, indicating multiple long-length matches. After filters, we excluded 665 of the 969 genes assembled. The remaining 304 genes were carried forward for phylogenetic analysis.

Phylogenetic Analysis Using Low-Copy Number Nuclear Genes—We analyzed a concatenated alignment of 304 nuclear genes for 60 individuals with maximum likelihood in RAxML-HPC v. 8.2.10 (Stamatakis 2014). We performed 1000 heuristic searches for the best tree applying a general time reversible model with the gamma parameter (GTR + G) separately to each gene. Bootstrapping was performed with 1000 replicates as part of the heuristic search with RAxML. The best tree was imported into FigTree v. 1.4.0 for further editing (Rambaut 2012). Samples with unstable topological positions (“rogues”) were identified using the maximum dropset size in RogueNaRok (Aberer et al. 2013), available as an online webserver ( http://rnr.h-its.org/). These individuals were excluded from the concatenated alignment and multispecies coalescent analyses.

Phylogenetic Analysis Using Coalescent-Based Methods—Coalescent analyses that accommodate ILS and hybridization were performed on the low-copy nuclear genes. Because estimating reticulation is computationally demanding, we used a nuclear gene tree as a guide for choosing a sub-sample of individuals representative of Pinus subsection Cembroides. The species tree was estimated with minimizing deep coalescences (MDC), a parsimony method (Maddison 1997; Than and Nakhleh 2009) and networks were estimated by maximum parsimony from gene trees. To study reticulation scenarios, we used PhyloNet v. 3.6.1 (Than et al. 2008) to analyze 22 sequences representing 22 taxa, of which 15 belong to subsection Cembroides. The remaining seven taxa were from the outgroup (one individual per species). We used the Perl script BeforePhylo.pl v. 0.9.0 ( https://github.com/qiyunzhu/BeforePhylo.git) to produce individual gene alignment files with the reduced representation of samples. The maximum likelihood tree for each of the 304 genes was inferred with RAxML using the GTR + G model and bootstrapping with the autoMRE option. We sequentially permitted up to three reticulation events under the MDC criterion (InferNetwork_MP). For each reticulation setting, we performed 10 independent searches of 5 to 20 (-x 5 to -x 20). The network with the lowest number of extra lineages was selected and displayed graphically with Dendroscope v. 3.0 (Huson and Scornavacca 2012). Based on the results from the 22-terminal analysis, we evaluated the effect of constraining P. remota and P. quadrifolia as hybrids. We included the taxa proposed to have been involved in past reticulation events with these species. For P. remota, these include P. cembroides, P. edulis, and P. fallax (Malusa 1992; Little 1966). For P. quadrifolia we included P. californiarum (Lanner 1974b).

Species and lineage tree inference was also performed with SVDquartets (Chifman and Kubatko 2014) in PAUP* v. 4.0a150 (Swofford 2002). SVDquartets is a robust method for multilocus data that is designed to build unrooted quartets that accommodate ILS but not reticulation and evaluates the correspondence of the nodes (see Chifman and Kubatko 2014). For this method we used a NEXUS input file that included a data set block specifying the 304 gene partitions and a taxon set assigning each individual to its corresponding species. The file included 58 sequences representing the 15 subsection Cembroides taxa and ten individuals from the outgroup, the same as in the concatenated analysis. A maximum of 1,000,000 randomly chosen quartets were analyzed and branch support was estimated with 1000 bootstrap replicates (bs). The trees were displayed graphically in FigTree v. 1.4.0 (Rambaut 2012).

In addition, we concatenated these 304 low-copy number nuclear genes with SNPs (see below) for 51 samples to infer a species tree with SVDquartets. The species trees resulting from the nuclear gene alignment and the concatenated SNPs analyses were displayed graphically as a tanglegram with Dendroscope (Huson and Scornavacca 2012).

The species tree was also estimated with ASTRAL-III v. 5.6.3 (Mirarab and Warnow 2015; Zhang et al. 2017). We used 58 individuals representing 22 taxa, removing two samples, one because it had an unstable position in the RAxML analysis of the concatenated alignment and the other because of the amount of shared ancestry between two species observed with SNPs. Individual trees from the 304 nuclear genes were estimated with RAxML using the GTR+G model and bootstrapping with the autoMRE option. The best tree from each analysis was concatenated and used as the input file for ASTRAL-III. Branch lengths in coalescent units and local posterior probabilities were estimated for the species tree (Sayyari and Mirarab 2016). We performed character state reconstruction by mapping of the number of leaves per fascicle as an ordered multistate character on the species tree inferred with ASTRAL-III to evaluate the origin of single-needle pinyons based on the likelihood ancestral state approach in Mesquite v. 3.6 (Maddison and Maddison 2018).

Sequence Read Alignment and SNP Calling—Single nucleotide polymorphisms constitute a valuable source of genetic variation to study evolutionary relationships in non-model organisms (Leaché and Oaks 2017). Although only 304 assembled genes were selected to carry out phylogenetic analyses, we observed that for other target genes some regions were assembled and could be used to identify genetic markers in spite of our gene filtering criteria. For this reason, we performed reference-guided SNP calling in an effort to increase the number of sites for phylogenetic analysis. This method provided an alternative way of interpreting the Illumina data compared to assembling with HybPiper and applying our ad hoc filters. It has the potential to identify more informative sites and reduce bias introduced by HybPiper, which only returns one allele per gene. Also, it might be more reliable for removing paralogs or improving homology across samples. Demultiplexed sequence reads from each of the 52 samples were evaluated using FastQC and MultiQC and filtered with Trimmomatic v. 0.36 to discard low quality and adapter sequences (Bolger et al. 2014; Ewels et al. 2016). The minimum length of reads to be kept was set to 50 bp. Paired cleaned sequences were mapped against the Pinus taeda genome v. 2.0 (Neale et al. 2014; Wegrzyn et al. 2014) using BWA-MEM with default parameters. The SAM files were converted to BAM files with SAMtools v. 0.1.19 (Li et al. 2009). Uniquely mapped and sorted reads were obtained for each alignment file using the view and sort routines from SAMtools. Potential PCR duplicates were discarded using the Picard tool MarkDuplicates ( http://broadinstitute.github.io/picard). SAMtools mpi-leup and BCFtools (with the options for biallelic variants SNPs/indels, no-BAQ, minimum mapping quality of 20, and minimum base quality of 25) were used for variant calling. SNPs were further filtered with VCFtools v. 0.1.13 (Danecek et al. 2011) for phylogenetic and admixture analyses.

SNP-Based Phylogenetic and Admixture Analyses—The SNPs were filtered if genotypes were not called across all samples (100%), minimum mean depth was below 5, minimum quality score was below 30, or minor allele frequency was lower than 0.05. The variant calling format (VCF) file was converted into a tab-delimited text file with VCFtools and subsequently SNPs were concatenated into a FASTA file including heterozygous sites. A multiple sequence alignment was generated with MAFFT v. 7.0 (Katoh et al. 2002) and used as input to infer a maximum likelihood tree in RAxML-HPC v. 8.0.26 (Stamatakis 2014). The maximum likelihood and bootstrap searches were carried out with the GTR+G model and a total of 1000 bootstrap replicates were computed.

For the admixture analysis, SNPs were filtered if genotypes called were below 80% across all samples, minimum mean depth was below 5, and minimum quality score was below 30. The VCF file was converted to an ordinary PLINK file (.ped) using PLINK v. 1.9 (Purcell et al. 2007). Maximum likelihood estimation of individual ancestries (population structure analysis) was carried out using ADMIXTURE v. 1.3.0 (Alexander et al. 2009). ADMIXTURE's cross-validation procedure (5-fold) was run for ancestral populations values (K) from 2–5. The Q-matrices generated for the different K values (1–4) were then clustered and plotted using CLUMPAK v. 1.1 (Kopelman et al. 2015).

Results

Filtering and Editing Sequences—Statistics for the 60 samples assembled in HybPiper are summarized in Appendix 2. Fifty samples corresponded to subsection Cembroides and 10 to the outgroups. Mean coverage of the 969 genes was 226.42×. We excluded 68.6% of genes as a result of the data filtering steps (Table 3), resulting in a final tally of 304 included gene alignments.

Table 3.

Results of the data after filtering. ^aOne or more samples were not assembled to the reference sequence (996 genes). ^bFewer than 50% of sites identical. ^cPairwise identity less than 93%. ^dThe HybPiper script identified contigs with lengths ≥ 85% of the reference sequence, indicating multiple long-length matches (see Johnson et al. 2016).

Phylogenetic Analyses Using Nuclear Genes—The alignment of 304 concatenated nuclear genes was 222,129 bp in length and included 16,503 parsimony informative sites and 13,639 variable but parsimony uninformative sites. The RAxML analysis recovered Pinus subsection Cembroides as monophyletic (100% bs). Rooting with section Quinquefoliae resulted in recovering subsections Nelsoniae and Balfourianae sister to one another and in turn sister to subsection Cembroides (Fig. 2). This section Parrya clade of exclusively North American taxa was recovered with high bootstrap support (100% bs). In the outgroup, the subsections Krempfianae and Gerardianae were united (100% bs). In the subsection Cembroides clade, eight of 13 taxa represented by multiple individuals were recovered as exclusive lineages (P. culminicola, P. johannis, P. fallax, P. maximartinezii, P. monophylla, P. pinceana, P. remota, and P. rzedowskii; Fig. 2). The bootstrap values for most branches in the phylogeny were high. The three large-cone pinyon pines with a restricted geographical distribution in Mexico formed a well-supported clade, in which P. maximartinezii and P. pinceana formed a monophyletic group (100% bs) sister to P. rzedowskii (100% bs). The large-cone clade was sister to a clade of the remaining (small-cone) species of subsection Cembroides. In the small-cone clade, P. johannis and P. discolor were sister to P. culminicola (94% bs). Five individuals of P. remota were sister to the clade of P. johannis, P. discolor, and P. culminicola (92% bs). Pinus cembroides subsp. cembroides, P. cembroides subsp. orizabensis, and P. lagunae were paraphyletic to the clade of P. discolor, P. johannis, P. culminicola, and P. remota, with individuals of P. cembroides subsp. cembroides from Chihuahua and San Luis Potosí forming a clade with two individuals of P. lagunae (Baja California Sur), and P. cembroides subsp. orizabensis (Puebla) sister to an individual of P. cembroides subsp. cembroides from Hidalgo. The two individuals of P. fallax formed a group with one sample of P. edulis (93% bs), all from Utah (USA). Pinus monophylla, P. californiarum, and P. quadrifolia were recovered as a clade with high bootstrap support (98% bs). Pinus monophylla was recovered as monophyletic and sister to P. quadrifolia, and both as sister to a single sample of P. californiarum from Baja California.

Only one sample, Pinus monophylla UT1 (DSG478), was identified as having an unstable position with RogueNaRok (dropset size 3.0). We did not include it in the subsequent coalescence analyses with MDC, SVDquartets, and ASTRAL-III, and it was removed from the concatenated analysis in RAxML. Also, one sample of P. californiarum (CA2; DSG403) was removed from the coalescence and concatenated nuclear genes analyses to avoid the probability of bias in the evolu-tiionary relationships due to the proportions of shared ancestry observed for this sample with SNPs (see Fig. 3A).

SNP-BasedPhylogeneticandAdmixtureAnalyses—Mapping the Illumina reads to the genome of P. taeda identified 26,499 SNPs genotyped in all 52 pines, including 7262 parsimony informative sites and 735 variable but parsimony uninformative sites. The resulting maximum likelihood tree inferred with RAxML (Stamatakis 2014; Fig. 3B) had a similar backbone topology as that inferred from the concatenated nuclear genes (Fig. 2), with Pinus subsection Cembroides monophyletic (Fig. 3B) and divided into two main clades comprising the large-cone pinyon pines (100% bs) and the small-cone pinyon pines (100% bs). In the small-cone clade, P. californiarum, P. monophylla, P. quadrifolia, P. fallax, and P. edulis formed a well-supported monophyletic group (90% bs). Pinus culminicola, P. johannis, P. discolor, P. lagunae, P. cembroides subsp. cembroides, P. cembroides subsp. orizabensis, and P. remota formed another well-supported clade (100% bs). In contrast to the analysis of the HybPiper exon assembly (Fig. 2), the maximum likelihood tree inferred from SNPs grouped the two subspecies of P. cembroides and P. lagunae together in an exclusive lineage (Fig. 3B). Pinus californiarum individuals were sister to P. monophylla + P. quadrifolia, but the two samples of P. californiarum were paraphyletic in the SNP-based phylogeny, suggesting significant genetic variation within this species (Fig. 3B)

Fig. 2.

Phylogeny of Pinus subsection Cembroides. Maximum likelihood tree inferred from the concatenated alignment (222,129 bp; 58 terminals and 304 genes). Bootstrap values > 50% are shown above the branches. The taxonomic subsections are represented by colors in some trees in this study (blue = Cembroides; green = Nelsoniae; navy blue = Balfourianae; orange = Gerardianae; violet = Krempfianae; black = Strobus). The sample names indicate the taxon and the state of collection. Locality codes are provided in Table 1.

Fig. 3.

Admixture analysis of a pinyon pine clade. A. Analyses were performed for K values ranging from 2 to 4 with a matrix containing 26,499 single nucleotide polymorphisms. Different colors represent different clusters. The combination of different colors in a bar indicates the degree of admixture. Samples in the admixture model are in the same topological order as in the maximum likelihood phylogeny which is shown at right for comparison. B. Maximum likelihood tree inferred from single nucleotide polymorphisms. Bootstrap values > 50% are shown above the branches.

The interesting phylogenetic relationships for the clade comprising P. californiarum, P. monophylla, P. quadrifolia, P. fallax, and P. edulis, together with previous evidence of hybridization and introgression between some of these species, led us to explore the genetic structure of the individuals from this clade using ADMIXTURE. In total, 67,938 SNPs corresponding to only these five species were used. The cross-validation errors for ancestral structured populations values increased as the number of populations increased. Therefore, although the lowest cross-validation error was for K = 2, we used a cluster range of 2–4 to explore the dynamics of classifications in the different number of populations. Interestingly, for the lowest K all the samples were mainly defined by one subpopulation with the exception of the sample corresponding to P. californiarum CA2, which showed the strongest signs of admixture (Fig. 3A). For K = 3, samples from P. quadrifolia defined a third substructure and admixture patterns for sample CA2 was maintained. Pinus edulis and P. fallax were grouped differentially only at the highest K value. Pinus californiarum CA2 and P. californiarum BC1 consistently showed evidence of admixture across all three values of K. Pinus monophylla individuals were mainly assigned to the same population regardless of which K value was used (Fig. 3A). Also, P. edulis and P. fallax showed consistency regarding their classification in the same population for the first two K values.

Phylogenetic Results for Coalescent-Based Methods—HybPiper exon data were used to infer a species tree under the MDC criterion with Phylonet on a subset of subsection Cembroides individuals (22). Pinus subsection Cembroides was recovered as monophyletic (8084 lineages; Fig. 4). The analyses sequentially permitting up to three reticulation events (Fig. 4B–D) always identified gene flow within subsection Cembroides; none was identified among the outgroups. Allowing for a single reticulation resulted in a reduction from 8084 to 7619 lineages and involved introgression from P. monophylla into P. edulis (Fig. 4B). Allowing for two reticulations resulted in a reduction to 7453 lineages (Fig. 4C). The first reticulation involved introgression from P. fallax into P. cembroides subsp. cembroides (Fig. 4C: H1) and the second involved introgression from P. edulis into P. lagunae (Fig. 4C: H2). Allowing for three reticulations resulted in a reduction to 7211 lineages (Fig. 4D). The first reticulation (Fig. 4D: H1) involved introgression from P. culminicola into P. quadrifolia, the second involved introgression from a possible extinct taxon into P. lagunae (Fig. 4D: H2), and the third involved introgression from a possible extinct taxon into P. cembroides subsp. cembroides (Fig. 4D: H3).

Specifying P. remota as a hybrid under MDC resulted in a reduction from 8084 lineages inferred in the species tree (Fig. 4A) to 7683 in the network (Fig. 5A). Pinus remota was sister to P. johannis and gene flow was inferred from P. fallax (28% inheritance probability). Specifying P. quadrifolia as a hybrid resulted in a reduction to 7627 lineages (Fig. 5B). In this network P. quadrifolia was sister to P. monophylla (54% inheritance probability) and gene flow was inferred from P. lagunae (46% inheritance probability) (Fig. 5B).

Pinus subsection Cembroides was monophyletic and relationships among the outgroups were well supported in the SVDquartets and ASTRAL-III trees (Figs. 6B–8). The Mexican pinyon pine P. nelsonii (subsection Nelsoniae) and P. aristata and P. longaeva were sister to subsection Cembroides in the SVDquartets lineage and species trees (Figs. 6–7) but was sister to subsections Cembroides + Balfourianae in the ASTRAL-III tree (Fig. 8). Pinus maximartinezii and P. pinceana were sister and in turn sister to P. rzedowskii in all trees (Figs. 2–8). In the lineage tree individuals of P. monophylla formed a group with P. quadrifolia and P. californiarum but individuals of the same species were not recovered as exclusive lineages. Individuals of P. edulis also were not recovered as exclusive lineages (Fig. 6B). Pinus lagunae and the taxonomic varieties of P. cembroides were not recovered as exclusive lineages (Fig. 6), whereas P. remota, P. culminicola, and P. johannis were each recovered as exclusive lineages (> 95% bs).

In both the SVDquartets and ASTRAL-III analyses, P. remota was sister to the P. culminicola, P. johannis, and P. discolor clade (Figs. 7–8). The relationships among P. discolor, P. culminicola, and P. johannis differed between ASTRAL-III and SVDquartets (Figs. 6–8). Whereas P. culminicola was sister to P. johannis and P. discolor in the SVDquartets species tree, in the ASTRAL-III tree P. discolor was sister to P. culminicola and P. johannis (Fig. 8). The position of P. fallax and P. edulis also differed in the two coalescence analyses. Moreover, the position of P. edulis as sister to both varieties of P. cembroides, P. culminicola, P. discolor, P. johannis, P. lagunae, and P. remota was not well supported in the SVDquartets species tree (65% bs) indicating a weak node (Fig. 7). Branch lengths in the tree inferred with ASTRAL-III were short within subsection Cembroides except for the branch subtending the clade with P. maximartinezii, P. pinceana, and P. rzedowskii.

The likelihood results for character state reconstruction on the ASTRAL-III species tree supported a common origin of reduction to single needles in P. californiarum, P. fallax, and P. monophylla followed by a single increase in the number of needles per fascicle in P. quadrifolia (-log 33.50).

Comparison Between SVDquartets Analyses from SNPs and Nuclear Genes—The concatenated nuclear gene alignment was 222,129 bp in length (16,503 parsimony-informative sites) compared to 26,499 characters for the SNP data set (7262 parsimony-informative sites). The analyses performed in SVDquartets recovered species trees that were topologically similar. Trees based on both SNPs and concatenated nuclear genes (Fig. 9) recovered subsection Cembroides as monophyletic with high branch support and the large- and small-cone clades were sister groups in both analyses. The relationships within subsection Cembroides differed for several small-cone species. The bootstrap values for the analyses based on concatenated genes and SNPs were similar, but in some positions the support was higher using concatenated genes (Fig. 9B). In the concatenated gene analysis, P. remota was sister to the P. culminicola + P. discolor + P. johannis clade (96% bs), whereas with SNPs, P. cembroides subsp. orizabensis was sister to the P. culminicola + P. discolor + P. johannis clade (61% bs) with SNPs. The clade of P. cembroides subsp. cembroides, P. cembroides subsp. orizabensis, and P. lagunae was recovered as monophyletic with the concatenated gene alignment (72% bs), but with SNPs it was paraphyletic due to the placement of P. cembroides subsp. orizabensis as sister to the clade of P. discolor, P. johannis, P. culminicola, and P. remota (61% bs). In the SNPs tree, P. fallax was sister to the single-needle pinyons together with P. quadrifolia (63% bs) whereas in the concatenated gene tree it was sister to all other small-cone species (100% bs). In both analyses of SVDquartets, P. monophylla was sister to P. quadrifolia rather than to P. californiarum. With SNPs, P. nelsonii + P. aristata were sister to subsection Cembroides (50% bs), whereas with the nuclear exon alignment only P. nelsonii was recovered as sister to subsection Cembroides (61% bs).

Fig. 4.

Results of minimizing deep coalescences (MDC) analyses with and without reticulation. Inheritance probabilities for the minor edge are indicated for the networks (B, C, and D), with the reticulations represented by lines. A. Best MDC tree with no reticulation. B. Best MDC network permitting one reticulation event. C. Best MDC network permitting two reticulation events. D. Best MDC network permitting three reticulation events.

Fig. 5.

Tests constraining P. remota and P. quadrifolia to be hybrids under the minimizing deep coalescences criterion (MDC). A. Best MDC network constraining P. remota as a hybrid. B. Best MDC network constraining P. quadrifolia as a hybrid. Reticulation events are represented by lines and inheritance probabilities are indicated at reticulations.

Discussion

Phylogeny of Pinus Subsection Cembroides—The phylogenies inferred with low-copy nuclear genes for Pinus subsection Cembroides recovered two main lineages (Figs. 2, 8). Nonetheless, some interrelationships vary among analysis and with previous studies. In studies with the nrDNA ITS region (Gernandt et al. 2001) and cpDNA (Gernandt et al. 2005; Parks et al. 2012; Ortiz-Medrano et al. 2016), P. rzedowskii is sister to the remaining species of Pinus subsection Cembroides. Here, P. rzedowskii is sister to P. pinceana and P. maximartinezii (100% bs). The same relationship of these large-cone pinyon pines was recovered with plastid DNA sequences by Gernandt et al. (2007) but did not receive high support (72% bs). In both the concatenated and coalescence-based analyses, the P. rzedowskii + P. pinceana + P. maximartinezii clade is always sister to the rest of the pinyon pines with strong support (100% bs). If P. rzedowskii forms a clade with P. maximartinezii + P. pinceana, instead of being the sister group of all the other species in subsection Cembroides, this implies either multiple gains of enlarged, functionally wingless seeds in subsection Cembroides, or a reversion to small winged seeds in P. rzedowskii (Gernandt et al. 2007). Maximum likelihood analyses of plastid DNA support the sister relationships between P. discolor + P. johannis and P. culminicola (Gernandt et al. 2007; Parks et al. 2012). The results with low-copy nuclear genes support that relationship with the concatenated alignment in RAxML and the SVDquartets species tree (Figs. 2, 7). In contrast, in the ASTRAL-III and SNPs trees (Figs. 3A, 8), P. discolor is recovered as sister to P. culminicola + P. johannis; this relationship was also supported with plastomes (Parks et al. 2012). Although Farjon and Styles (1997) treated P. discolor and P. johannis as a single variety of P. cembroides (as P. cembroides var. bicolor Little), the stomata of P. discolor, P. johannis, and P. culminicola are limited to the adaxial surfaces and the seed megagametophyte is white rather than pink. Our results with low-copy nuclear genes from the concatenated alignment, SVDquartets, and ASTRAL-III, agree with analyses of plastid DNA (Gernandt et al. 2007) in recovering P. johannis + P. discolor + P. culminicola as the sister group of P. remota, rather than grouping these species with P. cembroides. This relationship was not recovered with cladistic analyses of morphological characters (Malusa 1992). The coalescence analyses at the species level strongly support the clade of pinyon pines that are distributed in the Sierra Madre Oriental (SMO), P. discolor (mainly distributed in the Sierra Madre Occidental and the Sky Islands of the United States, but also present in the southern part of the Sierra Madre Oriental), P. johannis, and P. culminicola (100% bs in SVDquartets and ASTRAL-III; Figs. 7–8). The monophyly of the three taxa has been attributed to the evolution from an ancestor that was resistant to the calcareous soils that predominate in the SMO (Malusa 1992). Therefore, limestone soils tolerance may be a synapomorphy for P. discolor, P. johannis, and P. culminicola.

Infraspecific taxonomies of P. cembroides and their relationship to P. lagunae have disagreed (Zavarin and Snajberk 1985; Passini 1987; Farjon and Styles 1997; Farjon and Filer 2013). Phylogenetic results based on three plastid DNA regions recovered P. cembroides subsp. cembroides, P. cembroides subsp. orizabensis, and P. lagunae as a single exclusive lineage (Gernandt et al. 2003). Our results from nuclear gene assemblies (but not from SNPs; Fig. 9A) coincide with Gernandt et al. (2003) in recovering both varieties of P. cembroides and P. lagunae as an exclusive lineage (Figs. 7–9). Previously, Whang et al. (2001) reported differences of the leaf internal cuticle among P. cembroides subsp. cembroides, P. cembroides subsp. orizabensis, and P. lagunae. For instance, P. cembroides subsp. orizabensis differs from P. cembroides subsp. cembroides and P. lagunae by the width of the epidermal cell apex (thick in P. cembroides subsp. orizabensis), continuity of cell walls, stomatal apparatus shape, and cuticular flange-guard cell. Pinus cembroides subsp. cembroides and P. lagunae share more characters of the internal cuticle with each other than with P. cembroides subsp. orizabensis.

Fig. 6.

Tanglegram of RAxML tree and SVDquartets lineages tree. A. Maximum likelihood tree based on concatenated alignment with 304 nuclear genes. B. Coalescent-based tree based on 1,000,000 quartets. Trees were estimated with a matrix of 222,129 bp and 58 terminals. Bootstrap values > 50% are shown above the branches. The colors of the branches follow Fig. 2.

Zavarin and Snajberk (1985) found that the populations of P. cembroides subsp. orizabensis from southern Puebla and northeastern Veracruz differ from Pinus cembroides subsp. cembroides and P. lagunae in their chemical composition of monoterpenes but are very similar morphologically. They suggested that the divergence of southern populations of P. cembroides (P. cembroides subsp. orizabensis) as an isolated taxon most likely resulted from climatic and geographic isolation (middle Miocene). They also suggested that the isolation of P. lagunae from Baja California was related to movement of the coastal region from California during the Miocene and the formation of the Sierra Madre Occidental, resulting in the separation of this population from the rest of the continental populations by the breach of the Gulf of California (Zavarin and Snajberk 1985). However, this vicariance event is probably too old to explain divergence between P. lagunae and P. cembroides. The entire small-cone pinyon pine clade was estimated to have diversified during the Miocene, ∼11 MYA (Gernandt et al. 2008; Saladin et al. 2017). More studies are needed to determine whether these three taxa represent independent evolutionary entities or the same species.

Fig. 7.

Species tree inferred with SVDquartets. Coalescent-based tree estimated from 1,000,000 quartets (222,129 bp, 58 taxa and 304 genes). Total weight of incompatible quartets = 12.89%, and total weight of compatible quartets = 87.11%. Bootstrap values > 50% are shown above the branches.

The coalescence-based analyses of concatenated genes and SNPs support the sister relationship between P. monophylla and P. quadrifolia and place both taxa as sister to P. californiarum. This relationship was not recovered with morphology, where P. quadrifolia is sister to P. johannis + P. discolor + P. culminicola, united by sharing resinous cones (Malusa 1992). In this study we observed two different topologies with respect to the origin of a single needle. In both the SVDquartets and ML analysis of the concatenated alignment, single-needle taxa occurred in independent lineages (Figs. 2, 7). In contrast, ASTRAL-III and analyses of SNPs grouped P. californiarum, P. fallax, and P. monophylla together but paraphyletic to P. quadrifolia (Figs. 8–9). We performed a character state reconstruction with the ASTRAL-III species tree and the results supported a single reduction to single needles and common origin but with one independent loss in P. quadrifolia. Cole et al. (2013) compared variation in needle number to environmental variation and found that the proportions of the number of needles in P. edulis and P. fallax depend on annual fluctuations in precipitation. In addition, it was shown that P. edulis and P. fallax occur in an area with monsoon precipitation extremes, whereas P. monophylla and P. californiarum occur in areas with high levels of winter precipitation (Cole et al. 2008).

Fig. 8.

Species tree inferred with ASTRAL-III. The tree was estimated under a species coalescent model using 304 nuclear gene trees with 58 individuals as input. Final quartet score = 73,455,073 and final normalized quartet score = 0.58 (very high incomplete lineage sorting). The bootstrap support values are provided above branches. The branch length represents coalescent units. The colors of the branches follow Fig. 2.

In morphology-based views of phylogeny, P. monophylla, P. californiarum, and P. fallax are recovered together by sharing resinous cones, single needles (predominantly), and thinner seed coats (Malusa 1992). With nuclear genes, Pinus monophylla, P. californiarum, and P. fallax are not recovered as monophyletic by the ASTRAL-III and SVDquartets analyses (Figs. 7–8). In fact, in our results P. fallax and P. californiarum are separate from Pinus monophylla, suggesting that P. fallax and P. californiarum are not taxonomic varieties of P. monophylla as proposed by Silba (1990). Taxonomic uncertainty between P. monophylla and P. edulis can be attributed to the existence of trees with both one and two leaves per fascicle (Tausch and West 1987). However, P. edulis is not sister to P. monophylla. Furthermore, environmental studies indicate that P. edulis is more similar to P. fallax, and P. monophylla is more similar to P. californiarum. Besides, P. fallax occurs in an area with moderate summer rains, similar to P. edulis (Malusa 1992). Our analyses support the separation of P. californiarum from P. monophylla. Bailey (1987) segregated P. californiarum from P. monophylla based principally on the length and amount of curl-back of fascicle sheaths, the number of leaf resin canals, and the number of rows of foliar stomata. The ASTRAL-III analysis recovered a clade with the species from the southwestern United States, P. monophylla, P. californiarum, P. quadrifolia, P. edulis, and P. fallax. The pinyon pines in this region are mainly allopatric or parapatric in distribution (Malusa 1992), but P. californiarum and P. quadrifolia co-occur in California and Baja California. Pinus monophylla occurs in California, extending east (and north) into the Great Basin in the states of Utah, Colorado, Arizona, and Idaho (Critchfield and Little 1966; Farjon 2005; Cole et al. 2008), and populations of P. fallax and P. edulis co-occur in Arizona and New Mexico (P. edulis reaches eastern Nevada and southeastern California). If P. edulis and P. fallax (adapted to early summer or periodic drought) are paraphyletic to P. monophylla (adapted to summer-autumn drought), P. californiarum, and P. quadrifolia (Figs. 2, 7), this relationship may be explained by vicariance or may have evolved in response to summer drought (Cole et al. 2008). In fact, these taxa are distributed widely in geographical regions with distinct precipitation regimes (Cole et al. 2008). Particularly, P. monophylla occurs in regions with different precipitation from P. fallax and P. edulis, which are characterized by high monsoon precipitation (Cole et al. 2008). The ecological similarities of P. fallax to P. edulis (Bailey 1987) rather than to P. monophylla, and its genetic distinctness from P. monophylla (Fig. 8), support our decision to treat P. fallax as a separate species from P. monophylla.

Fig. 9.

Comparison between nuclear genes-based and single nucleotide polymorphisms-based trees inferred in SVDquartets. A. Single nucleotide polymorphisms tree based on 26,499 characters with a total weight of incompatible quartets = 5.35%, and total weight of compatible quartets = 94.65%. B. Coalescent-based tree based on 304 nuclear gene trees with a total weight of incompatible quartets = 14.05%, and total weight of compatible quartets = 85.95%. The trees were estimated from 1,000,000 quartets in 51 terminals. Bootstrap values > 50% are shown above the branches.

Hybridization and Introgression—Both natural and artificial hybridization have been well-documented in conifers (Saylor and Smith 1966; Lanner 1974a, 1974b; Delgado et al. 2007; Wachowiak et al. 2011; Zhou et al. 2017). Overlapping phenology and weak reproductive barriers can influence the direction of pollen-mediated gene flow in natural populations (Hamilton et al. 2013). In conifers, plastid DNA is paternally inherited, and higher gene flow for plastid DNA has been attributed to the high migration capacity of wind-dispersed pollen (Rieseberg and Soltis 1991; Petit et al. 2005).

In pines, plastid introgression has been observed at a low or moderate frequency in sympatric or parapatric populations (Dounavi et al. 2001; Delgado et al. 2007; Zhou et al. 2017). In closely-related Pinus species, shared genetic variation can be the result of introgression following secondary contact, with genetic differentiation in parapatric populations lower than in allopatric populations (Zhou et al. 2017). Using the minimizing deep coalescence criterion, we explored reticulation in subsection Cembroides. No gene flow was detected between subsection Cembroides and other closely related lineages (Figs. 4–5). This result supports the observation by Mirov (1967) that species of subsection Cembroides do not form hybrids with species from other pine subsections. In fact, Pinus subsection Cembroides may have diverged from other subsections (particularly Balfourianae) relatively early (Axelrod 1986).

The majority of taxa involved in the reticulation events detected here have geographic distributions that are somewhat close to one another (P. californiarum, P. edulis, P. fallax, P. monophylla, P. quadrifolia, and P. remota; only P. lagunae and P. cembroides subsp. orizabensis are geographically isolated). This coincides with other studies where gene flow has been reported (Edwards-Burke et al. 1997; Delgado et al. 2007; Zhou et al. 2017). Our results coincide with some hypotheses proposed by Lanner (1974a). We detected gene flow in P. edulis in only one reticulation scenario (Fig. 4B). This suggests that P. edulis is introgressed with P. monophylla. Some populations of P. edulis and P. monophylla occur in sympatry in the eastern Great Basin where trees of P. edulis with both single needles and two needles per fascicle have been observed (Lanner 1974a). Pinus edulis also occurs in Arizona and New Mexico (Cole et al. 2013), and sympatric populations of P. edulis and P. monophylla have been reported in in the Mojave Desert in southeastern California (Munz and Keck 1959; Critchfield and Little 1966), western Utah, and Nevada (Farjon and Filer 2013). For this reason, it would not be unusual if P. edulis is introgressed with P. monophylla where their populations are in contact. The direction of gene flow inferred from P. monophylla to P. edulis is consistent with the prevailing winds, which are from west to east. Overlapping phenology could facilitate introgression from P. monophylla to P. edulis, since both disperse pollen in a short period in the spring (Lanner 1970; Farjon and Styles 1997). Furthermore, the distribution range of P. edulis may have been more extensive in the past (Cole et al. 2008), resulting in more widespread contact with P. monophylla.

The MDC method detected reticulation in taxa for which gene flow had not been suspected. Currently Pinus fallax and P. cembroides subsp. cembroides have widely separated distributions. One possible explanation for reticulation between the two species is that this inference is incorrect (Fig. 4C). As an alternative explanation, we suggest studying species distribution and demographic history to test whether these two taxa came into contact in the past. The potential distribution and demographic history for P. fallax suggests that it may have been more widely distributed in southwestern California, southern Nevada, throughout Arizona and extending beyond into Utah, Colorado, and New Mexico (Cole et al. 2008). Another detected reticulation event that was unexpected was between P. lagunae and P. edulis (Fig. 4C). These species also are widely separated geographically. The genetic diversity shared by these allopatric taxa seems more likely to be influenced by the retention of ancestral polymorphism. However, ancient introgression events may have been interrupted by migration of the populations of P. edulis northward in present-day USA during the Holocene (Cole et al. 2008).

We detected reticulation in P. quadrifolia and P. culminicola (Fig. 4D) although their populations are allopatric and gene flow had not been suspected. This inference may be incorrect, or long-distance pollen dispersal could have resulted in introgression. It would be interesting to study the past distribution and demographic history of P. culminicola, P. cembroides subsp. cembroides, and P. lagunae to test whether they were formerly in contact with other species.

The Phylonet analysis also detected reticulation in taxa such as P. lagunae and P. cembroides subsp. cembroides (Fig. 4D) for which gene flow had not been suspected. The origin of these reticulations implies the existence of an extinct taxon. It is also possible that this inference is incorrect or not significant. Copetti et al. (2017) explained the origin of a reticulation with the existence of an extinct or unsampled taxon in cacti, but the authors do not discuss the result.

Specifying P. quadrifolia as a hybrid under MDC (Fig. 5B), our results did not indicate that P. quadrifolia is introgressed with P. californiarum as reported by Lanner (1974b). Pinus californiarum was recovered as the sister to P. quadrifolia and P. monophylla, but no reticulation was detected between P. californiarum or P. monophylla. We did detect reticulation in P. quadrifolia and P. lagunae for which gene flow had not been suspected (Fig. 5B). Long-distance pollen dispersal could have resulted in introgression. Likewise, Pinus lagunae may be a relictual population left behind from a time when its ancestors had a range that extended northwest into what is today Baja California and southern California.

Another species reported as a possible hybrid is P. remota. Our results do not support the proposal of Little. Little (1968) suggested the possibility that in the past its populations had been in contact with those of P. edulis, permitting introgressive hybridization. Nonetheless, reticulation with P. remota was inferred from P. fallax and not from P. edulis (Fig. 5A; Little 1968). No contemporary P. fallax populations come into contact with P. remota, but past secondary contact could have resulted in introgression. Some populations of P. fallax (western New Mexico) occur in proximity to P. remota (Texas and northeast Mexico). According to the packrat middens record for the Late Quaternary, P. remota appear to have expanded into the south of Edwards Plateau to the west or southwest (Van Devender 1990). Likewise, P. fallax in the Sonoran Desert may have expanded from California chaparral to Arizona across the Pleistocene-Holocene boundary (Betancourt et al. 1990).

Although some populations of P. californiarum and P. monophylla are found in limited sympatry in California (San Bernardino Co.), we did not detect gene flow in any direction, but the admixture analysis provided additional information about the relationship between P. californiarum and P. monophylla (Fig. 3A). Admixture analyses should be interpreted with care since it is difficult to estimate the real number of clusters, especially for species with long generation times like conifers; however, with the increased use of sequencing technologies and massive generation data, our capacity to detect admixture has substantially improved (Pritchard et al. 2000). The admixture analysis we carried out on the pinyon pines subclade from SW US and Baja California provided a picture of possible interbreeding in this group of pines. Distribution of ancestry fractions indicate that P. californiarum is introgressed but this scenario was not recovered in Phylonet. The admixture results indicate that P. quadrifolia shares little genetic variation with P. monophylla (s. s.). In addition, it stands out that P. quadrifolia is clustered in a singular population with K > 2. This pattern of ancestral structure suggests that P. quadrifolia accumulates particular genetic diversity and could be a valid species and not a hybrid. For K = 5 (data not shown) P. monophylla was not clustered into a separate population. It is also important to mention that according to clustering patterns P. fallax seems to share more genetic variation with P. edulis than P. monophylla. Clustering of P. fallax into a new well-defined population at a particular value of K could reflect enough genetic differentiation from P. edulis to support its isolation as a species. Although admixture patterns support in most cases our phylogenetic results, increased sampling is required for a more complete perspective of the admixture events for these complex and widely geographically distributed populations. Despite the reduced sample size, we were able to provide a preliminary insight into the admixture events occurring in this group of pines.

In conclusion, using target enrichment to characterize 304 nuclear genes and 26,499 SNPs, we corroborated the monophyly of Pinus subsection Cembroides in all analyses. The inferred phylogenies also corroborate other relationships previously recovered with plastid DNA. The results suggest that P. fallax and P. californiarum could be considered as valid species rather than as infraspecific taxa of P. monophylla or P. edulis. Also, our admixture results suggest that P. quadrifolia accumulates particular genetic diversity and could be a valid species and not a hybrid. The single-needle pinyons were recovered as a non-monophyletic group, with character reconstructions consistent with a single derivation and subsequent loss of the single-needle condition. The ASTRAL-III tree was consistent with the presence of ILS (very high) in the group of pinyon pines with small cones based on the short length in coalescent units of internal branches. Respecting reticulation events, we identified P. remota as having genes introgressed from P. fallax, and P. quadrifolia having genes introgressed from P. lagunae. Some hybridization scenarios were unexpected and not reported in the literature. Finally, further study is needed to determine the relative roles of ILS and introgression in explaining shared genetic diversity in Pinus subsection Cembroides.