The evening primrose family (Onagraceae) comprises 664 species in 22 genera (Wagner et al. 2007; Wagner and Hoch 2021). Taxa in this family, and the large tribe Onagreae in particular, have a diverse array of plant-insect interactions that have been central to investigations into the evolution of floral traits (e.g. Raven 1979; Kawano et al. 1995; Raguso and Pichersky 1995). Within tribe Onagreae, the most species-rich genus is Oenothera L., which is subdivided into 18 sections including Oenothera sect. Pachylophus (Spach) W.L.Wagner, the focus of ongoing studies of plant-insect co-evolution and floral fragrance (Artz et al. 2010; von Arx et al. 2012; Skogen et al. 2016; Rhodes et al. 2017).

The most recent taxonomic revision of Oenothera sect. Pachylophus (Wagner et al. 1985) restricted the section to the five species and five subspecies that are currently recognized. Taxa in this section (Fig. 1) are distinguished from the rest of Oenothera by white petals, tuberculate capsules, seeds with unique hollow seed collars, lignified and compressed mesotesta, and an exotegmen that is (2–)3(–4) cells thick (Wagner 2005; Wagner et al. 2007). Oenothera sect. Pachylophus has a geographic range spanning western North America from southern Alberta and Saskatchewan, Canada to northern Chihuahua, Mexico and contains perennial or rarely annual, vespertine flowering plants (Wagner et al. 1985; Wagner 2005). Of the five recognized species, four (Oenothera cavernae, Oenothera brandegeei, Oenothera psammophila, and Oenothera harringtonii) are edaphically, geographically, or ecologically specialized narrow endemics. The fifth species, Oenothera cespitosa, is much more widespread (Fig. 2) and highly polymorphic, with subspecies partitioning the intermountain west of North America in large part by soil type. As such, O. cespitosa has been divided into five subspecies (O. cespitosa subsp. cespitosa, O. cespitosa subsp. crinita, O. cespitosa subsp. macroglottis, O. cespitosa subsp. marginata, and O. cespitosa subsp. navajoensis). Three of the five species in O. sect. Pachylophus are self-incompatible, obligate outcrossers (O. cespitosa, O. harringtonii, and O. psammophila; Wagner et al. 1985). These species have large petals and a strong fragrance that attracts long-tongued, crepuscular hawkmoths and facilitates pollination. The remaining two species, O. brandegeei and O. cavernae, are self-compatible and autogamous with smaller, putatively unscented flowers (Wagner et al. 1985).

Wagner et al. (1985) conducted the most thorough study to date of the systematics and evolutionary relationships in O. sect. Pachylophus, inferring phylogenetic relationships from 36 morphological traits. Results suggested that three of the four narrow endemic species were nested within the O. cespitosa species complex, whereas the fourth, O. psammophila, was hypothesized to be the result of reticulate evolution via hybridization between O. cespitosa subsp. cespitosa and O. cespitosa subsp. marginata (Wagner et al. 1985). This phylogeny, along with analyses of morphology, hybridization, cytology, and natural history, is the basis for current taxonomic circumscriptions within O. sect. Pachylophus. However, resolution of evolutionary relationships from these data was limited, with phylogenetic signal possibly obscured by processes such as ancient hybridization and convergent evolution. More recently, the evolutionary history of O. sect. Pachylophus was reconstructed using molecular data from two nuclear gene regions (internal and external transcribed spacer regions), but relationships remained weakly resolved (R. Levin unpubl. data) due to limited taxon sampling and insufficient phylogenetic signal.

Fig. 1.

Oenothera sect. Pachylophus. A. Oenothera cavernae. Habit with flower and immature capsule. Lincoln Co., Nevada (image by James M. André in 2006). B. Oenothera cespitosa subsp. crinita. Habit with flower. Inyo Co., California (image by Posnerk in 2018). C–D. Oenothera cespitosa subsp. navajoensis. Emery Co., Utah (images by W. L. Wagner in 2004). C. Stem with leaves, flower, and dehisced capsules from previous years. D. Close-up of stem with dehisced capsules. E–F. Oenothera cespitosa subsp. macroglottis. Boulder Co., Colorado. E. Habit with flowers (image by R. A. Raguso in 2012). F. Immature capsules (image by W. L. Wagner in 1979). G–H. Oenothera cespitosa subsp. marginata. G. Habit with flowers. Pima Co., Arizona (image by W. L. Wagner in 2003). H. Dehisced capsule. Sandoval Co., New Mexico (image by W. L. Wagner in 1979). I. Oenothera harringtonii. Habit. Fremont Co., Colorado (image by R. Bunn in 2018). J. Oenothera psammophila. Habit with flowers. Fremont Co., Idaho (image by R. A. Raguso in 2014).

Fig. 2.

