Samples for all Echinacea species—E. angustifolia DC., E. atrorubens (Nutt.) Nutt., E. laevigata, E. pallida (Nutt.) Nutt., E. paradoxa Britton, E. purpurea, E. sanguinea Nutt., E. simulata McGregor, and E. tennesseensis—were included in our phylogenomic analyses ( Figures 2 , 3 ) along with available sequences for three closely related species in the Heliantheae Cassini tribe—Helianthus annuus L., H. argophyllus Torr. & A.Gray, and Parthenium argentatum A.Gray (Mandel et al., 2019), serving as outgroup taxa for phylogenomic analyses. Sequence data for the three outgroup samples were extracted from NCBI’s Sequence Read Archive (SRA) database ( Appendix 1 ).

DNA was extracted from silica-dried, snap-frozen, and herbarium leaf tissue samples ( Appendix 1 ) using QIAGEN Plant Pro Kits, and a Qubit 2.0 fluorometer broad-range assay was used to assess DNA concentrations. Shotgun sequencing libraries were constructed using Roche KAPA HyperPlus Library Kits with universal Y-yoke stub adapters (30 Mm) and dual-indexed iTru primers (Glenn et al., 2019). Double-stranded DNA was fragmented, aiming for an average fragment size between 500 and 600 bp for sequencing preparation. The DNA was then end-repaired using an A-tailing reaction for 30 min at 35°C before ligation. Universal stub adapters were then ligated to the A-tailed overhangs in an overnight incubation at 20°C, followed by a post-ligation 0.8X KAPA HyperPure bead cleanup. The cleaned ligation product was then amplified for six PCR cycles with the Roche Company KAPA HiFi HotStart Ready Mix and dual-ended primers (i5, i7), followed by a post-amplification 1x bead cleanup using KAPA HyperPure Beads (Roche, USA, Indianapolis IN) to remove free oligonucleotides and smaller fragments. Library fragment size distributions were measured using Bioanalyzer High Sensitivity DNA Kits (Agilent Technologies), and DNA concentrations (mean nM) were quantified through real-time PCR (qPCR), using KAPA Library Quantification kits and KAPA SYBR Fast qPCR Master Mix (Roche, USA, Indianapolis IN). Concentration estimates of each library were size-corrected using the following formula for the 452bp qPCR Standards:

They were then converted from nM to ng/μL using the following formula (660 g/mol being the approximate weight of a base pair):

Libraries were pooled (10–13 libraries per pool), with a target input of 150–200 ng from each library, and pools were bead-cleaned using 0.8x using KAPA HyperPure Beads (Roche). Pools were eluted from the beads in 100 μL of H₂O and SpeedVac™ (Savant) concentrated to 7 μL as described in the Arbor Biosciences Hybridization Capture protocol. Each pool was incubated at 60°C for 24 h with 5.5 μL of the Compositae-ParaLoss-1272 target capture baits (Moore-Pollard et al., 2023). Following bait capture, hybridization reactions were amplified for 14 PCR cycles to enrich for targets, followed by a 0.8× bead cleanup for purification. Fragment size distributions for target-enriched libraries were assessed using Bioanalyzer High Sensitivity DNA chips (Agilent Technologies) and qPCR. Hybridization pools were combined for a final concentration of 10 mM for sequencing. Sequencing on the Illumina NovaSeq X Plus platform was performed through the SeqCenter (Pittsburgh, PA) sequencing service provider.

All Hyb-Seq reads were quality-trimmed using Fastp v.0.23.2 (Chen et al., 2018). Reads that were shorter than 21 bp after trimming were removed. The HybPiper v.2.1.6 pipeline (Johnson et al., 2016) was used to map the filtered reads to target sequences in the Compositae target file (Mandel et al., 2019) and create “supercontig” assemblies for each target locus for each sample.

Multiple sequence alignments for each locus were constructed for each target locus using MAFFT v.7.505 (Katoh and Standley, 2013) with the “—auto” flag, instructing MAFFT to choose the best alignment strategy to best fit the data. Misaligned sequences were identified and trimmed using trimAl v.1.4.1 (Capella-Gutiérrez et al., 2009). Maximum likelihood gene trees with support values were estimated using IQTree v.1.6.12 (Nguyen et al., 2015) with the “– mfp and -bb “ flags to implement ModelFinder optimization (Kalyaanamoorthy et al., 2017) and the ultrafast bootstrap (Minh et al., 2013), respectively. Species tree estimation was performed using the summary method ASTRAL v.5.7.8 (Mirarab et al., 2014) with the unrooted gene trees as input. The “-t” flag was used to report quartet scores for each node. Quartet frequencies were visualized as pie diagrams using the R packages GGTree, Ape, and TidyVerse.

Using a plastome gene target file, plastome-encoded genes were extracted from the trimmed target-enriched sequence data as bycatch. HybPiper v.2.1.6, MAFFT v.7.505, and Trimal v.1.4.1 were used as described above to capture, assemble, and align plastome sequences. Multiple sequence alignments for 76 plastome-encoded genes were concatenated in Geneious v.2023.2.1 (Kearse et al., 2012) to construct a super matrix. A maximum likelihood tree and ultrafast bootstrap support values were estimated on the concatenated alignment using IQTree v.1.6.12 (Nguyen et al., 2015).

