Our phylogenomic analyses included a total of 58 samples, 56 samples from 36 Euphrasia species, and two outgroup species, Bartsia alpina L. and Neobartsia chilensis Uribe-Convers & Tank. This material included a combination of newly sequenced samples and previously generated data, with full sample information provided in Supplementary Table S1. To investigate global evolutionary relationships, broad-scale phylogenetic conflict, and diversity in modes of ploidy across the genus, we selected a shallow sample of taxa that maximized representation of geographic regions (including Northern and Southern Hemisphere taxa), and to capture taxonomic diversity and anticipated evolutionary divergence times. To investigate phylogenetic relationships and species cohesion in postglacial northern European Euphrasia, particularly in Britain, we used available sequences from a range of different studies to maximize species coverage, and where possible to include multiple individuals per species.

For the broad-scale analysis, we sequenced herbarium material from 17 Euphrasia species. Herbarium samples were obtained from the University of Copenhagen (C), the Royal Botanic Garden Edinburgh (E), and Oslo University Herbarium (O). The herbarium samples spanned 1861–2019 and included a broad range of collection localities covering 15 countries including Canada, New Zealand, and Sweden (Supplementary Table S1). DNA was extracted from samples using the Qiagen DNEasy Plant Extraction kit. Extractions were quantified using the Qubit 2.0 Fluorometer (Applied Biosystems).

Library building was performed using Copenhagen University’s EvoGenomics’ in-house BEST protocol (Carøe et al., 2018), which is a PCR-based, short-insert library preparation method designed to maximize historical DNA potential by accounting for low extraction yields. Purification steps were used both at the extraction stage and the library building stage and included both SPRI magnetic beads (Beckman Coulter) and membrane filter MinElute PCR Purification spin columns (Qiagen). Subsequently, Illumina dual indexes (8 bp) were used to facilitate multiplexing of samples. Sequencing was outsourced to NovoGene EU, using the NovaSeq 6,000 with 150 bp paired-end (PE) sequencing.

For the analysis of phylogenetic relationships and monophyly in postglacial northern European Euphrasia, sequencing data for 41 individuals were sourced from three previous studies. First, we integrated data for 18 Euphrasia samples previously used in a population genomic study of British Euphrasia, with a focus on tetraploid species on Fair Isle, Scotland (Becher et al., 2020). This study generated a reference genome of the tetraploid species E. arctica (described below) and high-coverage short-read data for 17 additional Euphrasia samples. All 18 samples had nuclear SNPs called relative to the reference genome, and plastid genomes and nrDNA arrays assembled de novo. Second, low-coverage short-read sequencing data of 12 samples: 10 other British Euphrasia, one Austrian Euphrasia, and an outgroup Bartsia alpina, were available from a study characterizing the landscape of genomic repeats (Becher et al., 2021), with the raw data reanalyzed here and used for de novo assembly of plastid genomes and nrDNA, and mapping to the reference genome. Finally, we included short-read data for 11 previously unpublished samples (Twyford, Unpublished Data) where data was available on the Sequence Read Archive (SRA; SRR17976421 - SRR17976431). These represent 10 diverse Euphrasia taxa and an outgroup Neobartsia chilensis. These low-coverage genome skims were generated from NEB Ultra PCR-based libraries sequenced with 125 bp PE sequencing on the Illumina HiSeq 2500 or 150 bp PE sequencing on the Illumina NovaSeq 6000 at Edinburgh Genomics.

Our final data set included 31 samples collected in Britain, including multiple samples for: E. anglica Pugsley (2 samples), E. arctica Lange ex Rostr. (7), E. confusa Pugsley (2), E. foulaensis Towns. ex Wettst. (5), E. micrantha Rchb. (7) and E. vigursii Davey (2). These also included two putative hybrids (E. confusa x E. foulaensis, E. arctica x foulaensis) and two species of putative hybrid origin (E. rivularis, E. vigursii, Yeo, 1956). We also include a sample of ‘Euphrasia fharaidensis’, a UK endemic awaiting formal description (French et al., 2008). Herbarium material from all newly sequenced samples are lodged at E.