Geographic range maps for taxa in Oenothera sect. Pachylophus. A. All taxa. B. Narrow endemic species (O. brandegeei, O. cavernae, O. harringtonii, and O. psammophila). C. Oenothera cespitosa subsp. cespitosa. D. Oenothera cespitosa subsp. crinita. E. Oenothera cespitosa subsp. macroglottis. F. Oenothera cespitosa subsp. marginata. G. Oenothera cespitosa subsp. navajoensis. Ranges are colored by taxon. Solid dots indicate exact sample locations, and unfilled dots indicate approximate sample locations in this study. Dots are labeled with the individual identifier for that sample (see Appendix 1).

Difficulties in inferring phylogenetic relationships among recently diverged species are not unique to Oenothera sect. Pachylophus. Recent studies have used high-throughput sequencing technologies to extensively sample nuclear gene regions to resolve historically recalcitrant nodes and generate more robust phylogenetic results (e.g. Gernandt et al. 2018; Stubbs et al. 2018; Couvreur et al. 2019; White et al. 2019). Specifically, hybridization-based target gene enrichment has developed into a valuable phylogenetic tool to target and sequence variable sites across the genome (e.g. Mamanova et al. 2010; Davey et al. 2011; Cronn et al. 2012; De Sousa et al. 2014; Mandel et al. 2014; Weitemier et al. 2014; Villaverde et al. 2018). This reduced-representation approach relies on the hybridization of genomic DNA libraries to biotinylated oligonucleotide probes, or baits, designed from genes of interest. In plants, highly conserved, single-copy exons are often targeted to allow probes designed for one species to capture orthologous regions in related species, reducing the cost and effort to obtain genome-wide data (Lemmon et al. 2012; McCormack et al. 2013; Johnson et al. 2019). Such techniques have been shown to effectively and accurately resolve evolutionary relationships at the order, family, and genus level within plant lineages (e.g. Mandel et al. 2014; Weitemier et al. 2014; Stephens et al. 2015; Heyduk et al. 2016; Sass et al. 2016). In the genus Inga (Fabaceae), Nicholls et al. (2015) were also able to resolve population-level relationships within one species, Inga umbellifera. Additionally, the recovery of flanking intronic regions can increase the number of phylogenetically informative sites available for analysis (McCormack et al. 2013; Weitemier et al. 2014; Heyduk et al. 2016; Johnson et al. 2019), which may allow for further elucidation of fine-scale relationships.

This study uses a target enrichment approach to reconstruct evolutionary relationships within Oenothera sect. Pachylophus. In doing so, we 1) explore the utility of combining sampling of multiple individuals within each taxon and target gene enrichment at low taxonomic levels, 2) assess the monophyly of taxa as they are currently defined, and 3) examine evolutionary relationships among taxa within Oenothera sect. Pachylophus. Further, the resulting phylogenetic framework is interpreted in light of geography (range size and overlap), edaphic characteristics, and life history (ecological and reproductive similarity) to better understand the drivers of speciation within Oenothera sect. Pachylophus.

Materials and Methods

Taxon Sampling—All species and subspecies in Oenothera sect. Pachylophus were sampled including 38 O. cespitosa individuals representing all currently recognized subspecies (11 O. cespitosa subsp. marginata, one O. cespitosa subsp. macroglottis, seven O. cespitosa subsp. navajoensis, 11 O. cespitosa subsp. cespitosa, and eight O. cespitosa subsp. crinita), four O. harringtonii individuals, one O. brandegeei individual, three O. cavernae individuals, and four O. psammophila individuals (Appendix 1). The number of individuals sampled was roughly proportional to the size of the geographic range of that taxon, and samples were chosen to encompass as much of the taxon range as possible (Fig. 2). One individual of O. triloba, which is outside Oenothera sect. Pachylophus, was included as the out-group based on genus-wide relationships (Overson et al. in mss.). All leaf tissue collected in the field was silica dried and stored at -20°C at the Chicago Botanic Garden until extraction. For individuals grown from seeds, fresh leaf tissue was used for extractions. Dried leaf tissue from herbarium sheets was also extracted for several taxa.

