A framework phylogeny of the diverse guinea-flowers (Hibbertia, Dilleniaceae) using high-throughput sequence data

Taxon sampling

Representatives were sampled from all four subgenera of Hibbertia, with 95 species selected to encompass the likely backbone of the genus based on morphological characteristics and prior phylogenetic studies (Table 1). Species were assigned to informal species groups by two of the authors (T.A. Hammer and K.R. Thiele) to aid ongoing taxonomic revisions within the genus. Subgenus Hemistemma (~250 spp.) was represented by 59 species sampled from 26 informal species groups and subg. Hibbertia (~100 spp.) by 32 species sampled from 15 informal species groups. The two smaller subgenera, subg. Pachynema (11 spp.) and subg. Adrastaea (1 sp.) were represented by three and one species respectively. Approximately 27% of the currently accepted and formally named species in Australia were included. As sequencing for this study was part of the GAP initiative, Hibbertia species outside of Australia (e.g. Madagascar and New Caledonia) were excluded and are planned to be included in a forthcoming study with expanded sampling instead. Four publicly available accessions of Dilleniaceae were included as outgroups (Table 1): Dillenia alata (R.Br. ex DC.) Banks ex Martelli, Tetracera daemeliana F.Muell., Curatella americana L. and Acrotrema sp. (the latter two having A353 data only, and new sequence generated for a Dillenia alata sample). In total, 100 (99 spp.) and 98 (97 spp.) genomic samples were included for nuclear and plastid markers respectively.

DNA extraction, library preparation and sequencing