Phylogenetic analyses were performed independently for plastid genomes, nrDNA arrays, conserved nuclear scaffolds and the complete set of sequence reads. For each data set, phylogenies were annotated with available ploidy and geographic information. Ploidy information came from Metherell and Rumsey (2018) and Becher et al. (2021) for British samples and Gussarova et al. (2008) for other taxa. Ploidy information was available for 35 individuals, with missing information for many newly sequenced non-British samples. Geographic areas were annotated on trees following the areas defined by Gussarova et al. (2008).

Previous genomic analyses of British Euphrasia inferred individual ploidy, and whether subgenomes are likely to be shared between individuals, based on sequencing coverage per genome scaffold relative to the E. arctica reference genome (Becher et al., 2020). In this previous analysis, while all samples had similar coverage across the conserved scaffold set (which are likely to be within a conserved subgenome present across diploid and tetraploid British Euphrasia), diploids had no (or very low) coverage in a large number of scaffolds restricted to the tetraploids. Here, we perform a similar mapping depth analysis across our global Euphrasia samples to understand whether genome structure is conserved across diverse species. The per-scaffold coverages were computed for each sample, retaining only scaffolds longer than 1 kb. In order to compare the samples within the study and to previously published data (Becher et al., 2020), the per-scaffold coverages were normalized by the average coverage across these scaffolds. Further, the analysis was restricted to ~10,000 scaffolds previously identified as conserved across the genus. Hierarchical clustering was performed on the resulting matrix of normalized per-scaffold coverages at these ~10,000 scaffolds. Coverages were then visualized per sample and per scaffold in a heatmap following Becher et al. (2020). To aid visualization samples were ordered based on relatedness inferred from hierarchical clustering of pairwise Manhattan distances between the samples’ mapping depth profiles (Becher et al., 2020).

We successfully assembled the plastid genome for 38 Euphrasia individuals and 2 outgroups and compared these to 18 previously assembled plastid genomes. Overall the assembly length was consistent across samples (range: 140,581–145,113 bp for Euphrasia species, up to 153,370 bp in Bartsia alpina), with a mean size of 144,792 bp across Euphrasia species. The average plastid genome GC content was 38.3%. Pairwise identity between Euphrasia samples was high at 99.1% and with 87.1% of sites identical across species, with these rising to 99.8 and 98.8%, respectively, across the 31 British samples. The final alignment of 56 Euphrasia and 2 outgroup plastid genomes was 157,898 bp in length.

Phylogenetic analyses of plastid genomes revealed a deep split between clades broadly corresponding to Northern Hemisphere taxa, and Southern Hemisphere taxa plus samples from Japan (Figure 1). The only exception to this biogeographic split were two samples of the European species E. cuspidata, which were placed on a long branch separate from all other samples, consistent with it being a morphological distinct diploid taxon belonging to a separate taxonomic section. Within clades, there was significant phylogenetic complexity, with some patterns of relatedness representing geography, ploidy, or species identities, though many relationships are hard to explain. Within the poorly sampled clade of largely Southern Hemisphere species (represented by eight samples) species do not cluster by geography, with species sampled from the same country (such as New Zealand) separated on the tree.

Within the Northern Hemisphere clade, there is a lack of discernable overall phylogenetic structure and the tree is characterized by extremely short terminal branches. However, some individual or species-level patterns of plastid haplotype sharing and relatedness emerge. These include clusters corresponding to: three samples of E. micrantha from Fair Isle; related tetraploid E. arctica and E. foulaensis from Fair Isle; four diploid samples from E. anglica and E. vigursii from England; and a cluster of distinct Eurasian species E. tatarica, E. pectinata, and E. regelii. In contrast to these clusters, many more patterns of relatedness appear more complex, for example four other samples of E. micrantha (excluding the three samples that cluster) are largely spread throughout the wider Palearctic and Alps group. This pattern of individuals being scattered throughout the Northern Hemisphere clade is also seen with E. arctica and E. foulaensis.