Library Construction, Enrichment, and Sequencing—Genomic DNA was extracted from leaf material using a modified CTAB protocol (Doyle and Doyle 1987). All genomic DNA was quantified using a Qubit 2.0 Fluorometer with a dsDNA HS Assay Kit (Invitrogen, Carlsbad, California) and stored at -20°C. Extracted DNA was fragmented using a Covaris M220 (Covaris, Woburn, Massachusetts), and 500 bp insert libraries were prepared using the TruSeq Nano HT DNA Library Preparation Kit (Illumina, Inc., San Diego, California) according to manufacturer instructions, but following fragmentation all reagents were used at half the recommended volume. Libraries were validated using an Agilent 2100 Bioanalyzer (Agilent Technologies, Inc., Santa Clara, California) and a Qubit 2.0 Fluorometer with a dsDNA HS Assay Kit (Invitrogen, Carlsbad, California). Target sequences were selected, and baits were designed ( https://doi.org/10.5281/zenodo.4437049) using two Oenothera species as described in Cooper et al. (in mss.); the baits were designed to target 322 conserved exon regions across the nuclear genomes of sampled taxa. In-solution target enrichment was performed using a MYbaits custom targeted enrichment kit (Arbor Biosciences, Ann Arbor, Michigan). Before hybridization, enriched libraries were normalized to 100 ng and pooled with 10–12 other samples. Manufacturer instructions were followed for hybridization; however, one quarter of the suggested bait volume was used, and the reaction volume was preserved using PCR grade water. Hybridization reactions were incubated at 65°C for 21 hrs, and 10 µl of the enriched library was PCR amplified for 16 cycles. Libraries were cleaned using a QiaQuick PCR Purification Kit (Qiagen, Hilden, Germany); excess adaptor, as detected with an Agilent 2100 Bioanalyzer (Agilent Technologies, Inc., Santa Clara, California), was removed using Ampure beads (Beckman Coulter, Beverly, Massachusetts). Libraries were pooled and diluted to a final concentration of 17 pM. Enriched libraries were sequenced on an Illumina MiSeq with paired-end 300 bp reads (2 × 300) and a 1% PhiX v3 sequencing control (Illumina, Inc., San Diego, California) at the Field Museum, Chicago, Illinois.

Sequence Processing and Alignment—Illumina data were processed with Trimmomatic to remove Illumina adaptor sequences and trim and filter reads by quality (Bolger et al. 2014). Any base below a Phred score of 10 was removed from either end of the sequence. The sliding window option was used to trim the sequence in areas where the average quality of four sequential base pairs was less than a Phred score of 20. Furthermore, any reads shorter than 20 bp were removed and not considered in further analyses.

The HybPiper pipeline was used to map reads, assemble contigs de novo, and identify targeted exonic and flanking non-coding regions ( https://github.com/mossmatters/HybPiper; Johnson et al. 2016). This generated a nucleotide sequence including both introns and exons, as well as an inferred protein sequence for the coding region of each gene (Johnson et al. 2016). We used a script (gene_recovery_heatmap.R) available within HybPiper to identify and remove any sequence that was less than 25% of the total reference (target) gene length. Following this pruning step, the remaining individuals retained 109–307 loci, with an average of 283 loci. MAFFT was used to align the inferred protein sequences across all samples using the L-INS-I option with default settings (Katoh and Standley 2013). Codon alignments were generated by forcing nucleic acids to the protein alignments using Pal2Nal with default settings (Suyama et al. 2006). Finally, trimAl was used to remove any position in the alignment that was missing a base in more than 50% of taxa (Capella-Gutiérrez et al. 2009).

Phylogenetic Analyses—Both concatenated and coalescent-based phylogenetic analyses were performed. For the concatenated phylogenetic reconstruction, the 307 aligned genes were concatenated into a single supermatrix with a Python script (fasta_merge.py) from HybPiper (Johnson et al. 2016). This concatenated supermatrix had a final aligned length of 473,892 bp, of which 35,652 bp were parsimony informative. Maximum likelihood analyses were performed using RAxML HPC v8.2.0 (Stamatakis 2014) on XSEDE through the CIPRES Science Gateway (Miller et al. 2010). A species tree was generated with 100 bootstraps using rapid bootstrapping, the GTRCAT model of heterogeneity for nucleotide substitution, and all other parameters as default. For coalescent-based phylogenetic reconstruction, each aligned locus was partitioned by gene and codon position, and RAxML HPC v8.2.0 (Stamatakis 2014) was used to generate gene trees with the same parameters specified above. The RAxML bipartition trees and bootstrap trees were input into ASTRAL-II v4.10.2 (Mirarab et al. 2014), which was run with 100 ASTRAL multilocus bootstrap replicates and default settings. ASTRAL generates two types of trees: a local posterior probability (LPP) based on the quartet scores calculated from gene trees, and a multilocus bootstrap (MLBS) that calculates ASTRAL trees from RAxML bootstrap pseudoreplicates generated for each gene, summarized with a greedy consensus tree (Sayyari and Mirarab 2016). To assess whether low support for some areas of the topology was related to poor signal within gene trees, we used SVDquartets (Chifman and Kubatko 2014), an alternative method for constructing species trees that are consistent with the presence of deep coalescence. SVDquartets takes a site-by-site approach; as such, the analysis was conducted with the concatenated supermatrix using the exhaustive method and 100 bootstrap replicates in PAUP^* v4.0a168 (Swofford 2003).

Structure Analysis—For SNP-discovery from the HybSeq data, raw reads were aligned against the assembled exons and flanking non-coding regions from O. harringtonii individual 44_CO using the alignment algorithm BWA-MEM (Li 2013). This alignment protocol assumed that each locus was independent. The MergeBamAlignment and MarkDuplicates tools in PICARD ( http://broadinstitute.github.io/picard) were used to adjust alignments and remove PCR duplicates, respectively. Exonic regions of each sample were then genotyped individually using the GATK HaplotypeCaller tool in cohort mode (GVCF; McKenna et al. 2010). All 50 individuals were considered jointly to call genotypes using GenotypeGVCFs in GATK (McKenna et al. 2010). Resulting variants were hard filtered using the GATK VariantFiltration tool following all GATK best practices manual suggestions except for Quality by Depth (QD), which was set to be greater than 5.0 (DePristo et al. 2011; Van der Auwera et al. 2013). These data were prepared for Structure using PLINK to remove all SNPs with more than 10% missing data (Chang et al. 2015), resulting in a total of 11,226 SNPs.

SNP data were analyzed with ParallelStructure v. 2.3.4 (Besnier and Glover 2013) on XEDE using the CIPRES Science Gateway v. 3.3 (Miller et al. 2010). Structure was run for K = 1–9 with the admixture model. For each value of K, Structure was run five separate times with a 50,000 burn-in and 100,000 iterations. Using Structure Harvester (Earl and vonHoldt 2012), results were examined across values of K, including the approach of Evanno et al. (2005). The five replicates for each K value were consolidated using CLUMPP to create Structure plots (Jakobsson and Rosenberg 2007).

Results

Concatenated Phylogenetic Reconstruction—Maximum likelihood analysis of concatenated data recovered a single species tree with 85% of nodes having bootstrap values ≥ 75 (Fig. 3). This topology supports the monophyly of three of the four narrow endemic species (sampling only one individual of O. brandegeei restricts conclusions as to species monophyly), whereas the fifth species, O. cespitosa, is paraphyletic. There is strong support for a sister relationship between O. brandegeei and O. cavernae (BS = 100). Together these are sister to all other taxa in O. sect. Pachylophus, within which there are several well-supported lineages. One of these lineages includes all O. cespitosa subsp. crinita individuals as well as two morphologically and ecologically intermediate O. cespitosa subsp. marginata individuals (17_NV and 24_CA), which have affinities for the geographically close O. cespitosa subsp. crinita individuals 20_NV and 19_CA, respectively. Sister to this lineage is a clade containing a monophyletic O. psammophila (BS = 100), which is nested within O. cespitosa subsp. cespitosa. Additionally, this clade contains one individual of O. cespitosa subsp. navajoensis (15_UT) that has a strong affinity for O. cespitosa subsp. cespitosa individual 35_ID. A third lineage contains a monophyletic O. harringtonii (BS = 100) that is strongly supported as sister to nine individuals of O. cespitosa subsp. marginata (BS = 100). Oenothera harringtonii + O. cespitosa subsp. marginata are sister to the single accession of O. cespitosa subsp. macroglottis (BS = 100), and sister to this group is a clade of O. cespitosa subsp. navajoensis (BS = 100).

Coalescent-Based Phylogenetic Reconstruction—The topologies inferred from the coalescent-based analyses are largely concordant with that inferred by concatenation, but with less resolution of relationships [ASTRAL (AST) topology, Fig. 4; topology inferred from SVDquartets (SVD) is not shown]. Both coalescent-based analyses strongly support the monophyly of three of the narrow endemic species (only one individual of O. brandegeei was sampled). Oenothera brandegeei and O. cavernae are supported as sister (AST MLBS = 100, SVD BS = 99). Within the rest of O. sect. Pachylophus, a strongly supported lineage (AST MLBS = 98, SVD BS = 96) includes all O. cespitosa subsp. crinita individuals, as well as two O. cespitosa subsp. marginata individuals (17_NV, 24_CA). As in the concatenated topology, these two individuals have affinities for O. cespitosa subsp. crinita 20_NV and 19_CA, respectively. This lineage is sister to a clade (AST MLBS = 89, SVD BS = 95) that includes a monophyletic O. psammophila (AST MLBS = 100, SVD BS = 100) nested within O. cespitosa subsp. cespitosa. These two lineages are sister to an individual of O. cespitosa subsp. navajoensis (15_UT). As in the concatenated analysis, a third lineage contains a monophyletic O. harringtonii (AST MLBS = 100, SVD BS = 100) that is strongly supported as sister to O. cespitosa subsp. marginata (AST MLBS = 100, SVD BS = 100). Oenothera harringtonii + O. cespitosa subsp. marginata are sister to O. cespitosa subsp. macroglottis (AST MLBS = 100, SVD = 100). These taxa are together weakly supported (AST MLBS = 81, SVD BS < 60) as sister to a clade of O. cespitosa subsp. navajoensis (AST MLBS = 100, SVD BS = 91).

Structure Analysis—Results using the Evanno et al. (2005) approach yielded two ΔK peaks at K = 2 and K = 4, with no improvement for K = 5 through K = 9; thus, we show K = 2 through K = 4 (Fig. 5). Across K = 2 through K = 4, similar genetic clusters correspond to Oenothera sect. Pachylophus taxa (Fig. 5). The main difference is that K = 3 separates O. cavernae and O. brandegeei into their own distinct genetic cluster. By contrast, K = 4 adds an additional genetic cluster (yellow) that may reflect admixture, but otherwise has no clear biological relevance, as this cluster occurs in all taxa except O. cavernae. Thus, comparison of K = 2 through K = 4 suggests that the samples are best explained by three genetic clusters (Fig. 5B). The three genetic clusters align with taxonomic designations, with taxa having a majority assignment to a single cluster. One genetic cluster is almost exclusively associated with O. cavernae and O. brandegeei, with all O. cavernae individuals showing high membership (65%) (gray; Fig. 5B). Oenothera brandegeei shows 39% membership in this same cluster, with 58% membership in a second genetic cluster (orange; Fig. 5B). Oenothera psammophila, O. cespitosa subsp. cespitosa, and O. cespitosa subsp. crinita have high membership (> 88%) in the second genetic cluster (orange; Fig. 5B); the two intermediate individuals of O. cespitosa subsp. marginata (24_CA and 17_NV) also show high membership (> 88%) to the same cluster. The third genetic cluster (blue; Fig. 5B) is associated primarily with O. harringtonii and O. cespitosa subsp. marginata, with over 58% membership for all individuals excluding those two O. cespitosa subsp. marginata (24_CA and 17_NV) mentioned above. All O. cespitosa subsp. navajoensis individuals and O. cespitosa subsp. macroglottis share high membership in the second cluster (> 59% orange; Fig. 5B), with 17–37% membership in the third genetic cluster (blue; Fig. 5B).

Fig. 3.

Concatenated majority rule consensus tree of Oenothera sect. Pachylophus inferred from 307 loci. Branches with bootstrap values < 60 are collapsed, and bootstrap values > 60 are noted. Samples are colored by taxon, and an individual identifier (unique number and US state where the material was collected; see Appendix 1) is in parentheses following the taxon name.

Fig. 4.

Coalescent-based ASTRAL multilocus bootstrap consensus topology of Oenothera sect. Pachylophus inferred from 307 loci. Branches with multilocus bootstrap values < 60 are collapsed. Bootstrap values > 60 are shown; above the branches are the ASTRAL multilocus bootstrap values, and below the branches are the SVDquartets bootstrap values. Samples are colored by taxon, and an individual identifier (unique number and US state where the material was collected; see Appendix 1) is in parentheses following the taxon name.

Fig. 5.

Structure plots showing the genetic structure of taxa across Oenothera sect. Pachylophus. The columns represent the genetic structure of each individual, and the y-axis indicates the percentage of estimated membership in each of the genetic clusters. A. K = 2. B. K = 3. C. K = 4.

Discussion

Utility of Target Gene Enrichment for Fine-scale Phylogenetic Inference—Genome-wide sequence data from in-solution target gene enrichment allowed for the inference of a robust sectional phylogeny, with much-needed resolution of relationships among taxa within Oenothera sect. Pachylophus. However, the weak support at several nodes in the coalescent-based analyses (Fig. 4) suggest that some gene trees are conflicting and poorly supported, which may explain the historical difficulty in resolving relationships within this section. This difficulty is expected at the low taxonomic level being explored, given that hybridization is suspected in regions where subspecies of O. cespitosa come into contact (Wagner et al. 1985), and introgression and incomplete lineage sorting (ILS) can complicate gene trees and blur species and subspecies delimitations (Galtier and Daubin 2008; Degnan and Rosenberg 2009; Smith et al. 2015). Indeed, reconstructing the genetic structure of the taxa sampled here (Fig. 5) suggests that there are many ancestral alleles shared among taxa. However, topologies inferred by both concatenation and coalescent-based analyses suggest similar relationships among taxa, supporting the robustness of these findings (Figs. 3–4). The ability of HybSeq to capture both exonic and intronic regions and distill their phylogenetic signals despite introgression and ILS, speaks to the power of this method to elucidate historically recalcitrant nodes, even within species. In particular, when combined with sampling multiple individuals per species across their ranges, this method allows for simultaneous capture of data that can be used in both phylogenetic and population genetics approaches that are complementary when considering the historical processes that account for discordance among gene trees.

Evolutionary Relationships within Oenothera sect. Pachylophus—The present study confirms the monophyly of Oenothera sect. Pachylophus and three of the narrow endemic species (O. cavernae, O. harringtonii, and O. psammophila; limited sampling of O. brandegeei restricts inference of monophyly), as well as the paraphyly of O. cespitosa. This concurs with previous results based on morphological data (Wagner et al. 1985).

Taxon relationships were largely concordant across the concatenated and coalescent-based topologies, with differences primarily in the degree of support; 85% of clades within the concatenated topology are supported with BS ≥ 75, whereas 59% of clades in the coalescent-based topology are supported at this level by at least one of the support measures. This reduced support may represent rapid speciation and subsequently high levels of ILS within Oenothera sect. Pachylophus, explaining the difficulties reconstructing relationships in previous studies. Although concatenation appears to offer improved resolution, there is debate regarding the legitimacy of concatenation approaches when there is notable gene tree discordance (e.g. Gatesy and Springer 2014; Mirarab et al. 2014; Liu et al. 2015). Kubatko and Degnan (2007) demonstrated that clade support with concatenation can be misleading, especially in cases with high occurrences of ILS, and others have confirmed this finding (Mirarab et al. 2014; Liu et al. 2015). Considerable ILS is expected in Oenothera sect. Pachylophus, especially among individuals of the five Oenothera cespitosa subspecies. The lower clade resolution in the coalescent-based analyses supports this assertion and indicates the potential for artificially inflated support in the concatenated topology. Thus, we have taken a conservative approach, and focused on relationships supported by both topologies in the discussion below (Figs. 3–4).

Oenothera cavernae and O. brandegeei—Oenothera cavernae and O. brandegeei are strongly supported as sister to each other and are together sister to the rest of Oenothera sect. Pachylophus (Figs. 3–4). The sister relationship between O. cavernae and O. brandegeei is further supported by analysis of genetic structure. Life history traits and morphological data reinforce this finding; flowers of O. brandegeei and O. cavernae clearly represent a morphological diminution of a larger out-crossing phenotype. Additionally, O. cavernae and O. brandegeei are the only autogamous taxa in the section and, as such, they share a “selfing-syndrome” (Wagner et al. 1985; Sicard and Lenhard 2011). This suggests a single transition to self-compatibility in the common ancestor of these two taxa, given that self-incompatibility is common in Onagraceae (Raven 1979), and all other taxa in Oenothera sect. Pachylophus are self-incompatible (Wagner et al. 1985). Transitions to self-compatibility are one of the most frequent evolutionary shifts in angiosperms, whereas the origin of self-incompatibility is rare (Busch 2005; Busch and Schoen 2008; Igic et al. 2008; Evans et al. 2011). Oenothera cavernae and O. brandegeei are annuals that thrive in extremely xeric habitats, some receiving as little as 100 mm of unpredictable rain each year (Wagner et al. 1985). Short reproductive life spans and harsh environmental conditions are hypothesized to limit reproductive assurance and drive selection for self-compatibility (Stebbins 1950; Baker 1955; Morgan et al. 1997; Morgan 2001; Evans et al. 2011; Busch and Delph 2012). This transition from self-incompatibility to self-compatibility may have been the mechanism driving isolation and divergence of the common ancestor of O. cavernae + O. brandegeei from the remainder of the section. Autogamy would greatly limit gene flow, encouraging rapid divergence and resulting in the distinct genetic structure of these taxa (Baker 1955).

The sister taxon relationship between the morphologically distinct (Wagner et al. 1985) O. cavernae and O. brandegeei may be surprising, given that their extant ranges are small and disjunct, isolated by more than 900 km. However, it is likely that they are paleoendemic species that arose following the range constriction of an ancestral taxon with a previously larger pre-glaciation range (similar to that of Hesperoyucca whipplei; see Segraves and Pellmyr 2001) and became isolated from one another due to specialization to different edaphic environments. Interestingly, the chromosomes of the derived annual O. brandegeei differ by at least 3–4 translocations as compared with those of O. cespitosa (Wagner et al. 1985). It remains unknown whether O. cavernae has a similar chromosomal rearrangement, but, if found, this would further support the close relationship between O. brandegeei and O. cavernae. Expanded population sampling would allow for a more robust assessment of the relative roles of self-compatibility, annual habit, edaphic specialization, and chromosomal rearrangements in the evolutionary history of these species.

Oenothera cespitosa subsp. marginata, O. harringtonii, O. cespitosa subsp. macroglottis, and O. navajoensis—The concatenated and coalescent-based analyses both support O. cespitosa subsp. navajoensis as sister to O. harringtonii + O. cespitosa subsp. marginata + O. cespitosa subsp. macroglottis (Figs. 3–4), which concurs with relationships inferred previously from morphological data (Wagner et al. 1985). However, one morphologically intermediate O. cespitosa subsp. navajoensis individual (15_UT) had affinities for other lineages within O. sect. Pachylophus (Figs. 3–4). By contrast, the genetic structure of O. cespitosa subsp. navajoensis, including individual 15_UT, is very similar to that of O. cespitosa subsp. macroglottis and corroborates the results from the concatenated and coalescent-based analyses in suggesting an affinity of O. cespitosa subsp. navajoensis with O. harringtonii + O. cespitosa subsp. marginata + O. cespitosa subsp. macroglottis (Fig. 5B).

Within this lineage, O. cespitosa subsp. macroglottis is sister to O. harringtonii + O. cespitosa subsp. marginata (Figs. 3–4), and this relationship is further supported by analysis of genetic structure. These taxa are distinguished from the remainder of Oenothera sect. Pachylophus based on high membership (> 58%; somewhat lower membership for the single O. cespitosa subsp. macroglottis individual) in the same genetic cluster (blue; Fig. 5B). The close relationship between O. cespitosa subsp. macroglottis and O. cespitosa subsp. marginata is consistent with previously suggested morphological affinities based on leaf length (15–30 cm), floral tube length (4.5–14 cm), fruit capsule shape (lance-cylindric to cylindric), and seed coat characteristics (Wagner et al. 1985). However, the evolutionary affinities of these three taxa conflict with an earlier suggestion that O. harringtonii was most closely related to O. cespitosa subsp. navajoensis (Wagner et al. 1985).

The newly hypothesized relationship among O. cespitosa subsp. macroglottis, O. cespitosa subsp. marginata, and O. harringtonii offers insight into the evolution of this narrow endemic species. Both O. cespitosa subsp. macroglottis and O. harringtonii have relatively limited ranges, with O. cespitosa subsp. macroglottis restricted to the western half of Colorado, and O. harringtonii limited to grasslands in eastern Colorado (Fig. 2). In the Colorado Springs area, the two species occur in close geographic proximity, but are ecologically distinct, with O. cespitosa subsp. macroglottis growing on igneous substrates at higher elevations, and O. harringtonii growing primarily in the silty clay soils of the arid shortgrass prairies at lower elevations (Wagner et al. 1985). As such, it may be somewhat surprising that O. cespitosa subsp. macroglottis and O. harringtonii are not each other's closest relatives. However, given the limited sampling (one individual) of O. cespitosa subsp. macroglottis in the present study, additional individuals are needed to confirm this finding. In contrast to the fairly limited geographic ranges of O. cespitosa subsp. macroglottis and O. harringtonii, O. cespitosa subsp. marginata has a wide geographic range, with the Rocky Mountains forming the eastern limit in Colorado (Fig. 2). Oenothera cespitosa subsp. marginata occurs on a diversity of substrates, although populations around the Colorado Plateau region and within and to the west of the Rocky Mountains are restricted to sandstone, unlike the silty clay soils east of the Rocky Mountains where O. harringtonii is found (Wagner et al. 1985). These taxa have non-overlapping, highly asymmetrical ranges, with the range of O. harringtonii being only 3% that of O. cespitosa subsp. marginata. Adjacent, asymmetrical ranges are typical of evolution by budding (e.g. Anacker and Strauss 2014; Grossenbacher et al. 2014), with habitat and edaphic shifts common in instances of budding speciation among angiosperms. Additionally, budding speciation appears to have been a common evolutionary mode in monkeyflowers (former genus Mimulus; currently Erythranthe and Diplacus), with asymmetrical ranges and “distinct” niches between sister species (Grossenbacher et al. 2014). The edaphic specialization of O. harringtonii to the unique silty clay foothills of the Colorado Front Range may have been the driver of isolation leading to divergence from an ancestral species, following genetic isolation-by-ecology (IBE; Shafer and Wolf 2013).

Oenothera cespitosa subsp. crinita and Introgression with O. cespitosa subsp. marginata—All individuals of O. cespitosa subsp. crinita comprise a single lineage that is sister to O. cespitosa subsp. cespitosa + O. psammophila. This placement of O. cespitosa subsp. crinita differs from previous analysis based on morphological data, which placed O. cespitosa subsp. crinita as most closely related to O. cavernae and O. brandegeei (Wagner et al. 1985). Interestingly, the present study indicates that two individuals of O. cespitosa subsp. marginata (17_NV, 24_CA) have affinities for O. cespitosa subsp. crinita (Figs. 3–4). These observed relationships are likely a result of introgression between the two subspecies, which have been found to hybridize when artificially crossed (Wagner et al. 1985). Additionally, natural hybridization based on observations of morphological intermediates in areas of range overlap has been suspected between subspecies of O. cespitosa and supported by similarities with known artificial hybrid phenotypes (Wagner et al. 1985). Close investigation of the specimen vouchers for both O. cespitosa subsp. marginata individuals 17_NV and 24_CA indicate that they are morphologically intermediate between O. cespitosa subsp. marginata and O. cespitosa subsp. crinita, suggesting introgression. This is supported by geography, as both O. cespitosa subsp. marginata individuals (17_NV, 24_CA) were collected in areas where the subspecies occur in close proximity at the edge of the Great Basin (Appendix 1). Indeed, O. cespitosa subsp. crinita has been previously collected only 0.3 km from O. cespitosa subsp. marginata individual 17_NV. The structure analysis further supports an introgressive origin for these two individuals, as they share high membership in the same genetic cluster (orange; Fig. 5B) as O. cespitosa subsp. crinita.

Oenothera cespitosa subsp. cespitosa and O. psammophila—Oenothera cespitosa subsp. cespitosa + O. psammophila + O. cespitosa subsp. crinita (Figs. 3–4) share a similar genetic structure, with > 88% membership in a single genetic cluster (orange; Fig. 5B). Oenothera cespitosa subsp. cespitosa is paraphyletic, with the narrow endemic O. psammophila nested within it (Figs. 3–4). Vegetative traits, such as the occurrence of glabrous leaves and stems, support the phylogenetic affinities between these two taxa (Wagner et al. 1985). Furthermore, O. psammophila is restricted to a sand dune region located entirely within the range of O. cespitosa subsp. cespitosa, and the ranges are highly asymmetric, with the range of O. psammophila being less than one one-hundredth of the size of the range of O. cespitosa subsp. cespitosa (Lee 1983; see also Fig. 2). As with O. harringtonii, this geographic pattern is strongly indicative of budding speciation (e.g. Anacker and Strauss 2014; Grossenbacher et al. 2014) due to genetic IBE (Shafer and Wolf 2013), which appears to be fairly common in directing gene flow (Sexton et al. 2014). Indeed, the pattern with O. psammophila is very similar to that between the serpentine endemic Layia discoidea (Asteraceae) and its progenitor L. glandulosa, where L. discoidea is phylogenetically nested within L. glandulosa (Baldwin 2005). Oenothera psammophila is likely a neoendemic species that diverged from O. cespitosa subsp. cespitosa as it specialized to a unique sand dune environment. The paraphyly of O. cespitosa subsp. cespitosa suggests a “progenitor-recent derivative” (P-D) relationship between these taxa, as described by Crawford and Smith (1982). As such, O. psammophila diverged from O. cespitosa subsp. cespitosa and accumulated many morphological changes and some genetic differences (Fig. 5B), while genetic and morphological characteristics of the parent species remained relatively static. Relatively few instances of P-D speciation have been reported in the literature, but include several in the genus Clarkia (Onagraceae) (Gottlieb 2003; Lopez et al. 2012). This type of relationship offers a unique opportunity to study speciation, because the direction of evolution and ancestral character states are clear. Thus, further exploration of character shifts between O. cespitosa subsp. cespitosa and O. psammophila could offer valuable insights into this unique type of speciation and the process of evolution.

In sum, even at the low taxonomic level being explored in this study, many historically recalcitrant nodes within Oenothera sect. Pachylophus were successfully resolved; target gene enrichment using both exonic and intronic data has proved a valuable tool to understand evolutionary relationships. We provide a robust phylogeny that elucidates species and intraspecific relationships, and facilitates the inference of drivers of evolution based on the integration of phylogenetic data with geographic and ecological information. Results of this study provide several key findings that extend our knowledge of Oenothera sect. Pachylophus compared to previous work (Wagner et al. 1985). Rather than being independently derived from O. cespitosa subsp. crinita (Wagner et al. 1985), O. brandegeei + O. cavernae are sister to the remainder of Oenothera sect. Pachylophus and likely represent relictual narrow endemics derived from a widespread common ancestor. Additional sampling of the rare O. brandegeei and analysis of chromosomal arrangements in O. cavernae are needed to better understand relationships between these two species. Although previously considered a hybrid between O. cespitosa subsp. cespitosa and O. cespitosa subsp. marginata (Wagner et al. 1985), O. psammophila appears derived as a peripheral isolate via an edaphic shift to a dune habitat from O. cespitosa subsp. cespitosa. Lastly, although Wagner et al. (1985) was equivocal as to whether O. harringtonii was most closely related to O. cespitosa subsp. navajoensis or O. cespitosa subsp. marginata, we confirm that O. harringtonii and O. cespitosa subsp. marginata are sister, with greater sampling needed within the closely related O. cespitosa subsp. macroglottis. Both O. harringtonii and O. psammophila are likely recent species that formed by budding of a more widespread sister species, providing further evidence for this important speciation mechanism in the western North American flora. Given that O. harringtonii and O. psammophila originated from O. cespitosa through budding, O. cespitosa is necessarily paraphyletic, accurately reflecting the evolutionary history of species within this group. Within O. cespitosa, reconsideration of current subspecies circumscriptions or elevation to species rank may be warranted, but requires further sampling.