INTRODUCTION {#s1}
============

Interspecies hybridization is relatively common across plant taxa and is thought to be the cause of several major speciation events ([@CIT0053]; [@CIT0001]). Hybridization tends to occur when species lack strong reproductive isolating mechanisms ([@CIT0005]), which normally act as barriers to reproduction and help in the process of speciation ([@CIT0050]). If species have loose genetic barriers (compatible zygotes) and at least partially overlap spatially and temporally, then hybrids are likely to form ([@CIT0065]). Intermediate forms can be detected in recent hybridization events ([@CIT0064]; [@CIT0013]); however, when backcrossing occurs in a hybrid system, morphological traits are present either as a continuum between the two parental species ([@CIT0027]), or more representative of one of the parents, resulting in cryptic introgression ([@CIT0046]; [@CIT0029]; [@CIT0045]). The use of morphological traits to determine hybridization events can be helpful when divergent features exist between parental species, but, when morphological limits are not well defined in taxonomically challenging groups, it becomes difficult to distinguish intermediacy from interpopulation variation ([@CIT0057]). For example, in some systems, hybridization leads to a mosaic of forms or even to extreme or novel characters ([@CIT0049]).

Hybridization in pines has been widely acknowledged ([@CIT0011], [@CIT0012]; [@CIT0062]; [@CIT0041], [@CIT0042]). The pinyon pine complex is an excellent system to study hybridization due to its lack of: strong reproductive isolating mechanisms ([@CIT0034]); taxonomic congruence ([@CIT0022], [@CIT0023]); and conclusive genetic studies ([@CIT0043]). Most pine species lack interspecific incompatibility mechanisms ([@CIT0011]) and have wind-mediated, long-distance pollen dispersals ([@CIT0060]) that can facilitate gene flow between allopatric groups (Wright, 1952). Most notably, the Parry pinyon pine, *Pinus quadrifolia* Parl. ex Sudworth, has been proposed to have a hybrid origin based on morphological features and geography ([@CIT0034]). *Pinus quadrifolia* commonly has four needles per fascicle, as its epithet suggests, and occurs from Riverside County, California, USA to northern Baja California, Mexico. Based on intermediate morphology, [@CIT0034]) hypothesized its origin as the result of interspecific hybridization between two species with similar pollen dispersal times ([@CIT0039]; [@CIT0017]) and partially overlapping distributions ([Fig. 1](#F1){ref-type="fig"}), *P. monophylla* Torrey & Fremont and *P. juarezensis* Lanner. At the time of his discovery, needle number was the main morphological trait used in classifying pinyon pines, so he reasoned that individuals with one needle per fascicle were one species (*P. monophylla*), individuals with five needles were another (*P. juarezensis*) and the three- and four-needled individuals (*P. quadrifolia*) were intermediate hybrids of the former two.

![Putative distribution map of *P. monophylla* (red), *P. californiarum* (blue) and *P. quadrifolia + P. juarezensis* (purple) ([@CIT0055]; [@CIT0009]). Sampling locations are represented by the overlaid shapes (circle = *P. monophylla*, triangle = *P. californiarum*, square = *P. quadrifolia + P. juarezensis*).](mcaa045f0001){#F1}

Several attempts have been made to resolve the hybrid nature *of P. quadrifolia*; however, they have provided inconclusive results. The idea of hybridization was explored using additional morphological features other than needle number, but these studies focused on the potential hybridization between *P. monophylla* and *P. quadrifolia* without considering *P. juarezensis* ([@CIT0064]; [@CIT0052]). These studies also explored chemical profiles, such as mono- and sesquiterpenoids, that could be used for species identification between *P. quadrifolia* and *P. monophylla* ([@CIT0064]); however, an intermediate chemical profile was detected in individuals in sympatry, suggesting hybridization. Additionally, the only genetic study including the three putative taxa (*P. juarezensis*, *P. quadrifolia* and *P. monophylla*) was done with internal transcribed spacer (ITS) sequences ([@CIT0022]), which have been proven to have paralogy in pines and other plants ([@CIT0024]; [@CIT0063]; [@CIT0019], [@CIT0018]). When using this nuclear marker, species were not monophyletic, so the authors then suggested a potential hybridization event between *P. monophylla* and *P. juarezensis* or incomplete lineage sorting as the cause. Therefore, the question of whether *P. quadrifolia* has a hybrid origin has remained unanswered for \>40 years.

Since Lanner's publication suggesting the hybrid origin of *P. quadrifolia* from the parental species *P. monophylla* and *P. juarezensis* ([@CIT0037]), several taxonomic studies have been conducted, some of them suggesting the splitting of *P. monophylla* into three different taxa based on morphology and niche variability: *P. monophylla*, *P. californiarum* and *P. fallax* ([@CIT0037]; [@CIT0004]; [@CIT0009]). One of these suggested taxa, *P. californiarum* D.K. Bailey, has more resin canals and occupies drier habitats in higher elevations than *P. monophylla*. In a recent study, hybridization between *P. californiarum* and *P. quadrifolia* was explored using a Hyb-seq method, yet no gene flow was detected ([@CIT0043]). One reason could be the lack of adequate sample size due to the taxonomic focus of the study which aimed to capture as many species as possible using only a few individuals. [@CIT0043] stated themselves that increased sampling would be required to draw accurate conclusions about admixture events in this system. In our study, we will be testing whether *P. juarezensis* hybridizes with *P. monophylla* or with the recently proposed taxa *P. californiarum.* In contrast, *P. juarezensis* has not been recognized as an independent species by some authors, but rather as a natural morphological variation of *P. quadrifolia*, able to produce 3--5 needles per fascicle. For example, [@CIT0017] challenged Lanner's conclusion, pointing to the extreme variation in needle number on a single tree and the year to year variation on trees reported by [@CIT0034]) himself. We will be testing the taxonomic validity of *P. juarezensis* as an independent species.

With advancements in next-generation sequencing (NGS), the ability to detect reticulated evolution through hybridization has drastically improved ([@CIT0047]; [@CIT0054]; [@CIT0040]). We are now able to revisit these hybridization studies with a finer scale, allowing us to examine the purported hybrid origin of *P. quadrifolia* and address the taxonomic standing of *P. juarezensis*. The goals of this study are: (1) to determine the taxonomic validity of *P. juarezensis* as a species; (2) to assess the hybrid origin of *P. quadrifolia*; and (3) to test for hybridization and intermediate morphology among species in this complex. We have included samples from the taxa of interest using nuclear genomic data as well as morphological data to assess these questions.

MATERIALS AND METHODS {#s2}
=====================

Sampling {#s3}
--------

Samples were taken from 13 locations across the Southwestern USA and Baja California, Mexico ([Table 1](#T1){ref-type="table"}; [Fig. 1](#F1){ref-type="fig"}), randomly sampling ten trees per putative species per site, at least 30 m apart from each other. The distribution of this complex extends from southern Idaho, USA to southern Baja California, Mexico, with a majority of *P. monophylla* occurring in western Utah, Nevada and central California, *P. californiarum* occurring in southern Nevada, southern California and down into Baja California, and *P. juarezensis* (five needle) and *P. quadrifolia* (3--4 needles) occurring in southern California and Baja California. These taxa are characterized by inhabiting small disjunct areas in southern California and Baja California. Ten centimetres of branch tips were cut from each tree, representing 2--3 years of growth and averaging approx. 50 fascicles per tree. Most collections were made from late 2017 to early 2019, with the exception of a few collections made previously (2011--2013). We collected a larger number of individuals in San Jacinto because [@CIT0034]) proposed this area as one of the primary hybrid zones between *P. juarezensis* and *P. monophylla.*

###### 

Sample locations with population co-ordinates and number of trees sampled per population

  Location                             Latitude    Longitude     No. of trees genotyped   Needle type
  ------------------------------------ ----------- ------------- ------------------------ --------------------------------------------------------------
  Markleeville, CA                     38.69894    --119.77082   10                       *P. monophylla*
  Manhattan, NV                        38.54101    --117.05981   6                        *P. monophylla*
  Mono Lake, CA                        37.92145    --119.06433   7                        *P. monophylla*
  Pine Mt, CA                          36.02982    --118.15251   9                        *P. monophylla*
  Mojave Desert, CA                    35.17397    --115.40707   13                       *P. californiarum* Hybrids
  San Jacinto, CA                      33.55887    --116.60994   37                       *P. quadrifolia* *P. californiarum* *P. juarezensis* Hybrids
  Laguna Mts, CA                       32.87554    --116.41017   9                        *P. quadrifolia, P. juarezensis*
  Jacumba, CA                          32.63317    --116.09395   15                       *P. quadrifolia* *P. californiarum* *P. juarezensis* Hybrids
  La Rumorosa, Baja CA                 32.521544   --116.04122   21                       *P. quadrifolia* *P. californiarum* *P. juarezensis* Hybrids
  San Salvador, Baja CA                31.742415   --115.97815   8                        *P. quadrifolia* *P. juarezensis*
  Lazaro Cardenas, Baja CA             31.257595   --115.59940   11                       *P. californiarum*
  Northern San Pedro Martir, Baja CA   31.034611   --115.46477   18                       *P. quadrifolia* *P. californiarum* *P. juarezensis*
  Southern San Pedro Martir, Baja CA   31.02131    --115.51488   10                       *P. quadrifolia* *P. juarezensis*

California (CA), Nevada (NV), Baja California (Baja CA).

Genomic clustering analyses {#s4}
---------------------------

In order to examine the taxonomic validity of *P. juarezensis* as a species, the hybrid origin of *P. quadrifolia* and hybridization among species of this complex, we first identified genetic clustering and admixed individuals. Nuclear DNA was extracted using the 2 % cetyltrimethylammonium bromide (CTAB) protocol ([@CIT0015]), quantified and sent to Diversity Arrays Technology (DArT), who produced a reduced library, and sequenced on an Illumina Hiseq2500 system (DArTseq). DArTseq™ represents a combination of DArT complexity reduction methods and NGS platforms ([@CIT0028]). Genome reduction is achieved by a combination of endonucleases that specifically target low-copy DNA areas, rather than repetitive DNA fragments ([@CIT0059]). This allows for detection of a high number of informative single nucleotide polymorphisms (SNPs) across the genome. The result is a genomic 'representation', comprising both constant and polymorphic fragments across individuals. NGS of these 'representations' reveals the sequence (approx. 70 bp) of an informative DNA fragment and each individual's state compared with all others, namely (1) homozygosity with the reference allele; (2) homozygosity with an alternative allele; or (3) heterozygosity, comprising both a reference and an alternative SNP allele. The technology was optimized for our taxa using combinations of enzymes (*Pst*I/*Hpa*II, *Pst*I/*Sph*I, *Sbf*I/*Hpa*II and *Sbf*I/*Mse*I) to select the most appropriate complexity reduction method, in terms of both the size of the representation and the fraction of a genome selected for assays. DArTseq™ has been successfully applied in genomic studies exploring species boundaries and hybridization in plants and animals ([@CIT0056]; [@CIT0040], [@CIT0051]). This technique enables genome-wide studies of non-model organisms, those for which there is limited genomic information. For our analysis, we selected *P. lambertiana* as the reference genome (GCA_001447015.2) due to its relatively close evolutionary relationship to the pinyon subsection *Cembroides*. A low-density DArTseq assay resulted in 18 518 SNPs. These SNPs have an index generated by reproducing the data independently, which is 'the proportion of technical replicate assay pairs for which the marker score is consistent' ([@CIT0026]). We filtered out all loci with reproducibility lower than 100 %, missing data lower than 15 %, all monomorphs, all loci departing from Hardy--Weinberg equilibrium and all but one locus where there was more than one locus per sequence tag, resulting in a final data set of 1868 loci. Input data were the metadata provided by DArTseq, saved as an xlsx file: '0' (homozygosity with reference allele); '1' (homozygosity with alternative allele); '2' (heterozygote) and '--', fragment missing in representation -- double null (absence of fragment with an SNP in genomic representation). The processed marker data were reformatted into appropriate file types for downstream analyses using the R program dartR ([@CIT0025]).

The 1868 loci were included in four complementary methods of genetic clustering: principal co-ordinates analysis (PCoA), discriminant analysis of principal components (DAPC), fastSTRUCTURE and ADMIXTURE, which also allowed us to detect admixed individuals between different genetic clusters. We then evaluated whether *P. juarezensis* (five-needled individuals) forms its own genetic cluster and whether *P. quadrifolia* (three- and four-needled individuals) is formed by admixed individuals from the parental species *P. monophylla* and *P. juarezensis* or whether there are other pairs of species hybridizing in this complex.

### PCoA {#s5}

A PCoA, which considers differences in allele frequencies between individuals, was performed in dartR (gl.pcoa.plot) to examine genetic distance among populations ([@CIT0026]).

### DAPC {#s6}

Population clustering was determined using a DAPC ([@CIT0032]) in R with the adegenet package ([@CIT0030]). DAPC uses a model-free k-means clustering algorithm on a principal components analysis (PCA)-transformed data set, which maximizes variation between groups while reducing the number of variables and computation time needed to identify existing genetic clusters ([@CIT0031]). Differing cluster solutions are compared using the Bayesian information criterion (BIC), with the lowest BIC score corresponding to the optimal cluster solution. The xvalDapc command, which uses a stratified random sampling to ensure one member of each population is represented in both a training set consisting of 90 % of the data and a validation set consisting of the remaining 10 % of the data, was used as a cross-validation method to determine the appropriate number of principal components to retain ([@CIT0031]).

### fastSTRUCTURE. {#s7}

A Bayesian analysis of population clustering was performed in the software fastSTRUCTURE ([@CIT0048]) using the logistic prior and five cross-validations. Output files include the mean Q value for each individual, defining the mean probability as belonging to any one of the populations K1 to Kx. Model complexity (*K*) was selected using the chooseK command built into fastSTRUCTURE. The resulting Q mean bar plots were visualized using the online application pophelper ([@CIT0020]). Results were confirmed using ADMIXTURE ([@CIT0002]) and a sub-structuring method, wherein samples with the same population identity (*K*) were reanalysed in fastSTRUCTURE using the logistic prior and five cross-validations to observe lower hierarchical structuring ([@CIT0048]).

Hybrid generation identification {#s8}
--------------------------------

### NewHybrids. {#s9}

In order to determine an individual's hybrid category, e.g. early F~1~ or advanced generation (F~2~ and backcross), we used a Bayesian model-based clustering method in the software NewHybrids 1.0 ([@CIT0003]). The program uses Markov chain Monte Carlo simulations to compute the posterior probability of an individual belonging to pre-defined ancestry categories, including pure, F~1~, F~2~ or backcrossed ([Table 2](#T2){ref-type="table"}). The program compares two parental genotypes at a time, so we created three data sets to represent each pairwise species cross by removing individuals with ancestry of the third genotype (e.g. for *P. monophylla* × *P. quadrifolia* analysis, we removed pure *P. californiarum* individuals and individuals with *P. californiarum* ancestry). Runs were initiated at different random starting points using the Jeffrey's prior for both theta and pi with a burn-in of 10 000 and 100 000 sweeps ([@CIT0010]).

###### 

Twelve genotype frequency categories input into NewHybrid analyses based on [@CIT0007]

  Genotype category                Cross type                       Expected ancestry proportions                   
  -------------------------------- -------------------------------- ------------------------------- ------- ------- -------
  Pure                             Species 1                        1                               0       0       0
  Pure                             Species 2                        0                               0       0       1
  F~1~                             Species 1 × Species 2            0                               0.5     0.5     0
  F~2~                             F~1~ × F~1~                      0.25                            0.25    0.25    0.25
  F~1~ Backcross 1                 F~1~ × Species 1                 0.5                             0.25    0.25    0
  F~1~ Backcross 2                 F~1~ × Species 2                 0                               0.25    0.25    0.5
  F~2~ Backcross 1                 F~2~ × Species 1                 0.5                             0.125   0.125   0.25
  F~2~ Backcross 2                 F~2~ × Species 2                 0.25                            0.125   0.125   0.5
  1 Backcross × F~1~ Backcross 1   Species 1 × (F~1~ × Species 1)   0.75                            0.125   0.125   0
  2 Backcross × F~1~ Backcross 2   Species 2 × (F~1~ × Species 2)   0                               0.125   0.125   0.75
  1 Backcross × F~2~ Backcross 1   Species 1 × (F~2~ × Species 1)   0.625                           0.125   0.125   0.125
  2 Backcross × F~2~ Backcross 2   Species 2 × (F~2~ × Species 2)   0.125                           0.125   0.125   0.625

Morphology {#s10}
----------

In order to determine which morphological features are consistent with the genetic clusters and to detect intermediate morphology associated with hybridization, we analysed a range of leaf morphological traits that have been used in most taxonomic studies of pinyon pines, including needle number, number of resin canals and number of rows of stomata ([@CIT0037]; [@CIT0004]; [@CIT0036]; [@CIT0039]; [@CIT0008]; [@CIT0009]; [@CIT0019]). All sample shoots were kept in the freezer at --20 °C and analysed individually at room temperature. All fascicles found in the 10 cm tip branches were visually examined for needle number. Most trees retain their needles, but some drop a few needles. To avoid any bias in needle number estimates, close examination of branch tips was done under a stereomicroscope to inspect fascicle scars, allowing us to detect dropped needles and accurately estimate the range of needle number per branch. Then we identified whether trees have the same number of needles across the branch (uniform number) or if they have a range of needle numbers (non-uniform). One fascicle was measured from branch tips with uniform needles per fascicle numbers. For branch tips with non-uniform needles per fascicle numbers, one fascicle of each needle number variant was measured. Clear nail polish was used when necessary to enhance the visualization of stomatal rows. A stomatal position index was created ranging from 1 to 5, where 1 means ventral position only, 2 is full ventral and partial dorsal, 3 is complete on both sides, 4 is full dorsal and partial ventral, and 5 is dorsal only. In addition to stomatal position, cross-sectional area and thickness of hypodermal layers were determined ([@CIT0019]). Every needle sample was cross-sectioned a few times with a razor blade in the midsection and rehydrated with water before measuring the cross-sectional area, dermal thickness and number of dermal layers. Measurements for cross-sectional area were traced in micrometres, and thickness of hypodermal layers was recorded as the average of four measurements around the needle. Variation in the number of resin canals in needles from different fascicles of the same branch was initially observed, thus two resin canal measurements were taken; however, the difference between these two measurements proved to be non-significant so only the first measurement was included in further analyses. All measurements were done using an ultra-high-resolution Nikon SMZ25 stereoscopic microscope zoom ×0.5--1.6 and NIS Elements software. We used Microsoft Excel to plot the means and standard errors of needle number by resin canals and stomatal rows by resin canals per population with and without genetic hybrids, following [@CIT0009]. This method allows us to examine differences among species while showing the intraspecific variation present among populations. It also allows us to examine the morphology of hybrids and compare them with general species characteristics.

RESULTS {#s11}
=======

Genetic clustering analyses {#s12}
---------------------------

The PCoA ([Fig. 2](#F2){ref-type="fig"}) shows individuals separating into three groups, with PCoA axis one representing 9.6 % of variation in genetic distance and PCoA axis two representing 7.9 %. While those two axes represent the two most informative axes, the subsequent axes incrementally represent less variation, with PCoA axis three representing 1.5 %, four representing 1.3 %, and so on. Each group represents individuals that have the morphology of *P. monophylla*, *P. quadrifolia*/*P. juarezensis* and *P. californiarum*, respectively. Individuals with allele frequencies in between *P. monophylla* and *P. quadrifolia*, *P. monophylla* and *P. californiarum*, and *P. quadrifolia* and *P. californiarum* can be seen falling out of the three main groups, suggesting possible admixture among the three taxa. It is important to note that a PCoA does not present statistical clusters, but instead functions as a tool to visualize differences in allele frequency and therefore must be interpreted with caution. Individuals on the PCoA graph may group together due to similarities in allele frequencies, but a PCoA does not provide a group assessment as in K-means clustering methods duch as DAPC and fastSTRUCTURE.

![Principal co-ordinates analysis showing the variation in genetic distance with three main groups circled by genetic identity. Each black dot represents an individual; dots outside of circles represent potential hybrids and hybrid derivatives.](mcaa045f0002){#F2}

The DAPC analysis ([Fig. 3](#F3){ref-type="fig"}) of the SNP data set resulted in a value of *K* = 3 by K-means clustering, retaining 20 principal components and two discriminant functions. All sample populations of *P. monophylla* clustered together as one group (coloured in red). All sample populations of *P. quadrifolia* and purported *P. juarezensis* clustered together as one group (coloured in purple). All sample populations of *P. californiarum* (coloured in blue) clustered together as one group.

![Discriminant analysis of principal components showing three distinct genetic clusters coloured by genetic identity. Each dot represents an individual connected to the genetic cluster's centroid by a line; inertia ellipses indicate their assignment to one of the three genetic clusters inferred by DAPC.](mcaa045f0003){#F3}

Our analyses of fastSTRUCTURE show three distinct genetic clusters using the chooseK method and visualized in pophelper ([Fig. 4](#F4){ref-type="fig"}; [Supplementary data Fig. S1](#sup1){ref-type="supplementary-material"}), supporting the DAPC and PCoA results. Again, *P. monophylla* forms one cluster (coloured in red), *P. quadrifolia* and purported *P. juarezensis* another cluster (coloured in purple) and *P. californiarum* the third cluster (coloured in blue). These results show no genetic distinction between individuals with morphological characteristics of *P. quadrifolia* (three to four needles) and *P. juarezensis* (five needle); this was further confirmed by the sub-structuring analysis of the cluster containing only the putative *P. juarezensis* and *P. quadrifolia* ([Supplementary data Fig. S2](#sup2){ref-type="supplementary-material"}). Moreover, this analysis did not support the hypothesized hybrid origin of *P. quadrifolia*. However, these results do show extensive admixture among the three taxa, with several examples of *P. monophylla* × *P. quadrifolia* (as in San Jacinto), *P. monophylla* × *P. californiarum* (as in Mojave) and *P. quadrifolia* × *P. californiarum* (as in Rumorosa, Mojave, San Jacinto and Jacumba) admixture events. These admixture events are evident in the fastSTRUCTURE plot ([Fig. 4](#F4){ref-type="fig"}) where individuals have more than one genetic cluster assigned to them. ADMIXTURE analyses supported three genetic clusters with the lowest cross-validation error at *K* = 3 ([Supplementary Data Fig. S3](#sup3){ref-type="supplementary-material"}). Sub-structuring using fastSTRUCTURE's logistic prior and ten cross-validations found no lower hierarchical structuring below *K* = 3.

![fastSTRUCTURE plot showing three genetic clusters (*K* = 3) coloured by genetic identity. Each line on the *x*-axis represents an individual, and the probability of those individuals belonging to a certain genetic cluster is represented by the proportion of colour on the *y*-axis. Lines with multiple colours represent individuals with admixture from multiple genetic clusters. Cartoons showing the number of needles per fascicle and cross-sectional pictures showing the resin canals are overlaid on their corresponding genetic groups.](mcaa045f0004){#F4}

Hybrid generation identification {#s13}
--------------------------------

Interestingly, the NewHybrids analyses ([Table 3](#T3){ref-type="table"}) show that most individuals initially identified as having hybrid ancestry in fastSTRUCTURE end up grouping as advanced generation backcrosses, suggesting that F~1~ hybrids are able to reproduce further. Several individuals had multiple category assignments with probability scores lower than 0.98, suggesting uncertainty in assignment ([@CIT0003]). These pure and mixed assignment individuals are possibly backcrosses of a more advanced generation than our genotype categories can detect. Individuals that appeared near a Q value of 0.4 on the fastSTRUCTURE plot were at least third- or fourth-generation hybrid derivatives, with the most recent hybrid ancestry generation of F~2~ × *P. quadrifolia* and no F~1~ or F~2~ individuals detected. A majority of the individuals were identified as fourth-generation backcrosses, having a genotype class of Backcross × (F~2~ × Backcross).

###### 

Levels of ancestry determined by NewHybrids for the apparent hybrid derivatives indicated in fastSTRUCTURE analyses.

  NewHybrids assignment                              Contributing species                   No. of individuals
  -------------------------------------------------- -------------------------------------- --------------------
  Pure *P. quadrifolia*                              *P. quadrifolia*                       61
  Pure *P. monophylla*                               *P. monophylla*                        32
  Pure *P. californiarum*                            *P. californiarum*                     53
  *P. quadrifolia* × (F~2~ × *P. quadrifolia*)       *P. monophylla*, *P. quadrifolia*      7
  *P. quadrifolia* × (F~2~ × *P. quadrifolia*)       *P. californiarum*, *P. quadrifolia*   4
  F~2~ × *P. Quadrifolia*                            *P. monophylla*, *P. quadrifolia*      1
  F~2~ × *P. Quadrifolia*                            *P. californiarum*, *P. quadrifolia*   1
  *P. californiarum* x (F~2~ × *P. californiarum*)   *P. monophylla*, *P. californiarum*    14

Morphological analysis {#s14}
----------------------

In order to properly assign morphological traits to species, characteristics were assigned to the genetic groups resulting from the genetic structuring analyses, particularly fastSTRUCTURE and NewHybrids. Individuals assigned exclusively to one genetic cluster were determined to be pure, while individuals of mixed assignment were run through NewHybrids and distinguished as hybrids (F~1~) and hybrid derivatives (F~2~, backcrosses or advanced generation backcrosses). All morphological results are listed in [Table 4](#T4){ref-type="table"} and distinguishing traits are plotted in [Fig. 5](#F5){ref-type="fig"}. Morphological results also supported three groups, with *P. monophylla* and *P. californiarum* having one needle per fascicle on average but separated by the larger number of resin canals in *P. californiarum.* In contrast, *P. quadrifolia + P. juarezensis* formed one group, having the highest number of needles but the lowest number of rows of stomata, and stomata just on the ventral surface. Intermediate forms were detected when hybrids and hybrid derivatives were plotted ([Fig. 5B](#F5){ref-type="fig"}, [C](#F5){ref-type="fig"}); however, in some cases, some individuals more closely resembled one of the parental species.

###### 

Morphological characteristics based on genetic identity

  Genetic identity     Needle per fascicle          Resin canals           Stomatal rows   Stomatal position          Area (mm^2^)                No. of dermal layers   Thickness of dermal layers (µm)
  -------------------- --------------------- ------ -------------- ------- --------------- -------------------------- --------------------------- ---------------------- ---------------------------------
  *P. monophylla*      1.06 ± 0.24           1--2   2.82 ± 1.08    1--6    24.48 ± 4.09    Complete                   1.9 × 10^3^ ± 7.0 × 10^2^   2.36 ± 0.65            59.08 ± 13.86
  *P. quadrifolia*     4.19 ± 0.77           2--5   1.91 ± 0.78    0--7    7.69 ± 1.80     Ventral only               7.9 × 10^2^ ± 2.5 × 10^3^   2.59 ± 0.55            36.76 ± 37.57
  *P. californiarum*   1.11 ± 0.31           1--2   12.30 ± 2.96   6--23   19.81 ± 4.02    Complete                   1.1 × 10^3^ ± 2.1 × 10^3^   2.49 ± 0.50            37.43 ± 25.25
  M × Q                3.71 ± 1.03           2--5   2.04 ± 0.26    2--3    8.71 ± 1.98     Ventral only or complete   5.6 × 10^2^ ± 3.6 × 10^2^   2.14 ± 0.64            40.37 ± 13.90
  Q × C                2.83 ± 1.40           1--5   4.22 ± 4.60    1--17   11.58 ± 3.90    Ventral only or complete   5.9 × 10^2^ ± 4.0 × 10^2^   2.66 ± 0.47            38.49 ± 18.42
  C × M                1.13 ± 0.33           1--2   8.94 ± 2.65    1--12   16.92 ± 2.56    Complete                   2.0 × 10^3^ ± 2.4 × 10^3^   2.97 ± 0.62            51.38 ± 12.86
  M × Q × C            1.50 ± 0.50           1--2   8.50 ± 0.50    8--9    14.00 ± 0.00    Complete                   8.8 × 10^2^ ± 9.5 × 10^1^   2.00 ± 0.00            51.38 ± 0.84

Note that *P. juarezensis* is treated as a synonym of *P. quadrifolia* because no genetic differences were found between them.

Letters in the hybrid cross rows represent the first initial of the species name (M = *P. monophylla*, Q = *P. quadrifolia*, C = *P. californiarum*).

![Morphological analyses showing in (A and B) needle number by resin canals and (C and D) stomatal rows by resin canals, with (A and C) showing pure species, and (B and D) showing pure species and hybrid derivatives. Each dot represents the mean of a population with standard error bars; each colour corresponds to a species (red circle = *P. monophylla*, blue triangle = *P. californiarum*, purple square = *P. quadrifolia*, diamonds = hybrid derivatives).](mcaa045f0005){#F5}

DISCUSSION {#s15}
==========

Genetic clustering analyses show three groups: *P. monophylla*, *P. californiarum* and *P. quadrifolia*. Five-needled individuals are not genetically distinguishable from three- and four-needled individuals, contrary to Lanner's original hypothesis. Individuals with hybrid ancestry were detected; however, the hypothesized hybrid origin of *P. quadrifolia* from the species tested was not supported. Morphological analyses showed a clear distinction between *P. californiarum* and *P. monophylla* in the number of resin canals, while both *P. californiarum* and *P. monophylla* were distinguishable from *P. quadrifolia* in the number of needles per fascicle. *Pinus quadrifolia* showed variation in the number of needles per fascicle throughout its range, with most individuals having a variable number of needles on the same branch.

Lanner's hypotheses on the existence of *P. juarezensis* and the hybrid origin *P. quadrifolia* proved to be incorrect. While *P. quadrifolia* does have variation in its needle number, this variation is not explained by interspecific hybridization. *Pinus juarezensis* is not a valid taxon, but rather a synonym of *P. quadrifolia*, as five-needled individuals genetically cluster with the three- and four-needled individuals. Hybridization was detected between *P. quadrifolia* and *P. monophylla*, as well as between *P. quadrifolia* and *P. californiarum*. Interestingly, a new admixture combination was detected between *P. californiarum* and *P. monophylla*. The morphological traits of recent hybrid derivatives seemed to be intermediate between the parental species, yet, as apparent backcrossing occurred, the advanced generation hybrid derivatives appeared morphologically indistinguishable from their main genetically contributing parent.

Genetic clustering analyses do not support *P. juarezensis* as a valid taxon {#s16}
----------------------------------------------------------------------------

Previous studies have either failed to include *P. juarezensis* in their analyses or found inconclusive results regarding its taxonomic validity ([@CIT0064]; [@CIT0052]; [@CIT0022]; [@CIT0043]). In our study, we included several populations with the five-needled pinyon individuals, including some from Lanner's original discovery ([@CIT0034]). Genetic analyses, including PCoA, DAPC, fastSTRUCTURE and ADMIXTURE, showed three distinct genetic groups, with no difference between *P. quadrifolia* and the purported *P. juarezensis*. This may initially appear to align with Lanner's hypothesis that the five-needled taxon would eventually be introgressed out of existence ([@CIT0034]); however, purely five-needled trees were sampled and did not have a distinct genetic structure from the three- and four-needled samples identified as *P. quadrifolia*, implying that they are all one species. With pine trees' long generation time and slow mutation rate ([@CIT0006]; [@CIT0014]), it is not plausible that a species could have disappeared from the gene pool since Lanner's observation ([@CIT0034]) as a result of genomic swamping. It is more likely that *P. quadrifolia* can express 3--5 needles per fascicle. Additionally, purely five-needled trees were so rare compared with the mixed three-, four- and five-needled individuals that they seem to be more of an exception, rather than the standard for a species. Lanner himself reported such rarity in his initial study, with only three individuals having 100 % five needle coverage for 1 year but then expressing three or four needles in subsequent years ([@CIT0034]). With all this evidence considered, we reduce *P. juarezensis* to a synonym of *P. quadrifolia*, supporting the hypothesis of [@CIT0017].

*Pinus quadrifolia* is not of *P. monophylla* × *P. juarezensis* hybrid origin {#s17}
------------------------------------------------------------------------------

The unique needle number variation (i.e. 3--5 needles per fascicle) present in *P. quadrifolia* populations has led many to question its origin. Hybridization events between other pinyon species, such as *P. monophylla* and *P. edulis* for example, produce individuals with mixed needle number ([@CIT0035]; [@CIT0036]; [@CIT0009]). Previous research explored the possibility that the same mechanisms could be producing such variation in *P. quadrifolia*. For the first time, our study used NGS methods to test for the hybrid origin of *P. quadrifolia* suggested by the intermediate morphology between five-needled and one-needled pinyon pines. Genetic analyses showed no evidence that *P. quadrifolia* was of *P. monophylla* × *P. juarezensis* hybrid origin, but instead has a unique genetic identity. FastSTRUCTURE results showed *P. quadrifolia* as a distinct genetic group, with most individuals lacking genetic admixture. Individuals with a combination of three, four and five needles were most common, yet showed no genetic distinction from the individuals with only four or only five needles. We therefore reject Lanner's declassification of the Parry pinyon to the hybrid epithet *P. × quadrifolia* and support its species status as *P. quadrifolia*. The morphological variation originally seen by Lanner and others is probably not due to hybridization events but may be explained by intraspecies variation or environmentally induced plasticity. In fact, in several tree stands where *P. quadrifolia* and *P. californiarum* grow together, the needle variation is not increased compared with isolated pure *P. quadrifolia* stands but instead this needle number range occurs throughout its distribution. This variation could be environmentally related, yet *P. californiarum* shares a large portion of its distribution with *P. quadrifolia* and does not seem to have nearly as much variation; however, the mode of monophylly inheritance is unknown ([@CIT0021]).

Tridirectional hybridization: detection of hybridization among the valid taxa *P. quadrifolia*, *P. californiarum* and *P. monophylla* {#s18}
--------------------------------------------------------------------------------------------------------------------------------------

Our data support the existence of three taxa: *P. quadrifolia*, *P. californiarum* and *P. monophylla.*[@CIT0034]) treated the one-needled individuals from Southern California and Baja California (San Jacinto Mountains, Rumorosa and San Pedro Martir) as *P. monophylla*. We found that these populations are actually *P. californiarum* based on the genetic data. Furthermore, the morphology of these populations matches the morphology of other *P. californiarum* populations, with a high number of resin canals. We found no evidence of *P. monophylla* in Baja California. While *P. quadrifolia* is not of *P. monophylla* × *P. juarezensis* hybrid origin, it hybridizes with *P. californiarum* and *P. monophylla*. Previous studies hypothesized that hybridization occurred between *P. quadrifolia* and *P. monophylla* ([@CIT0064]; [@CIT0052]; [@CIT0033]); however, we also discovered genetic admixture between *P. quadrifolia* and *P. californiarum*. These two species share a significant amount of their range and even have stands where they grow a few metres apart. It appears that there are no strong genetic/post-zygotic barriers preventing hybrids of these species from forming and proliferating within natural populations. However, more extensive pre- and post-zygotic studies examining shifted phenology and outbreeding depression are needed to understand why some populations have more hybrid derivatives than others. A few instances of *P. quadrifolia* hybridizing with *P. monophylla* were seen in San Jacinto, and, while these species' ranges do not overlap geographically, their nearest distributions are only 30 miles apart, a feasible distance for both pollen and seeds dispersed by corvids to travel ([@CIT0058]; [@CIT0060]). This southern direction of gene flow would correspond to the Santa Ana winds that occur through October to May ([@CIT0044]), well within the pollen dispersal times of *P. monophylla* and *P. quadrifolia* around May ([@CIT0039]; [@CIT0017]). Given the lower frequency of *P. monophylla* genes present in *P. quadrifolia*, these individuals may be relics from the early Holocene when *P. monophylla* had a more restricted southernly range, a pattern of distribution that has been detected in some populations of *P. monophylla* and its close relative *Pinus edulis* ([@CIT0038]; [@CIT0016]). This contraction in distribution has been considered the source of hybridization in the latter ([@CIT0016]).

Surprisingly, a majority of the introgression detected in this study was between two species not previously thought to hybridize, *P. monophylla* and *P. californiarum*. Originally considered synonymous, the potential of these two species for hybridization has long been overlooked, possibly due to their identical needle number. They share a significant portion of their range in southern Nevada and southern California, but differ quite drastically in their number of resin canals. Curiously, in the site where hybrid derivatives are present (Mojave National Park), only one *P. californiarum* individual was detected, indicating how readily these two species hybridize. This could additionally imply heterosis, with increased fertility of the F~1~ generation allowing for frequent backcrossing, or selection potentially favouring hybrid phenotypes, promoting the proliferation of individuals with hybrid ancestry. Future studies on the fitness of *P. monophylla* × *P. californiarum* crosses must be done to explore these hypotheses.

Morphological traits in pure species vs. hybrids {#s19}
------------------------------------------------

Our study delimited species by genetic clusters first, without any *a priori* bias towards morphological traits. We then generalized the average characteristics of each genetic group to assign morphological qualities to each species. In doing so, our analyses were not influenced by strict morphological cut-offs and we were able to more accurately identify individuals with hybrid ancestry and intraspecific variation compared with other studies in pinyon pines. It is important to note that while some distinguishing morphological patterns were found, these should be used as suggestions for species identification, rather than strict species guidelines/boundaries.

Interestingly, the morphology of all three hybrid types largely depends on the species contributing more of its genotype. For example, in *P. quadrifolia × P. californiarum* hybrid derivatives, individuals with a majority of *P. californiarum* genotype tend to have one to two needles per fascicle with over ten resin canals; whereas individuals with a majority of *P. quadrifolia* genotype tend to have two to four needles per fascicle and two to four resin canals. This implies that needle morphology is largely genetically controlled and may be less influenced by environment than previously thought. The only exceptions to this were in recent hybrid derivatives, whose traits were more intermediate between their parents. For example, one *P. monophylla* × *P. quadrifolia* hybrid derivative found in San Jacinto had exclusively two needles per fascicle, 12 stomatal rows and two resin canals, resembling *P. edulis*. It is possible that as these individuals of hybrid ancestry further backcross, the intermediate traits could be recombined or introgressed out to more closely represent one of the parental species. This unfortunately makes the identification of individuals with advanced generation hybrid ancestry quite difficult without using genetic methods because these hybrid derivatives are often indistinguishable from their main genetically contributing parent, concordant with other hybrid systems ([@CIT0049]).

This study provides an insight into the morphological and genetic outcomes of species undergoing gene flow. We plan to further explore these outcomes and patterns of hybridization in pinyon pines, especially in our newly discovered admixture combination of *P. monophylla* and *P. californiarum*. While our study addresses long-standing taxonomic issues within a small group of pinyon pines, the NGS methods and genetics first approach we used can be applied to a variety of taxa with controversial standings and cryptic introgression.

SUPPLEMENTARY DATA {#s20}
==================

Supplementary data are available online at <https://academic.oup.com/aob> and consist of the following. [Figure S1](#sup1){ref-type="supplementary-material"}: fastSTRUCTURE bar plots of a growing number of genetic clusters, with the software-determined best K framed. [Figure S2](#sup2){ref-type="supplementary-material"}: the sub-structuring of *P. quadrifolia* and *P. juarezensis* samples in fastSTRUCTURE, showing one genetic cluster. [Figure S3](#sup3){ref-type="supplementary-material"}: the ADMIXTURE analysis of all samples, showing three genetic clusters corresponding to *P. monophylla*, *P. californiarum* and *P. quadrifolia* + *P. juarezensis*.

###### 

Click here for additional data file.

###### 

Click here for additional data file.

###### 

Click here for additional data file.

FUNDING {#s0104}
=======

This work was supported by the Hispanic-Serving Institutions Education Grants (HSI) Program \[grant no. 2018-38422-28614/project accession no. 1016839\] from the USDA National Institute of Food and Agriculture, by the National Science Foundation \[grant no. DEB-0816675\] and by the University Grants Program of San Diego State University.

We thank the Forest Service, Mojave National Park, San Pedro Martir National Park, Drs Sula Vanderplank, Jose Delgadillo, Rhett Michelson, Amy Whipple, Nuri Benet-Pierce, Rhett Michelson, Jamie Lamit, Hillary Cooper, Jean-Guillaume Lonjaret, Alexandra McElwee-Adame, Amy Orduño-Baez, Kyle Gunther, Mario Galvan, Stephanie Stragier and Brian Myers for facilitating plant collections, Katya Geissler for her help in leaf morphology, The Binational Species Group, the Flores-Renteria lab for feedback on the manuscript, and the two referees Vikram Chhatre and anonymous for their insight and helpful comments. R.B. and L.F.R. designed the study, collected plant material, conducted lab work, analysed and interpreted data and wrote paper; S.H. and A.H. collected morphological data and helped write the Materials and Methods section. All authors provided feedback on the manuscript.