Overall, the plastid phylogeny highlights the phylogenetic complexity present in the genus, with only weak clustering by geography, and with multiple samples within species not showing clear taxonomic coherence.

The nrDNA array was successfully assembled for all 38 samples except E. gothlandica and E. micrantha sample M4, which consistently failed. The new nrDNA assemblies were aligned with the 18 existing assemblies (Becher et al., 2020) to produce a conserved nrDNA alignment of 5,845 bp in length. Per sample conserved nrDNA array lengths varied between 5,571–5,828 bp across the alignment of 55 Euphrasia individuals and 2 outgroup samples. Pairwise sequence identity between Euphrasia individuals was 99.2, and 91.8% of sites were identical across the alignment. Total GC content was 54.3%. Most sites were unambiguously identified, with only 423 sites in the alignment (0.13% sites) coded as ambiguous. Across samples 21 individuals had no ambiguous sites (36.8%). These samples were either Southern Hemisphere taxa, all of which had no ambiguities (apart from octoploid E. antarctica), known Northern Hemisphere diploids (e.g., E. alpina and E. vigursii), or Northern Hemisphere species of unknown ploidy clustering with the diploids (E. bottnica, E. canadensis, see below). Exceptions to this finding are diploid E. rostkoviana and the putative diploid hybrid species E. rivularis, which had some ambiguous sites. All known tetraploid samples from the Northern Hemisphere had at least two ambiguous sites. Most ambiguous sites were not random in their position across the alignment and instead were common at sites segregating for two alleles, indicative of the retention of multiple nrDNA copies rather than assembly errors.

Phylogenetic analyses of the nrDNA array confirmed the presence of a clear and moderately well-supported biogeographic break (BS = 86%) largely corresponding to Northern vs. Southern Hemisphere taxa (Figure 2). The southern clade is characterized by long branches typical of older and more divergent lineages, but also reflect artifacts related to poorer sampling. In the Northern Hemisphere clustering is largely by ploidy, with clear separation of most diploid taxa on a well-supported (BS = 100) long branch from most known tetraploids. Within tetraploids, there is a clade comprised of Siberian E. tatarica and Chinese E. regelii and E. pectinata, and a spanish sample of E. willkommii. The better-sampled clade of predominantly North Western European tetraploids lacks strong geographic or taxonomic structure, with species and geographic locations being largely intermixed. For example, nrDNA sequences of E. micrantha are scattered across the clade and with one sequence on a long branch. However, there is evidence of some geographic structure and clustering by taxonomy, for example with four samples representing three species and one hybrid found on Fair Isle possessing identical nrDNA sequences.

Our species tree of conserved nuclear scaffolds (Figure 3), and associated quantification of discordance with DiscoVista (Figure 4), revealed that 66% (36/55) of internal nodes had one (and only one) topology across more than 33% of the gene trees. For the remaining (34%) internal branches, there is an equal weighting for a second or third possible topology. Across the species tree, there is notable variation in patterns of discordance. There is generally a lower likelihood of the alternative topology in branches within the clade of Southern Hemisphere taxa and species from Japan (nodes 33–38), which are newly sequenced here, as well as in British diploids (nodes 20–24), or where there is geographically cohesive species sampling of tetraploids (such as three sample of the selfing species E. foulaensis sampled on Fair Isle, node 12). In contrast, there is a near equivalent representation of two, or all three, possible relationships in many other branches across the phylogeny, particularly those involving British tetraploids. While the general trend is of greater discordance in branches connecting recent species relationships in European tetraploids, there are also some early diverging nodes with discordance where an alternative topology is frequent (e.g., nodes 40 and 47), showing complexity across the Euphrasia phylogeny.