Target capture efficiency ranged between 9% and 24% for the Compositae-Paraloss-1272 bait set, resulting in 1,234–1,272 recovered genes per sample ( Figure 4 ), and gene trees for all targeted loci were estimated for the ASTRAL species tree analysis. The species tree inferred from the analysis of 1,272 nuclear gene trees placed a clade of E. purpurea samples as sister to all other Echinacea species ( Figure 5 ). Most species were found to be monophyletic in the species tree analysis with local posterior probabilities ranging from 0.8 to 1 ( Figure 5 ). At the same time, high quartet frequencies for alternative resolutions of each node (Q2 and Q3) implicated a high degree of gene tree–species tree discordance, suggesting rapid speciation and incomplete lineage sorting, particularly along the spine of the species tree ( Figure 5 ). Alternative quartet frequencies (Q2 and Q3) are generally balanced, as would be expected with ILS ( Figure 5 ). Two species, E. pallida and E. angustifolia, were not recovered as monophyletic. An E. angustifolia individual appeared in a clade with E. paradoxa, while E. pallida samples were scattered among clades containing E. sanguinea, E. atrorubens, E. angustifolia, and E. simulata ( Figure 5 ).

Plastome genes were assembled using off-target reads from the bait capture sequence data. The plastome tree revealed extensive polyphyly for each species ( Figure 6 ), and the branch lengths for internodes along the spine were very short ( Figure 6 ). Most bootstrap support values along the spine of the plastome tree were low despite the massive length of the concatenated plastome gene alignment, including 67,743 bases ( Figure 6 ). Interestingly, samples E_purpurea_CJ18, E_angustifolia_CJ23, and the outgroup species exhibited much longer branch lengths compared to all other ingroup samples.

Here, we analyzed Echinacea species relationships based on an analysis of single-copy target capture genes (Mandel et al., 2019; Moore-Pollard et al., 2023). The large number of recovered genes for this genus indicates that the Compositae-Paraloss-1272 bait set is effective for sunflower family taxa beyond those originally tested (Moore-Pollard et al., 2023).

Previous plastome phylogenomic analysis with single samples per species (Zhang et al., 2017) provided a misleading view of Echinacea relationships and diversification. The branching patterns evident in our plastome ( Figure 6 ) and ASTRAL ( Figure 5 ) trees suggest that the evolutionary history of Echinacea includes an early burst of rapid speciation and some degree of interspecific hybridization. In contrast to the plastome tree presented by Zhang et al. (2017), our species tree ( Figure 5 ) places the cosmopolitan E. purpurea as sister to a clade with the remainder of the genus, and most southeastern species form a grade with a clade dominated by midwestern species arising later in the evolutionary history of Echinacea. A rapid rate of early branching (i.e., speciation) is indicated by the short internodes along the spine of the Echinacea species tree ( Figure 5A ) and extensive signal for incomplete lineage sorting as seen in the high frequency of alternative quartets in the species tree analysis ( Figure 5B ). The lack of monophyly among samples of each species in the plastome tree ( Figure 6 ) may also be a consequence of ILS. At the same time, polyphyly of E. pallida and E. angustifolia in both the ASTRAL ( Figure 5 ) and plastome ( Figure 6 ) trees implicates hybridization and interspecific gene flow.

Nearly even quartet frequencies for many nodes in the species tree estimation, shown in Figure 5B as pie charts, indicate high levels of gene tree–species tree discordance. Given the low posterior probabilities for many nodes on the spine of the clade dominated by midwestern species ( Figure 5A ), we cannot discount the possibility of polytomies in the Echinacea species tree (i.e., ancestral species spawning more than two daughter species). In support of this hypothesis, a number of nodes exhibit quartets with equal frequencies for all three resolutions of nodes on the spine of the Midwest-dominated clade ( Figure 5B ). This pattern would be expected in the face of rapid radiation and rampant ILS. At the same time, the polyphyly of some species implicates post-speciation interspecific gene flow.

Echinacea angustifolia and E. pallida are polyphyletic in the species tree based on analysis of 1272 nuclear genes. While most E. angustifolia samples were grouped together in the species tree ( Figure 5 ), one individual was placed in a clade with E. paradoxa. Echinacea pallida samples exhibited a higher degree of polyphyly, with samples being placed in several clades with E. sanguinea, E. atrorubens, E. angustifolia, and E. simulata, all with ranges in the Midwest. Previous work has documented that all species that have overlapping ranges can hybridize (Kindscher, 2016; McGregor, 1968; McKeown, 1999). We hypothesize that more extensive range overlap among midwestern species has resulted in higher rates of introgressive hybridization. Additional sampling across the ranges of E. simulata and the midwestern species ( Figures 3 , 5 ) is required to gain a deeper understanding of the influence of interspecific gene flow among these species.

Most nodes on the plastome tree have bootstrap robust support values, but as seen in the species tree analysis, some nodes on the spine of the tree are poorly supported, suggesting rapid early divergence among sampled plastome lineages. Whereas Zhang et al. (2017) inferred distinct clades for eastern and western species in their whole-plastome analysis, this interpretation of the plastome history does not hold up when we sampled more than a single individual per species. Whereas Zhang et al. (2017) identified distinct eastern and midwestern species clades, our plastome tree ( Figure 6 ) suggests rapid speciation early in the history of the genus, and multigene species tree estimation ( Figure 5 ) suggests that E. purpurea, with the broadest range of all Echinacea species, is sister to a clade with the remainder of the genus. The lack of species monophyly in the plastome tree and the poor correspondence between the species and the plastome trees ( Figure 5 vs. Figure 6 ) are consistent with the hypothesis that species lineages retained ancestral variation in plastome haplotypes as they were diversifying. Hybridization resulting in interspecific exchange of plastomes may have also contributed to the polyphyly of the species tree ( Figure 5 ) in the plastome phylogeny ( Figure 6 ).