DNA extraction, library preparation and sequencing were performed at the Australian Genome Research Facility (Melbourne, Australia) as part of the GAP initiative by Bioplatforms Australia (Sydney, 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 quantity and quality were assessed using 1% E-gel with Sybr Safe dye (ThermoFisher) and concentrations assessed using Quantifluor dsDNA assay (Promega). DNA samples were fragmented enzymatically as part of the NEBNext Ultra II FS library preparation workflow. Angiosperms353 and OzBaits nuclear 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 of ~350 bp. Leftover DNA was sent back to the University of Adelaide and OzBaits plastid libraries were prepared again 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 perform an additional chloroplast bait enrichment and genome skim to capture a wider range of genomic components (Zeng et al. 2018).

Libraries prepared by AGRF were enriched using both the Angiosperms353 (A353; Johnson et al. 2019) 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. The libraries prepared in the University of Adelaide laboratory were enriched (16 plex) with the OzBaits_CP (Waycott et al. 2021) chloroplast probe kit (catalogue number 300196R.v5) at 65°C using the V5 chemistry. Post-capture libraries were amplified with indexed primers and cleaned up with Ampure (1.0×). Additionally, the unenriched library was amplified (eight cycles) with indexed p5 and p7 primers to generate a genome skim with the intention to sequence high copy loci such as 16S, 18S and ITS (e.g. Weitemier et al. 2014). Both Cp and genome skim libraries were quantified using a 4150 TapeStation System (Agilent) with High Sensitivity DNA Screen tapes. Libraries were pooled equimolar and size-selected (300–700 bp) on a Pippin Prep (Sage Sciences) 1.5% agarose 250–1.5-kb gel cassette. All sequencing was performed on a NovaSeq. 6000 (Illumina Inc., San Diego, CA, USA) with v1.5 chemistry and 150-bp paired-end reads.

Data processing

High-throughput 150-bp paired-end reads were imported into CLC Genomics Workbench (ver. 20.0.2, QIAGEN, see https://digitalinsights.qiagen.com/) for demultiplexing and trimming using a quality score limit of 0.05 (Phred score ~13). Reads for each individual were randomly sampled to 2.5 million paired-end reads. The sampled reads were de novo assembled by MEGAHIT (ver. 1.2.9, see https://github.com/voutcn/megahit; Li et al. 2015) using the ‘assemble’ function in the CAPTUS bioinformatics pipeline (ver. 1.0.0, see https://github.com/edgardomortiz/Captus; Raza et al. 2023) with the CAPSKIM default –k-list within CAPTUS (see https://edgardomortiz.github.io/captus.docs/assembly/assemble/) and with ‘–min-count’ and ‘–prune-level’ both set to three.

We used the CAPTUS ‘extract’ function to extract target gene regions using custom nuclear and plastid references. Dilleniaceae reference sequences of OzBaits gene targets were derived from the transcriptome assembly of Dillenia indica L. and Hibbertia grossulariifolia (Salisb.) Salisb. (sample codes EHNF and OKEF respectively that were generated as part of the OneKP project; Leebens-Mack et al. 2019) accessed by local BLAST searches against the transcriptome sequences using the OzBaits target gene region from Arabidopsis thaliana (Araport11, see https://phytozome-next.jgi.doe.gov/info/Athaliana_Araport11; Cheng et al. 2017). References for the A353 gene targets were sourced from H. grossulariifolia (as for OzBaits) and the Kew Tree of Life Explorer (Baker et al. 2022), consisting of Acrotrema sp. (ERR4180187), Dillenia alata (GAP 80947) and Tetracera daemeliana (ERR7599814). Plastid references were sourced from GenBank and included coding regions extracted from the plastid genome sequences of D. indica (NC_042740), D. turbinata Finet & Gagnep. (NC_062798) and Tetracera sarmentosa (L.) Vahl (NC_065058). We used a minimum score of 0.2, minimum identity of 75 and minimum coverage of 30 (for definitions, see https://edgardomortiz.github.io/captus.docs/) for both the nuclear and plastid data to match the contigs to the target references using Scipio (ver. 1.4, see http://www.webscipio.org; Keller et al. 2008). Dillenia alata (GAP 80947) and Tetracera daemeliana (GAP 80948) were previously sequenced by GAP for only A353 and were represented in the present plastid dataset by extracting off-target plastid genes.

We used the CAPTUS ‘align’ function to generate multiple sequence alignments (MSAs) for the extracted nuclear and plastid gene data. All extracted markers were aligned with MAFFT (ver. 7.511, see https://mafft.cbrc.jp/alignment/software/; Katoh and Standley 2013) using the ‘mafft_auto’ algorithm. We specified a maximum of three paralogs per sample for paralogy inference implemented through CAPTUS and used the ‘informed’ filter method to generate alignments of putative orthologs for each gene region, using ‘–tolerance 2’. Alignments were trimmed with ClipKIT (ver. 1.3.0, see https://github.com/JLSteenwyk/ClipKIT; Steenwyk et al. 2020) using default settings in CAPTUS. Subsequently, we used the Paragone (ver. 1.0.0, see https://github.com/chrisjackson-pellicle/ParaGone) pipeline to only refine alignment quality. Paragone uses 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) to remove sequencing and alignment errors, and rogue taxa from MSAs. Alignments with fewer than 20 sequences (≥80% missing samples) were removed from downstream analyses.

As output from CAPTUS, ‘NT’ (coding sequences in nucleotides) output was used for A353 and ‘genes-flanked’ (complete gene sequences, including exons and introns along with flanking regions) was used for OzBaits, noting that the OzBaits baits are generally designed on a single exon per gene and thus the flanking regions are potentially introns. Genes recovered by both the A353 and OzBaits bait sets were identified, and only one version of the gene was retained that was selected based on the largest number of included sequences. The resulting genes from both bait sets were subsequently combined to create the full nuclear dataset. We used the ‘NT’ sequences output from CAPTUS for the plastid dataset.

Phylogenetic analyses

The nuclear and plastid alignments were each partitioned and analysed using IQ-TREE (ver. 2.2.0, see http://www.iqtree.org/; Nguyen et al. 2015; Chernomor et al. 2016) to generate a concatenated maximum likelihood tree. ModelFinder, as implemented in IQ-TREE (MF + MERGE, see http://www.iqtree.org/ModelFinder/), was used to find the best nucleotide substitution model for each partition and merge like partitions (Kalyaanamoorthy et al. 2017). To save computational time, only the top 10% of partition merging schemes was examined by using the relaxed hierarchical clustering algorithm (Lanfear et al. 2014). Support for the topology 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. Interpretation of support values follows Minh et al. (2013) and recommendations in the IQ-TREE manual (see http://www.iqtree.org/doc/iqtree-doc.pdf). Tetracera daemeliana F.Muell. was used as the outgroup for all phylogenetic reconstructions given the resolution as sister to the remaining Dilleniaceae (Horn 2009). Individual locus trees were produced in IQ-TREE from the concatenated alignment with the approximate Bayes test to generate branch support values (Anisimova et al. 2011). The initial phylogenies were used as input for genesortR (see https://github.com/mongiardino/genesortR, accessed 30 January 2024; Mongiardino Koch 2021), an R script that sorts and subsamples phylogenomic datasets based on properties that quantify phylogenetic usefulness. The sorted top 300 and 60 loci were retained for downstream analyses for the nuclear and plastid datasets respectively. The previous steps generating the concatenated and loci trees using IQ-TREE were re-run using the newly subsampled datasets.

Locus trees were analysed using wASTRAL hybrid in Accurate Species Tree Estimator (ASTER, ver. 1.4.2.3, see https://github.com/chaoszhang/ASTER) that 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) for both genomic datasets. We followed the author’s recommendations (see https://github.com/chaoszhang/ASTER/blob/90f475fbdc79151b86c27f0881c549a1c6d230a4/tutorial/astral-hybrid.md) for maximum and minimum values for local Bayesian supports. The concatenated and coalescent (ASTRAL) trees were both analysed using Quartet Sampling (QS) (see https://github.com/fephyfofum/quartetsampling; Pease et al. 2018). QS indicates how well the data support the tree topology of the phylogeny. Branch support for the QS analyses was assessed by the Quartet Concordance (QC), Quartet Differential (QD) and Quartet Informativeness (QI) values that were described in Pease et al. (2018). Branches were considered concordant with the displayed topology if QC ≥ 0.5, indicating that half or more of all quartet trees are concordant with the focal branch. High QD means that neither of the discordant topologies are favoured versus a skew to one discordant topology at low QD values. QD can therefore provide evidence as to whether discordance is due to incomplete lineage sorting (ILS) or introgression, with the former indicated by high QD and the latter by low QD. QI reports the proportion of the replicates that are informative for the quartet and can be used to distinguish lack of support from conflict (low to high QI respectively). Additionally, terminal nodes (taxa) were assessed by Quartet Fidelity (QF) values. QF reports the proportion of replicates that included the taxon and resulted in a concordant topology (Pease et al. 2018). The suitability of using coalescent methods on the plastome has been questionable (e.g. Doyle 2022) but we have nevertheless included this analysis for comparison.

Dataset creation

Assembly of the nuclear dataset recovered 351 loci using the A353 bait set and 66 loci using OzBaits, of which 5 and 10 were respectively removed due to having fewer than 20 sequences or being overlapping loci. Before subsampling with genesortR, the nuclear dataset contained 346 A353 loci and 56 OzBaits loci, resulting in the combined alignment containing 402 loci and an overall concatenated length of 227,848 bp. Assembly of the plastid dataset recovered 79 loci that were all retained for the genesortR step and totalled 61,776 bp in length. Hibbertia dentata R.Br. ex DC. and Hibbertia fasciculiflora K.R.Thiele sequenced poorly for all bait sets and were excluded from all subsequent analyses. The H. dentata and H. fasciculiflora samples were respectively represented in the nuclear dataset by only 3 and 16% of loci and 1 and 7% of the total alignment length, and in the plastid dataset by 5 and 3% of loci and 1 and 0.6% of total alignment length.

The subsampled nuclear alignment of 300 loci and 98 sequenced samples had a concatenated length of 188,421 bp, and the subsampled plastid alignment of 60 loci and 96 sequenced samples had a concatenated length of 58,370 bp. The nuclear dataset, after subsampling, consisted of 258 loci derived from A353 and 42 loci derived from OzBaits. Output figures from genesortR are available as Supplementary Fig. S1 and S2. Gene recovery for samples in the nuclear dataset was reasonably high, with 81/98 samples represented in 80% or more of the loci (mean 82%), whereas 90/98 samples had an ungapped length of ≥50% of the total alignment (≥94,210 bp; mean 75%). In the plastid dataset, 71/96 samples were represented in 80% or more of the loci (mean 85%) and 82/96 samples had an ungapped length of ≥50% of the total alignment (≥29,185 bp; mean 69%). Coverage information for the genomic datasets is available in Supplementary Table S1.

Phylogenetic analyses

Topologies from all phylogenetic reconstructions were highly similar but not perfectly congruent, especially with relationships that were poorly supported. The IQ-TREE concatenated nuclear and chloroplast datasets exhibited overall high congruence in the phylogenies, with most incongruences observed in a reasonably small number of branches lacking robust support (Fig. 2). Likewise, the phylogenies reconstructed through a coalescent approach implemented by ASTRAL were mostly congruent with concatenated trees of the corresponding dataset (Supplementary Fig. S3, S4).

Fig. 2.

IQ-TREE analysis of the concatenated datasets for nrDNA (a) and cpDNA (b). Nodes are coloured based on values from approximate likelihood ratio test (ALRT) and ultrafast bootstrap (UFBS) supports. Nodes with red circles indicate that both values fail to support the node (ALRT < 80 and UFBS < 95) and orange nodes indicate that only one value supports the node. Unlabelled nodes are supported with ALRT ≥ 80 and UFBS ≥ 95. Taxa are coloured based on the geographic distribution within Australia, with regions indicated on the inset map. Hibbertia scandens (dark green) occupies the northern and eastern Australian regions. Branches of the four subgenera of Hibbertia are labelled i–iv. Major monophyletic clades of subg. Hibbertia and subg. Hemistemma represented in all analyses for both datasets are labelled clades A–N; finer dashed, unlabelled lines are species without clade relationships.

In all analyses of both genomic datasets, Hibbertia and the four subgenera were monophyletic with very high support, with subg. Adrastaea and subg. Hibbertia placed as sisters, and subg. Hemistemma and subg. Pachynema as sisters. The concatenated nuclear phylogeny was mostly well supported along the backbone apart from one node in subg. Hemistemma and two nodes in subg. Hibbertia (Fig. 2a). The concatenated plastid phylogeny was particularly poorly resolved on the backbone of subg. Hemistemma with very short branch lengths and had one poorly supported node on the backbone of subg. Hibbertia (Fig. 2b). Nevertheless, we were able to use these phylogenies to identify clades that were consistently monophyletic in every analysis to aid in the discussion of the relationships for which we can be reasonably certain. In subg. Hemistemma and subg. Hibbertia, ‘major clades’ have been labelled as A–I and J–N respectively, with the two smaller subgenera labelled as ‘subg. Adrastaea’ and ‘subg. Pachynema’ (Fig. 2). Hibbertia hermaniifolia subsp. recondita Toelken and H. scandens (Willd.) Gilg were consistently resolved outside of the other clades. Several of the relationships between the major clades were incongruent between the two genomic datasets, especially where the trees were poorly supported. This poor support was reflected in both concatenated phylogenies in many of the same areas, such as in subg. Hibbertia with the relationships between clades J, K and L, and along the whole backbone of subg. Hemistemma due to the many unresolved nodes in the plastid trees (Fig. 2, Supplementary Fig. S3, S4). There was genuine incongruence between the trees in some well-supported relationships, such as the different placements of clades D and M, H. scandens, H. hermanniifolia subsp. recondita, and placement of H. gracilipes Benth. within clade A. The clade having the highest number of nodes with poor support was consistently clade F, which included most of the species of subg. Hemistemma from eastern Australia (Fig. 2, Supplementary Fig. S3, S4).

The QS analyses of the concatenated phylogenies of Fig. 2 revealed more discordant nodes than those found poorly supported with UFBS and ALRT methods (Fig. 3). The QS analyses of the ASTRAL trees (Supplementary Fig. S5, S6) were similar to the other analyses. Hibbertia and the subgenera were well supported in the QS analyses. Although all labelled clades were well supported in the prior analyses, clades A and M on the nuclear tree were discordant in the QS analysis (Fig. 3a), with QC values of 0.46 and 0.087 respectively; the QD scores for these nodes were reasonably high (0.98 and 0.71 respectively), indicating that no particular topology was favoured. All other major clades in the nuclear tree were concordant in the QS analyses and all major clades were concordant on the plastid tree. In both genomic datasets, the same general patterns were observed in the QS analyses as for the IQ-TREE and ASTRAL analyses, with the backbones of the two largest subgenera being poorly supported or discordant, leading to uncertainty in the relationships between the major clades within subgenera. In both genomic datasets, there were a few discordant nodes in clades E, F and J, and also within clade H for the nuclear tree. Clade F had the most discordant nodes of any clade, with six discordant nodes on the nuclear tree (Fig. 3a) and four discordant nodes on the plastid tree (Fig. 3b). In total, 22 of the 32 discordant nodes on the nuclear tree had low QD scores (<0.7) indicating that an alternate topology is favoured for these relationships. In general, there were more discordant nodes for the nuclear tree than the plastid tree, possibly indicating more conflict between genes in the former dataset. Reasonably low QF scores (<0.7) were found for five terminals on the nuclear tree (Fig. 3). QF scores were lower in the plastid tree (14 terminals with <0.7) that could indicate more uncertainty in taxon placement due to lower variation relative to the nuclear data. In general, terminals with lower QF scores (<0.8) corresponded to relationships that show some degree of discordance.

Fig. 3.

Quartet Sampling analysis on the nuclear (left) and plastid (right) phylogenies of the IQ-TREE concatenated analysis (shown in Fig. 2). Phylogenies are presented as cladograms for readability. Internal nodes are coloured based on their Quartet Concordance (QC) and Quartet Differential (QD) values; nodes with <0.5 QC have the QC, QD and Quartet Informativeness (QI) values reported respectively. Terminal nodes are coloured based on Quartet Fidelity (QF) values. Clades are labelled as for Fig. 2.

Performance of the methods

The phylogenies presented here are the first published that focus on Hibbertia and the first to implement phylogenomic methods within Dilleniaceae, being based on 300 nuclear and 60 plastid loci for ~100 species. Implementation of the Genomics for Australian Plants (GAP) Phylogenomics protocols for the family was successful in generating a significant dataset with high quality phylogenetic resolution. Only two samples (H. dentata and H. fasciculiflora) failed among the 96 that were newly sequenced for this study, with most of the remaining samples having high recovery of loci for all three bait sets (Supplementary Table S1). Both nuclear and plastid genomic data were successful in reconstructing phylogenies that are largely congruent (Fig. 2), yet we have identified areas where both datasets consistently fail to satisfactorily resolve relationships due to a lack of phylogenetic signal or from conflicting loci. All four subgenera and all 14 of the identified major clades of subg. Hibbertia and subg. Hemistemma have been found to be consistently monophyletic with robust support in both datasets.

Subgenera and diagnostic features

The phylogenies presented here agree in large part with those of Horn (2005, 2009). As in Horn (2005, 2009), the two largest subgenera, subg. Hemistemma and subg. Hibbertia, are each sister to the smaller subgenera, subg. Pachynema and subg. Adrastaea respectively (Fig. 2). These relationships are interesting in that the two largest subgenera are more similar in most morphological characters to each other than either is to the smaller subgenera and vice versa for the small subgenera.

All species in Hibbertia subg. Pachynema are distinct in the genus in having mature plants that are seemingly leafless, the stems, that in some species are distinctly flattened, being the primary photosynthetic organ at maturity. Broad leaves are only present at the base of these plants at the juvenile or regrowth stages and are absent at maturity, the mature stems of the plants having scale-like leaves that are similar in appearance to floral bracts (Hammer 2022). Hibbertia subg. Adrastaea comprises the single species H. salicifolia (DC.) F.Muell. that despite the leafy habit is reasonably similar in floral morphology to the H. conspicua (J.Drumm. ex Harv.) Gilg species group of H. subg. Pachynema that consists of H. conspicua, H. goyderi F.Muell. and H. triquetra T.Hammer (represented here by H. conspicua). These four species have a cone-like androecium comprising two or three whorls of stamens enclosing the two-carpellate gynoecium at anthesis (Fig. 1i), occasionally with a ring of staminodes on the outside in some species. Other members of subg. Pachynema usually have two distinct whorls of stamens (occasionally the inner reduced to staminodes). Petal colour in subg. Pachynema is the most varied in Hibbertia, ranging from yellow or white to pink or dark red; besides three orange-flowered species in south-western Australian subg. Hibbertia, all other Hibbertia species have yellow (sometimes very pale) petals.

Floral morphology in subg. Hibbertia and subg. Hemistemma is varied (especially in the latter) but both subgenera contain species that have the stamens arranged all around the carpels (but not in distinct whorls as in the smaller subgenera); the carpels usually number three or five in subg. Hibbertia and two, three or five in subg. Hemistemma (Fig. 1b, d, j–l). In subg. Hibbertia, stamens may be free in a single or multi-layered ring around the carpels (Fig. 1k), free in gaps between the carpels (Fig. 1l) or fused in bundles between the carpels (Fig. 1j). In subg. Hemistemma, stamens are most commonly placed on one side of two carpels, where these are either more or less erect with the styles erect to spreading (Fig. 1e, f) or curved over the carpels as a hand of bananas with the styles short and curved beneath the stamens (Fig. 1a, g, h). Most Hibbertia species have solitary flowers but a few distinct groups of species in subg. Hemistemma have multi-flowered inflorescences as pseudocincinnae (Hammer and Thiele 2022), including the H. polystachya Benth. species group in south-western Australia (Fig. 1c; Thiele 2019b), the H. banksii (R.Br. ex DC.) Benth. species group (Fig. 1h; Toelken 2023) in the Australian monsoon tropics, and the non-Australian groups in Madagascar and New Caledonia (e.g. Veillon 1990). The two former groups are not closely related here (Fig. 2), indicating that the multi-flowered inflorescence is independently derived. Vestiture is another key feature that varies between these otherwise similar subgenera, with species in subg. Hibbertia having only simple hairs, whereas those in subg. Hemistemma have hairs that range from simple to hooked or fasciculate, the latter often being star-like to scale-like in appearance (Toelken 2010).

Clades and species groups

All informal taxonomic species groups represented by more than one sample were found to occur in the same major clade (Fig. 2); examples of these include H. lineata Steud. and H. verrucosa (Turcz.) Benth. from the H. lineata species group in clade A (Thiele 2017) and H. hypericoides (DC.) Benth. and H. silvestris Diels from the H. hypericoides species group in clade H (Thiele and Cockerton 2015). In some cases, major clades comprise a single species group; examples of these include the H. aspera DC. group (e.g. H. cinerea Toelken and H. radians (Toelken) T. Hammer; Hammer 2023b) corresponding to clade G, the H. banksii group (e.g. H. candicans (Hook.f.) Benth., H. hooglandii J.R.Wheeler and H. ledifolia Benth.; Toelken 2023) corresponding to clade I, and the H. commutata species group (e.g. H. montana Steud., H. ovata Steud. and H. pilosa Steud.; Thiele 2019a) corresponding to clade K. Some of these clades (e.g. J, K, L and N) are also reflected in the relationships recovered by Horn (2005).

Revisionary work in Hibbertia has until recently been split between eastern and northern Australia on the one hand, and south-western Australia on the other, due to there being little overlap between the species that occur in these regions. We show, however, that species or clades of Hibbertia are often closely related or nested within clades that occur elsewhere on the continent. A noteworthy example is H. ferox (from eastern Queensland) that is robustly supported as sister to H. drummondii (Turcz.) Gilg (from south-western Western Australia; Fig. 2). An emphasis on vegetative rather than floral characters and geographic bias has led previous authors to assume that the sharply pointed leaves of H. ferox indicate a close relationship to other pungent-leaved species (such as H. exutiacies and other members of the H. acicularis group from eastern Australia; Jackes 2018; Toelken 2024). Hibbertia drummondii and H. ferox, although differing markedly in morphological leaf characteristics, share important and overlooked morphological floral characteristics, with two glabrous carpels surrounded by a whorl of stamens (Fig. 1d; Hammer et al. 2019, p. 5), an unusual combination within subg. Hemistemma. Although H. drummondii does not have sharply pointed leaves, this trait is characteristic of some species in clades A and C–F, and has clearly evolved multiple times in subg. Hemistemma.

None of the a priori species groups are found to cross two or more major clades; however, not all species groups are monophyletic. For example, clade F (Fig. 2) includes a large, taxonomically challenging group of species from eastern Australia that have been informally classified by H. Toelken into the H. acicularis (Labill.) F.Muell., H. sericea (R.Br. ex DC.) Benth., H. stricta (DC.) R.Br. ex F.Muell. and H. vestita A.Cunn. ex Benth. species groups (e.g. Toelken 2000, 2013, 2024). However, the H. acicularis species group (represented by H. exutiacies N.A.Wakef. and H. obtusibracteata Toelken) is clearly not monophyletic. Clade F also contains a small group of Western Australia species (i.e. H. charlesii J.R.Wheeler and H. pachyphylla J.R.Wheeler) that are well resolved as sister to the eastern Australian species in every analysis.

Species in clade F are morphologically diverse, covering a broad diversity of floral and vegetative character states in Hibbertia. The eastern component of clade F, although represented by only 13 species in this study (Fig. 2), is likely to include ~100 current species and ~100 putative new species (H. R. Toelken, pers. comm.), making this the most morphologically diverse and species-rich radiation of Hibbertia in eastern Australia. The study by Raheem (2012) treated the eastern Australian species of subg. Hemistemma as monophyletic but we resolve some subg. Hemistemma species from subtropical and temperate eastern Australia in clades other than clade F (e.g. H. aspera group, H. hermaniifolia subsp. recondita and H. ferox Jackes). The difficulties encountered by Raheem (2012) in finding support for morphological species groups and resolving trees based on ITS and trnL-F alone are not surprising given our results using far more data.

Possible causes of discordance

Biological factors affecting phylogenetic signal within genomic datasets is a subject of intense study, with many groups of organisms having recalcitrant relationships even with large datasets (e.g. Roycroft et al. 2019; Wikström et al. 2020; Nge et al. 2021; Yang et al. 2021; McLay et al. 2023). This may also be the case within Hibbertia given that both genomic datasets show poor support and discordance along the backbones of the two largest subgenera and within a few complex groups (e.g. clade F; Fig. 2 and 3). The QS analysis indicates that some of this discordance could be due to non-tree-like evolution (e.g. due to ILS and introgression) associated with short branches (e.g. Paetzold et al. 2019; Kong et al. 2022). Internal nodes with low QC and high QD (coloured red on Fig. 3) show discordant topologies with a low or no skew towards a particular topology that can indicate ILS (Pease et al. 2018; Paetzold et al. 2019). Internal nodes with low QC and low QD (coloured red on Fig. 3) show discordant topologies with high skew towards a particular topology that can indicate the presence of a secondary evolutionary history (e.g. due to introgression; Pease et al. 2018; Paetzold et al. 2019). Clade F, for example, contains nodes showing both of these patterns, that in combination could significantly blur the evolutionary relationships and account for discordance and poor support values in this clade (e.g. McLay et al. 2023). A complex evolutionary history of clade F, possibly including introgression and rapid radiation, could also account for the difficulties encountered with working on the complex taxonomic classification of the group. In comparison, species groups in subg. Hemistemma in south-western Australia, although very diverse, have not had the same taxonomic issues and are reasonably well defined with readily identifiable characters. These groups are also represented in clades of the phylogeny with usually longer branches and higher support values (Fig. 2 and 3). The difference between clade F and these clades from south-western Australia could be due to differing diversification dynamics. Many other plant groups also show this pattern between south-eastern and south-western Australia that may be attributed to the historically more stable environment in the south-west (Nge et al. 2020). Future work will attempt to elucidate biological factors influencing relationships recovered by phylogenies and explore the historical complexities of the evolution of Hibbertia.