This species tree supports the plastid and nrDNA phylogenies in recovering an early diverging clade predominantly composed of Southern Hemisphere taxa, with the South American taxa placed sister to the Australia-South Asian taxa. The better-sampled Northern Hemisphere clade has the Japanese taxon E. maximowiczii included within an otherwise exclusively North American-European clade. Within the Northern Hemisphere clade, there are clearly resolved early diverging relationships between mainland European, Japanese, and North American taxa. British diploid species are resolved as monophyletic, while British tetraploids fall in three large clades, each with a small number of European or North American taxa. Of particular note is one clade where mainland British samples are placed sister to a monophyletic sub-clade of species found on Fair Isle. Euphrasia micrantha is found in all three clades including British taxa, including one clade of three E. micrantha samples with European E. gracilis.

The neighbor-joining tree based on distances estimated by MASH sketches broadly mirrors the topology of the species tree built using trees inferred from the 3,454 conserved sequence scaffolds (Supplementary Figure S1), albeit with some short branches. The only sample showing discordant placement relative to the species tree is Euphrasia gothlandica (E29), which suffers from low coverage.

Our tanglegram analyses revealed extensive phylogenetic discordance between genomic regions. Given the large amount of informative sites and high support, we focus on comparisons to the nuclear species tree analysis. Incongruence is seen in comparisons with the plastid genome phylogeny (Figure 5), as the plastid analysis does not recover any major clades apart from the early diverging predominantly Southern Hemisphere group. Most notably, the plastid phylogeny does not recover a group of diploid taxa or resolve any overall geographic structure within the Northern Hemisphere group (Figure 5). When the nuclear species tree is compared to the nrDNA tree (Figure 6), there are some similarities including groups corresponding to ploidy, though its placement within the phylogeny differs, appearing on a long branch sister to all Northern Hemisphere tetraploids in the nrDNA tree, and being placed in a more derived clade in the nuclear species tree. Interestingly however, both the nuclear species tree and the plastid tree have some consistent individual-level relationships, such as three samples from Fair Isle, whereas this group is not recovered in the nrDNA analysis. Otherwise the plastid and nrDNA phylogenies are largely incongruent except some early diverging Southern Hemisphere lineages (Supplementary Figure S2). One interesting case of incongruence is diploid E. rivularis, a species of putative cross-ploidy hybrid origin. This is the only diploid species clustering with a group of tetraploids in the plastid tree, but yet it clusters with other diploids in the nrDNA tree, consistent with its proposed origins.

We assessed sample ploidy and subgenome relationships based on short-read sequence coverage relative to scaffolds present in the genome of tetraploid E. arctica. While the mean mapping depth of ~14X is generally sufficient to infer presence/absence and estimate copy number, there was notable variation, and 5 samples had below 5X mapping which partly obscures patterns (Supplementary Table S2).

Euphrasia bottnica (E1, Finland) demonstrated a coverage pattern similar to previously analysed British diploids, with a large set of scaffolds having no (or nearly no) reads mapping, with these absent scaffolds corresponding to the divergent tetraploid subgenome (Figure 7). Similarly, E. calida (E11, Iceland), E. cebennensis (E13, France), E. disjuncta (E21, Canada), E. latifolia (E27, Canada), and E. gracilis (E30, Sweden) show broadly similar coverage patterns to British tetraploids, albeit with higher variance. The other samples had distinctly different mapping depth patterns. However, four samples from the Southern Hemisphere: E. townsonii (E5, New Zealand), E. cuneata (E17, New Zealand), E. lasianthera (E22, Australia), and E. erecta (E25, New Guinea), show starkly different coverage patterns for a subset of the conserved scaffolds including both Euphrasia-wide conserved scaffolds (3454) and tetraploid-only conserved scaffolds (~7,000). For this subset of scaffolds, these samples show double the average genome coverage, suggesting that these samples are either octoploid or have partial genome duplication post-polyploidization. The inference of polyploid history for samples E3, E12, and E29 is more ambiguous, though these represent tetraploids divergent from the the reference.

