Phylogenomics of Australian sundews (Drosera: Droseraceae)

Taxon sampling

We included sequence data from 92 Drosera samples (Table 1) in collections held at the State Herbarium of South Australia (AD), the Queensland Herbarium (BRI) and the Beadle Herbarium (NE), supplemented with publicly available sequence data at GenBank (see https://www.ncbi.nlm.nih.gov/genbank/) and the Kew Tree of Life (see https://treeoflife.kew.org/tree-of-life) (Baker et al. 2022; Zuntini et al. 2024). Taxon sampling included representatives of all subgenera and sections occurring in Australia and selected taxa from other continents, following the classification system of Fleischmann et al. (2017a) and Mohn et al. (2023). These included type species acting as representatives for subgenus Arcturia (2 spp.), subgenus Drosera (~110 spp.), subgenus Ergaleium (~150 spp.) and subgenus Regiae Seine & Barthlott (1 sp.). Type species of the recently synonymised (Fleischmann et al. 2017a) sections Annuerecta Lowrie, Ergaleium, Erythrorhiza (Planch.) Diels, Stolonifera DeBuhr, Macrantha Lowrie and Luniferae Lowrie, as had been circumscribed by Lowrie (2013a, 2013b, 2013c), were included.

To fill sampling gaps, herbarium vouchers were made from cultivated D. regia Stephens and D. hamiltonii C.R.P.Andrews held in living collections at The University of Adelaide. Dionaea muscipula (SRR2807634), Aldrovanda vesiculosa (ERR7599885) and Nepenthes madagascariensis Poir. (ERR3672062) nuclear sequence data from the Kew Tree of Life (Baker et al. 2022; Zuntini et al. 2024) were included in our nuclear analyses as outgroups. Plastid sequence data for D. muscipula (NC_035417) and A. vesiculosa (NC_035416) (Nevill et al. 2019) from GenBank were included in our plastid analyses as outgroups.

All described taxa within Drosera subgenus Regiae (1 sp.), D. subgenus Drosera sections Prolifera (4 spp.), Stelogyne (1 sp.), D. subgenus Ergaleium sections Phycopsis Planch. (1 sp.) and Coelophylla Planch. (1 sp.) were sampled. Other species groups, i.e. D. subgenus Arcturia, D. subgenus Drosera sections Arachnopus and Psychophila Planch., D. subgenus Ergaleium sections Bryastrum, Ergaleium and Lasiocephala were not sampled exhaustively. Taxon sampling for this study represents ~50% of the species diversity in Australia; for a comprehensive sample list, see Table 1.

DNA extraction, library preparation and sequencing

Material was removed from herbarium specimens using forceps sterilised by immersing in a 6% sodium hypochlorite solution, followed by a RO water and 70% ethanol rinse. The forceps were physically dried with a Kimwipe (Kimtech Science) before progressing to the next sample. DNA extraction, library preparation and sequencing were performed at the Australian Genome Research Facility (AGRF, Melbourne, Vic., Australia) as part of the GAP initiative supported by Bioplatforms Australia (Sydney, NSW, Australia). Dried plant tissue (~20–30 mg) was ground using a TissueLyser II (Qiagen) with tungsten carbide beads for simultaneous disruption and homogenisation of the sample, as per the manufacturer’s instructions. Genomic DNA was extracted using the DNeasy Plant mini kit (Qiagen) as per the manufacturer’s instructions on a QIAcube Connect (Qiagen). DNA quality was assessed using 1% E-gel with Sybr Safe dye (ThermoFisher) to visualise fragment size range and concentrations were measured with a Quantifluor dsDNA assay (Promega). Libraries were prepared using the NEBNext Ultra II FS Library Prep Kit (New England Biolabs, Ipswich, MA, USA), following the manufacturer’s instructions with inserts on average having a length of ~350 bp. As part of this workflow DNA was enzymatically fragmented to the required size. Libraries prepared by AGRF were enriched using both the Angiosperms353 (A353; Johnson et al. 2018) and OzBaits (Waycott et al. 2021) probe kits. Pooled libraries (12–16 plex) were enriched using the A353 probe kit by hybridising at 65°C with the Arbor Biosciences MyBaits Expert Plant A353 v1 bait set with v5 chemistry (catalogue number 308108.v5) and Arbor Biosciences MyBaits custom OzBaits_NR set (catalogue number 300496R.V5) set with V5 chemistry and 64°C hybridisation following the manufacturer’s instructions.

Remaining genomic DNA was sent to the Advanced DNA, Identification and Forensic Facility at The University of Adelaide and additional libraries were prepared based on similar conditions to the GAP libraries but following the two-step library preparation protocol from Waycott et al. (2021). This was done to enrich ~60 chloroplast loci using the OzBaits_CP bait set. Libraries were pooled 16 plex (6 pools) and enriched with the OzBaits_CP v2.0 (Waycott et al. 2021) chloroplast probe kit (catalogue number 300196R.v5) at 65°C using V5 chemistry. Post-capture libraries were amplified with indexed primers and purified with 1× AMPure XP beads (Beckman Coulter).

To obtain additional plastid and ribosomal sequence coverage, a subsample of the 6 un-enriched pooled libraries was used for a genome skim (i.e. Zeng et al. 2018); 1 µL of library was amplified (8 cycles, 65°C annealing temperature) in a 25-µL reaction with indexed primers using the KAPA HiFi Hotstart ReadyMix (Roche). All libraries were purified with AMPure XP beads (Beckman Coulter) and quantified on a 4150 TapeStation using a High Sensitivity DNA screen tape (Agilent). Both plastid and genome skim libraries were pooled to equimolar concentrations and size-selected (300–700 bp) on a Pippin Prep (Sage Sciences) 1.5% agarose 250–1500-bp gel cassette. The resulting cleaned up product was subsequently sequenced on a NovaSeq. 6000 (Illumina Inc., San Diego, CA, USA) with v1.5 chemistry and 150-bp paired-end reads.

Data processing

Paired-end Illumina reads (150 bp in length) were imported into CLC Genomics Workbench (ver. 20.0.2, see https://digitalinsights.qiagen.com/products-overview/discovery-insights-portfolio/analysis-and-visualization/qiagen-clc-genomics-workbench/) for demultiplexing and adapter trimming. Reads with a Phred score (Ewing et al. 1998) >13 (corresponding to a quality score cutoff of 0.05) were retained. To reduce computational processing time, individuals were randomly sampled to 2.5 million paired-end reads. Sampled reads were then de novo assembled with MEGAHIT (ver. 1.2.9, see https://github.com/voutcn/megahit; Li et al. 2015), using the ‘assemble’ command in CAPTUS bioinformatics pipeline (ver. 1.0.0, see https://github.com/edgardomortiz/Captus; Raza et al. 2023) with settings left as default with the recommended ‘–min_count 3’ and ‘–prune_level 3’ arguments.

We used the CAPTUS ‘extract’ function to extract target gene regions using custom nuclear and plastid references. A set of reference sequences of OzBaits gene targets was derived from sequences of 43 genera within the Caryophyllales that were sourced from the 1KP data set (Carpenter et al. 2019), Genbank (Drosera capensis genome: GCA_001925005; Butts et al. 2016) and the transcript sequences of Aldrovanda vesiculosa and Dionaea muscipula available at https://www.biozentrum.uni-wuerzburg.de/carnivorom/resources (Palfalvi et al. 2020). References for the A353 gene targets were sourced from the Kew Tree of Life Explorer (Baker et al. 2022), consisting of Aldrovanda vesiculosa (ERR3639225), Dionaea muscipula (SRR2807634), Drosera dichrosepala Turcz. (ERR3672058) and Nepenthes madagascariensis (ERR3672062). Plastid references were sourced from GenBank and included coding regions extracted from the plastid genome sequences of Aldrovanda vesiculosa (MK397911, NC_035416), Dionaea muscipula (MK397918, NC_035417), Drosera erythrorhiza (KY651214, NC_035241), Drosera indica L. (MK397919), Drosera regia (KY679199, NC_035415) and Drosera rotundifolia (KU168830, NC_029770). As parameters for the CAPTUS ‘extract’ function, we used a minimum score of ‘0.2’ and minimum identity of ‘80’ for the nuclear data, and a minimum score of ‘0.3’ and minimum identity of ‘90’ for the plastid data.

The CAPTUS ‘align’ command was used to generate multiple sequence alignments for the extracted nuclear and plastid loci using the ‘genes’ format. Contigs were aligned with MAFFT (ver. 7.520, see https://mafft.cbrc.jp/alignment/software/; Katoh and Standley 2013) with the recommended ‘mafft_auto’ argument and the following paralog filtering settings: ‘–max_paralogs 3’, ‘–filter_method informed’ and ‘–tolerance 1’ (refer to the CAPTUS user manual for a full description of the settings for details, see https://edgardomortiz.github.io/captus.docs/). A maximum of three paralogs was included in the pre-filtered alignments based upon initial runs that indicated values above three led to the retention of short and highly divergent copies indicative of contamination. Alignments were trimmed in CAPTUS using the default ‘gappy’ trimming mode in ClipKIT (ver. 2.1.1, see https://github.com/JLSteenwyk/ClipKIT; Steenwyk et al. 2020). Subsequently, the Paragone (ver. 1.0.0, see https://github.com/chrisjackson-pellicle/ParaGone; Yang and Smith 2014; Morales-Briones et al. 2021) pipeline was used to refine alignment quality. Paragone removes sequencing and alignment errors, and rogue taxa from alignments using HmmCleaner (ver. 0.243280, see https://metacpan.org/dist/Bio-MUST-Apps-HmmCleaner; Di Franco et al. 2019), TreeShrink (ver. 1.3.9, see https://github.com/uym2/TreeShrink; Mai and Mirarab 2018) and TrimAl (ver. 1.5.0, see http://trimal.cgenomics.org/; Capella-Gutiérrez et al. 2009). Aligned loci represented by fewer than 20 samples were removed from downstream analyses. Nuclear OzBaits and A353 targeted loci were combined to create the nuclear dataset for use in the subsequent analyses. Any duplicate loci sequenced by both nuclear bait sets were represented only once (i.e. the longest aligned locus was retained) in the final data set.

Phylogenetic analyses

Concatenated nuclear and plastid datasets were each analysed separately with a maximum likelihood approach in IQ-TREE (ver. 2.2.0, see http://www.iqtree.org/; Nguyen et al. 2015; Chernomor et al. 2016), with partitioning for individual loci. The best nucleotide substitution model for each partition was found using ModelFinder (see http://www.iqtree.org/ModelFinder/), as implemented in IQ-TREE (MFP + MERGE) and like partitions were subsequently merged (Kalyaanamoorthy et al. 2017). Only the top 10% of partition merging schemes were examined by using the relaxed hierarchical clustering algorithm to save computational time (Lanfear et al. 2014). Topological support was assessed using 1000 ultrafast bootstrap replicates (UFBS; Minh et al. 2013) and the approximate likelihood ratio test (ALRT; Guindon et al. 2010), as implemented in IQ-TREE. We interpreted support values following Minh et al. (2013) and recommendations in the IQ-TREE manual (see https://iqtree.github.io/doc/iqtree-doc.pdf). Additionally, we used IQ-TREE to produce individual locus trees of the nuclear dataset as input for a coalescent-based analysis using wASTRAL–hybrid (ver. 1.4.2.3, see https://github.com/smirarab/ASTRAL) in Accurate Species Tree Estimator (ASTER, hereafter referred to as ASTRAL) using the ‘likelihood’ preset (i.e. the ALRT value is used to weight branch support using default values). This method reimplements ASTRAL-III and accounts for phylogenetic uncertainty by considering signals derived from branch length and branch support in gene trees (Zhang et al. 2018; Zhang and Mirarab 2022). While a general concern for phylogenomic analyses is the presence of undetected gene copies, species tree inference using ASTRAL is considered robust to the presence of paralogs (Yan et al. 2022). Nevertheless, we conducted a further analysis using ASTRAL-Pro (ASTRAL for paralogs and orthologs, see https://github.com/chaoszhang/A-pro; Zhang et al. 2020) that attempts to identify the most likely species tree given a set of paralog gene trees. Here, we used the ‘unfiltered’ alignments from the CAPTUS align output to infer paralog gene trees using IQ-TREE (analysis settings as above) as input to ASTRAL-Pro and analyses were conducted using default settings.

The concatenated and coalescent trees of the nuclear dataset were both analysed using Quartet Sampling (QS, see https://github.com/fephyfofum/quartetsampling; Pease et al. 2018). QS is a quartet-based method that provides a wealth of information on how well the data support the topology and when support is poor, can indicate which underlying biological factor may be contributing to support values. Branch support for the QS analyses was assessed by the associated Quartet Concordance (QC), Quartet Differential (QD) and Quartet Informativeness (QI) values that were described in Pease et al. (2018). Branches with a QC value ≥0.5 were considered concordant, indicating that half or more of all quartet trees are concordant with the focal branch. QD values indicate whether neither of the discordant topologies are favoured (high QD) or if there is a skew to one discordant topology (low QD). QD can therefore provide evidence as to whether discordance is due to incomplete lineage sorting (ILS) or introgression, with the former suggested by high QD and the latter by low QD (Pease et al. 2018; Paetzold et al. 2019). QI reports the proportion of the replicates that are informative for the quartet and can be used to distinguish lack of support (low QI) from conflict (high QI). Additionally, terminal nodes were assessed by the associated Quartet Fidelity (QF) values. QF reports the proportion of replicates that included the taxon and resulted in a concordant topology and therefore can identify ‘rogue taxa’ (Pease et al. 2018) that may otherwise decay support values.

Dataset creation

Following assembly, 431 nuclear gene region alignments were recovered, including 343 A353 loci and 88 OzBaits loci. Prior to tree building, 39 loci were excluded from further analyses; these included alignments with fewer than 20 sequences (12 A353 loci) and loci overlapping between OzBaits and A353 (27 loci). The final nuclear dataset included 392 loci (331 A353 and 61 OzBaits loci) for Drosera and outgroups, with a concatenated length of ~259,000 bp, for a total of 95 terminals. Gene recovery for nuclear alignments was high, with a mean of ~229 loci represented for each sample. Across all samples, 50–380 nuclear loci (Supplementary Table S1) were represented with an ungapped sequence length of ~17,000–189,000 bp, with 62/92 samples having an ungapped length of ≥50% of the total alignment length (≥94,697 bp), excluding outgroups.

The final plastid dataset included 59 loci, with a concatenated length of ~41,000 bp across 80 terminals. Between 2 and 57 plastid loci were retained in the final dataset (Supplementary Table S1), with an average of ~35.7 loci represented for each sample. Samples with an ungapped length ≥50% of the total plastid alignment length (≥20,551 bp) were represented by 45/78 samples, excluding outgroups. Full recovery statistics for the nuclear and plastid datasets are available as Supplementary Table S1.

Phylogenetic analyses

IQ-TREE analyses for each of the concatenated nuclear and plastid datasets recovered trees with largely congruent topologies and with well-supported relationships between the major clades (Fig. 2). The ASTRAL tree of the nuclear dataset (Supplementary Fig. S1) and the ASTRAL-Pro nuclear species tree topology generated from the paralog gene alignments (Supplementary Fig. S2) were also found to have well-supported relationships between major clades that are consistent with those found by the IQ-TREE analyses. Differences between the trees are confined to relationships that are poorly supported amongst analyses. We confine our discussion to the IQ-TREE and ASTRAL results. Tree topologies that are mostly consistent with the infrageneric classification of Fleischmann et al. (2017a) were recovered by our analyses, with all subgenera and sections resolved as monophyletic in both nuclear and plastid datasets (Fig. 2). Unlike previous studies by Rivadavia et al. (2003) and Fleischmann et al. (2017a), subgenera Arcturia (D. arcturi) and Regiae (D. regia) are resolved as sister in IQ-TREE and ASTRAL analyses of the nuclear dataset (Fig. 3), with weak support for this relationship from the IQ-TREE analysis (Fig. 2; ALRT 45.7 and UFBS 84) but high support in the ASTRAL tree (Supplementary Fig. S1; Local Posterior Probability (LPP) 0.99). Drosera arcturi was not included in the plastid dataset due to poor locus recovery (Supplementary Table S1). Subgenus Ergaleium was well supported in the nuclear tree but was unsupported in the plastid tree, with low support (ALRT 74.2 and UFBS 84, Fig. 2) for the node subtending section Coelophylla from the remaining subgenus Ergaleium. Some internal nodes within sections Lasiocephala and Bryastrum were moderately to poorly supported in the nuclear trees (Fig. 2, 3). In the plastid tree (Fig. 2), many relationships within these sections were very poorly supported with low branch lengths, including clades encompassing D. patens Lowrie & Conran to D. trichocaulis (Diels) Lowrie & Conran (section Bryastrum), D. platystigma Lehm. to D. nitidula Planch. (section Bryastrum), and D. dilatatopetiolaris K.Kondo to D. subtilis N.G.Marchant (section Lasiocephala). In all analyses, section Phycopsis (D. binata) was well resolved as sister to section Ergaleium and, apart from a few nodes, there was high support for most relationships within section Ergaleium in both the nuclear and plastid IQ-TREE analyses (Fig. 2).

Fig. 2.

Comparative nuclear (left) and plastid (right) phylogenies using IQ-TREE analysis. Ultrafast bootstrap support (UFBS) and approximate likelihood ratio test values (ALRT) shown adjacent to nodes where values are <100. Subgeneric and sectional classifications follow Fleischmann et al. (2017a).

Fig. 3.

Comparative IQ-TREE (left) and ASTRAL (right) phylogenies from recovered nuclear sequence data using Quartet Sampling. Phylogenies are presented as cladograms for readability. Internal nodes are coloured based on the associated Quartet Concordance (QC) and Quartet Differential (QD) values; nodes with <0.5 QC have the associated QC, QD and Quartet Informativeness (QI) values reported, respectively. Terminal nodes are coloured based on the associated Quartet Fidelity (QF) values.

Concordance and discordance on the nuclear phylogenies analysed with QS analyses show a large degree of congruence between the IQ-TREE and ASTRAL analyses (Fig. 3). In some cases, discordant nodes in the QS analyses correspond to those highly supported by ALRT, UFBS or LPP (Fig. 2, Supplementary Fig. S1). In particular, the Drosera crown node has a very low QC value (−0.43) and a low QD value (0.35) (Fig. 3) that indicates that an alternate topology is favoured for this node. The node segregating section Coelophylla from the remainder of subgenus Ergaleium is also discordant (QC 0.36 and 0.38 on the IQ-TREE and ASTRAL trees respectively; Fig. 3). Interestingly, the sister relationship between subgenera Arcturia and Regiae was highly concordant in the QS analyses (Fig. 3), in contrast to the low support found by ALRT and UFBS (Fig. 2).

Both QS analyses showed low QC values for internal nodes within sections Bryastrum (10 IQ-TREE and ASTRAL nodes) and Lasiocephala (4 IQ-TREE and 7 ASTRAL nodes), and subgenus Drosera at the internal node bifurcating sections Arachnopus and Thelocalyx from section Prolifera. Comparatively few nodes with low QC values were found within section Ergaleium (8 IQ-TREE and 9 ASTRAL nodes). High discordance within clades corresponds to low QF scores found at terminal nodes, with lowest being in section Lasiocephala.

Biogeographic patterns within Australian Drosera

Patterns of disjunct sister clades occurring in South Africa (D. regia), south-eastern Australia (D. arcturi and D. murfetii) and New Zealand (D. arcturi) appear similar to other prominent Gondwanan Southern Hemisphere angiosperm families, e.g. Proteaceae (tribes Petrophileae and Persoonieae) (Barker et al. 2007) and Myrtaceae (tribes Myrteae and Syzygieae) (Thornhill et al. 2015). However recent molecular phylogenies (Thornhill et al. 2015; Lamont et al. 2024) highlight complex biogeographic histories, including both vicariance and intercontinental long distance dispersal events for these families. As support values are low for the sister relationship between D. arcturi and D. regia in the nuclear tree (Fig. 2), plastid sequence data were not recovered for D. arcturi (Fig. 2), taxon sampling is incomplete for subgenus Arcturia and molecular clock methods are not applied to our analyses, we interpret these results only as moderate support for an ancient Gondwanan relationship between subgenera Arcturia and Regiae, with a great deal of additional work remaining.

The highest levels of diversity and endemism at both sectional and species levels of Australian Drosera occur in the Southwest Australian Floristic Region (SWAFR), in part due to the combination of long-term, stable but heterogeneous habitats occurring on ancient, highly infertile soils (Groom and Lamont 2015; Brundrett 2021), with sections Bryastrum and Ergaleium representing most of the species diversity in this region. The single representative of section Bryastrum from eastern Australia and New Zealand represented in our analyses (D. pygmaea from South Australia) is sister to D. nitidula in the SWAFR in our nuclear phylogenies (Fig. 2, 3), and in a polytomy with D. nitidula and D. platystigma in the plastid phylogeny (Fig. 2). There is a second region of Drosera diversity for sections Lasiocephala and Arachnopus across the seasonally dry monsoonal tropics of northern Australia, particularly in the Kimberley regions, with section Prolifera confined to wetter regions of the wet tropics of far north Queensland (Lowrie 2013a, 2013b, 2013c; Nunn and Bourke 2017). The close molecular affinity between sections Arachnopus, Stelogyne and Prolifera first reported by Rivadavia et al. (2003) is supported in our analyses (Fig. 2); these taxa occur respectively in northern and central Australia, the SWAFR and in the wet tropics of Queensland but are poorly sampled in our analyses (section Arachnopus) or species depauperate (sections Prolifera and Stelogyne). Within section Prolifera, a close affinity between D. adelae and D. schizandra has been implied by phytochemical analyses (Schlauer et al. 2019), aligning with Fig. 2 and 3 that show a close affiliation between D. adelae, D. buubugujin and D. schizandra, with D. prolifera more distantly related. There is no overlap in the modern distribution of each species (Nunn and Bourke 2017; Mathieson and Thompson 2020) and despite D. schizandra and D. prolifera hybrids (D. × ‘Andromeda’) existing in cultivation, our QS analyses find no evidence of ILS or introgression between section Prolifera taxa (Fig. 3). These factors support a hypothesis of extended geographic isolation between section Prolifera taxa within the wet tropical region, with limited opportunity for introgression between species.

The dryland growing members of section Ergaleium are most diverse in the Mediterranean climatic regions of southern Australia, especially in the SWAFR, they (and Drosera species generally) are largely absent from chenopod-dominated vegetation types across this region (Specht and Specht 1999), although some taxa in the section extend to New Zealand and Asia (Lowrie et al. 2017a, 2017b). Our analyses support three major east–west disjunctions across the Nullarbor Plain within section Erglaeium. A close sister relationship is resolved by our analyses between D. macrantha ssp. planchoni in South Australia and D. macrantha in the SWAFR (Fig. 3). As D. macrantha ssp. planchonii occurs in both south-eastern Australia and the arid eastern extent of the SWAFR (Lowrie 2013a, 2013b, 2013c), this close trans-Nullarbor affiliation is possibly due to a comparatively recent dispersal event from the SWAFR to south-eastern Australia; however, denser sampling is required to further resolve the D. macrantha complex. The D. whittakeri complex (Fig. 3; clade extending from D. whittakeri to D. praefolia) represents a second radiation of section Ergaleium in south-east Australia; the closest relatives to that complex in our analyses are D. macrophylla, D. monantha and D. platypoda from the SWAFR but with a comparatively large molecular divergence, low support values and discordance for internal nodes (Fig. 2, 3). A third, well-supported radiation within the predominantly eastern Australian D. peltata complex (Fig. 3; clade extending from D. gracilis to D. auriculata) is resolved as a distant sister clade to D. yilgarnensis (Fig. 2, 3), a species endemic to the SWAFR correctly considered by Lowrie (2013a, 2013b, 2013c) and Lowrie et al. (2017a) to have a close affiliation with the D. peltata complex. This relationship potentially supports a model of ancient dispersal or vicariance for the broad D. peltata clade; however, additional sampling of D. bicolor in the SWAFR and D. lunata from eastern Australia and overseas is required to further resolve relationships within this complex. Most subgenus Ergaleium species are confined to either the eastern or western sides of the Nullarbor Plain, however, D. binata Labill., D. glanduligera Lehm., D. pygmaea DC. and D. stricticaulis (Diels) O.H.Sarg. (Conran and Lowrie 2007; Lowrie 2013a, 2013b, 2013c) also extend across the east–west continental divide. The degree of molecular divergence between populations across this significant biogeographic barrier (Crisp and Cook 2007) remains uncharacterised here for these taxa; however, duplicate sampling across the Nullarbor Plain has been conducted for many of these species (L. T. Williamson et al., unpubl. data) and will be discussed in a future publication.

Potential drivers of phylogenetic discordance

Despite differences in support between datasets, presumably reflecting differences in alignment length and evolutionary rate between these genomic compartments (Rokas and Carroll 2005; Asar et al. 2024), topologies for the nuclear and plastid trees are broadly congruent. Support values within sections Bryastrum and Lasiocephala are, however, consistently low in both nuclear and plastid trees (Fig. 2, 3). Similar patterns of poor within-clade support and discordance have been documented in other angiosperm groups; these are typically driven by ancient (Lin et al. 2019; Zhou et al. 2022) or modern (Scharmann et al. 2021) introgression events, ILS (Cai et al. 2021; McLay et al. 2023) and gene tree estimation errors (Cai et al. 2021).

The comparatively poor support values and inconsistent topologies resolved by IQ-TREE and ASTRAL analyses within subgenus Ergaleium sections Bryastrum and Lasiocephala may be driven by ILS (section Bryastrum) and introgression (sections Bryastrum and Lasiocephala) (Fig. 3). The pattern of introgression implied by low QC and QD values (Fig. 3) within section Lasiocephala is supported by common field observations of hybridisation between natural populations (Lowrie 1991). Introgression zones where D. stipularis Baleeiro, R.W.Jobson & R.L.Barrett and D. petiolaris R.Br. ex DC. hybridise are known in the wild (B. Ng, pers. comm.), and numerous other section Lasiocephala hybrids have been documented where preferred habitats and soil types co-occur (Lowrie 1991). Widespread hybridisation between section Lasiocephala taxa is also reported by carnivorous plant cultivators, with 30 known hybrids listed by International Carnivorous Plant Society (2024) in private collections. Although further study is needed, this potentially represents an example of introgression playing a major role in the evolution of a carnivorous plant genus, similar to Nepenthes (Nauheimer et al. 2019; Scharmann et al. 2021) and Sarracenia (Baldwin et al. 2023). Unlike the widespread introgression within Nepenthes, potential introgression appears to be more common in two sections in our analyses and within section Drosera, as reported by other authors (Wood 1955; Hoshi et al. 1994; Schnell 1995; Nakano et al. 2004; Pearman 2004; Schlauer and Fleischmann 2016).

Taxa within section Bryastrum also hybridise in both natural conditions and cultivation (Lowrie 2013a, 2013b, 2013c; International Carnivorous Plant Society 2024), with hybrid populations potentially able to persist indefinitely by clonal gemmae reproduction. Interpreting relationships among section Bryastrum taxa is limited by the same issues present for section Lasiocephala, namely poorly supported relationships (Fig. 2) reflecting extensive gene tree conflict underlying the species tree and the plastid topology (Fig. 3). These observed patterns of poor support and discordance within sections Bryastrum and Lasiocephala are confounding factors that require additional analyses (e.g. networks) and denser (population-level) sampling for resolution.