Phylogenetic discordance has been observed in numerous plant studies and is increasingly considered the norm (Rose et al., 2021). Here, we confirm that signals of phylogenetic discordance observed with Sanger sequencing of few loci at the regional (Wang et al., 2018), and the global scale (Gussarova et al., 2008), are indeed detected with genomic data. The lack of clear geographic, taxonomic or ploidy-related structuring in the plastid genome phylogeny, coupled with incongruence to the well-supported nuclear genomic phylogeny, shows the plastid does not track other loci and variation is shaped by other evolutionary forces. In particular, frequent hybridization as well as self-fertilization may lead to the geographically restricted fixation of plastid haplotypes giving a local geographic signal conflicting with species boundaries.

While there is less discordance between nrDNA and the nuclear genome there is still conflict, which may be due to the stochasticity underlying a single gene tree, or a specific outcome of nrDNA being maintained in multiple copies and subsequently experiencing concerted evolution (Xu et al., 2017). Perhaps most importantly, nucleotide ambiguity in the tetraploids suggests the maintenance of multiple nrDNA variants following polyploidy, a finding not observed with direct Sanger sequencing in Euphrasia (Wang et al., 2018). While intra-individual variation is problematic for reconstructing phylogenies, the broad concordance of major clades between nrDNA and the nuclear species tree suggests this issue does not obscure signal present in the data. Future work will look to clarify evolutionary relationships using phased nrDNA sequences as has been done with polyploid taxa in the Asteraceae (Fehrer et al., 2021), in the hope that this will reveal the currently unknown second progenitor of British tetraploid Euphrasia species (Becher et al., 2020; discussed below).

In addition to challenges with the retention of multiple nrDNA copies, there is also ample evidence consistent with incomplete lineage sorting, with individuals with multiple samples combining variation greatly predating speciation, including variants across the tetraploid clade which was previously estimated to be 7.3 million years old (Gussarova et al., 2008). Regardless of the evolutionary processes shaping variation, these results highlight that an “extended barcode” (Coissac et al., 2016) based on plastid and nrDNA from genome skimming will fail to identify species in Euphrasia, as well as a range of other complex groups like willows (Wagner et al., 2021). Here, researchers should instead look to sample the nuclear genome, using methods such as whole-genome resequencing or sequence capture with target enrichment probes (Johnson et al., 2019).

Our analyses used extensive sampling of British species, coupled with representatives from the geographic range and phylogenetic diversity of Euphrasia. This allows us to investigate the colonization and speciation history of recent postglacial species divergence. Previous work has shown that diploid and tetraploid Euphrasia diverged long before recent pleistocene glaciation events and colonized Britain independently (Wang et al., 2018). Our higher resolution nuclear genomic data further show British tetraploids are present in at least three clades mixed with mainland European samples, and therefore likely represent at least three waves of colonization. These results mirror many other plant phylogeographic studies, where different genetic lineages present in Britain have colonized recurrently from continental Europe and co-exist (Wood et al., 2018).

In addition to broad-scale phylogeographic patterns, our work also reveals fine-scale insights, such as showing that the geographically widespread small-flowered selfing taxon Euphrasia micrantha is polyphyletic in analyses of the nuclear genome, plastid and nrDNA array. This finding is surprising, as it is one of the most morphologically distinct Euphrasia species, characterized by purple leaves, stems, and flowers, and has a distinctive ecology being predominantly found in heather moorland (Stone, 2012). Our more limited sampling of a number of other taxa reveals a similar lack of species cohesion. These results are in line with previous genetic and genomic studies showing genetic variation in Euphrasia often clusters by geography rather than by species, at least within a ploidy level (French et al., 2008; Becher et al., 2020). While phenotypic plasticity and taxonomic confusion are profound challenges for studies of Euphrasia (Brown et al., 2020), these are of limited concern at least for E. micrantha, which maintains its morphological distinctiveness under a range of conditions and is unlikely to be confused by Euphrasia experts. This leaves a number of non-mutually exclusive explanations underlying the origin and maintenance of the species.

Firstly, our analyses here either looked at largely non-recombinant single genomic regions (plastid or nrDNA), or aggregated regions with partially independent evolutionary histories (conserved nuclear scaffolds). In particular the use of scaffolds conserved across individuals means we have only investigated genomic relatedness in one subgenome of the tetraploid. As such these results may have overlooked or masked more subtle genomic signatures at individual nuclear loci. It may be that E. micrantha and other Euphrasia species are monophyletic at specific nuclear regions underlying species differences that were not analyzed separately here, but experience homogenizing gene flow, such as via hybridization across the rest of the genome (Twyford and Friedman, 2015). This remains a distinct possibility given weak reproductive barriers between Euphrasia species. Such an explanation would be consistent with either a single origin of the species followed by hybridization, or multiple origins at different sites perhaps from a shared pool of standing genetic variation (i.e., combinatorial speciation sensu Marques et al., 2019). However, the maintenance of such high genetic diversity and divergent haplotypes within species, and within a single sampling location, particularly in a selfing taxon, may also point toward the presence of cryptic species. We are currently pursuing these hypotheses using further genomic sequencing of population samples, where we aim to quantify the extent of hybridization between individuals from different geographic areas, and with contrasting ploidy and mating systems.

The previous study of Becher et al. (2020) identified an allotetraploid origin of British Euphrasia, with one subgenome closely related to British diploids. Our analysis of sequencing coverage revealed that Euphrasia bottnica sampled from Finland possesses a similar genome to extant British diploids, and this species clusters with diploids in the nrDNA phylogeny, suggesting a shared genomic affinity across this region of postglacial recolonization. Similarly, a number of European tetraploid taxa, such as E. calida and E. cebennensis, as well as Canadian E. disjuncta and E. latifolia, have similar sequence coverage patterns to British tetraploids, consistent with a shared allopolyploid origin for these Northern Hemisphere taxa. Perhaps more notable are cases such as E. cuneata and E. townsonii (New Zealand), E. lasianthera (Australia), E. erecta (New Guinea), and E. cuspidata (Austria), which are characterized by an extremely different coverage patterns, with some scaffolds in the E. arctica reference having no coverage and others having two-fold coverage. Firm conclusions of the ploidy and subgenome constituents of these taxa are hard to make given they are likely to show substantial divergence from the British E. arctica reference genome, and are also characterized by low mapping coverage. We have attempted to further characterize sequence reads from these individuals using k-mer based approaches including KAT and Tetmer (Mapleson et al., 2017; Becher et al., 2020), but failed to retrieve a clear signal (Unpublished Results), likely due to low sequencing coverage and potential DNA degradation from these herbarium specimens. Regardless, the finding that these divergent species possess a different genome structure to other Euphrasia warrants further study, and suggests recurrent polyploidy in the genus, in line with known tetraploids and hexaploids being scattered across the Euphrasia phylogeny.

Taxonomically complex groups frequently pose the joint challenges of taxonomic issues, where species definitions may be uncertain and monographic work is often sorely needed, and systematic/phylogenetic issues, where molecular phylogenies show complex patterns of relatedness. Both issues are relevant to Euphrasia. Much has been learnt about species limits since the world monograph of the genus by von Wettstein (1896) and the European revision by Yeo (1978), not least the extent of phenotypic plasticity that taxa exhibit in response to host species and ecological conditions (Karlsson, 1984; Zopfi, 1997; Brown et al., 2020). Our general view is that Euphrasia species as currently described, particularly in Britain, may have been too finely divided, and future monographic work may look to “lump” a range of species where trait differences are minimal or prove unreliable. Despite these taxonomic issues, we note that even the most distinct species, such as E. micrantha, are not monophyletic (discussed above), showing phylogenetic complexity will persist even following taxonomic realignment of species. This is unsurprising given the nature of these species, showing recent speciation, rampant hybridization, and selfing or mixed-mating systems.

One source of samples that has proved particularly useful in our study has been verified material present in herbarium collections. The search for genomic tools that reliably recover information from herbarium specimens is driven by the incredible amount of historic plant diversity contained within these collections (Särkinen et al., 2012; Buerki and Baker, 2016; Brewer et al., 2019). Accessing herbaria’s genomic data will allow researchers to (figuratively) travel through time and space to study extinct taxa and changes in genetic diversity over time. Many studies have now demonstrated the efficacy of genome skimming or target capture to recover genomic data from historical samples (Zeng et al., 2018). Both these approaches rely on a form of “enrichment,” with genome skimming analyzing “naturally enriched” regions (i.e., those at high copy number) while target capture enriches regions homologous to target baits. Here, we show that non-enriched, direct whole-genome sequencing can be successfully used for degraded herbarium material. As well as being used to assemble the plastid and nrDNA array, we were able to map sufficient reads to the E. arctica reference genome to infer sample ploidy. While useful, however, there were notable issues with these analyses, particularly due to low sample mapping depths. This may either be a consequence of species divergence or contamination, with previous work showing over 70% of sequence reads from herbarium material may be contaminants (Bieker et al., 2020). This issue, combined with DNA error profiles of dried plant tissue, prevented us performing further characterization of these genomes, and suggests future work must oversample herbarium DNA to ensure sufficient data post bioinformatic filtering, or use silica dried plant tissue where available (Brewer et al., 2019). Despite these concerns, we found neighbor-joining trees generated from raw sequence reads from herbarium samples, inferred using MASH, mirrored the topology of our more rigorous scaffold-based nuclear phylogenetic analyses, suggesting sample degradation and contamination do not obscure the main signal of genome-wide relatedness.

Phylogenomic analyses of taxonomically complex groups are often made difficult due to reticulation coupled with polyploidy. Here, we circumvented a number of issues by analyzing haploid plastid genomes, though our hope that nrDNA would have been homogenized within an individual appears not to be the case. We similarly focused our nuclear genome analyses on conserved disomically inherited scaffolds, allowing us to compare across diverse ploidies and to represent evolutionary relationships using a species tree analysis. We chose not to further interrogate genomic relationships within putative subgenomes due to the uncertain homology across these diverse species of differing ploidy and with potentially different parental progenitors. Future work in Euphrasia, and other taxonomically complex groups, may look to long-read sequencing and pangenome analyses to better represent structural genomic variation across diverse taxa without reference bias, and to provide robust sorting of homoeologs between subgenomes (Bayer et al., 2021). More integrated polyploidy-aware phylogenomic networks, such as alloPPnet, are also likely to prove fruitful, particularly in the future if this or other methods are developed that are less computationally demanding and allow larger multi-sample data sets as well as more diverse ploidy levels (Rothfels, 2021).

Studies of the extent of discordance in phylogenies have given important insights into a range of topics, including hybridization (Patterson et al., 2012; Martin et al., 2013) and hybrid speciation (Köhler et al., 2021), evolutionary conflict (Hedtke and Hillis, 2010), horizontal gene transfer (Davis et al., 2005) and rates of phenotypic innovation (Parins-Fukuchi et al., 2021). Our study shows taxonomically complex Euphrasia represent a genus where phylogenetic discordance is extensive, at both shallow and deep nodes in the phylogeny. This discordance is likely to be driven by the interaction of different processes, including recurrent rounds of polyploidy, rampant hybridization, and recent postglacial species divergence. Future work will look to estimate the contribution of these processes to phylogenetic conflict in chromosome level genome assemblies.