Using RADseq to resolve species boundaries in a morphologically complex group of yellow-flowered shrubs (Geleznowia, Rutaceae)

ABSTRACT The morphologically complex and charismatic genus Geleznowia (Rutaceae) is endemic to south-western Australia and faces existing and potential conservation issues associated with land clearing, climate change and commercial harvesting. Two species are currently recognised in the genus, but horticulturally recognised forms and phrase-named taxa reflect additional suspected species diversity. The genus exhibits complicated and subtle patterns of morphological variation that have historically inhibited delimitation of taxonomic entities and, as a result, precluded effective conservation assessments. Here we used ddRAD data from 25 populations across the range of Geleznowia to elucidate genomic diversity in the group in conjunction with morphological re-assessment so as to delimit species and revise the taxonomy. Our analyses consistently identified seven entities that maintain genomic distinctiveness even in sympatry with other entities, supporting the inference of reproductive barriers and lineage divergence. Morphological assessment of more than 300 specimens corroborated these seven taxa. Consequently, we recognise seven species of Geleznowia in Western Australia, retaining G. amabilis K.A.Sheph. & A.D.Crawford, recircumscribing G. verrucosa Turcz., reinstating G. calycina (J.Drumm. ex Harv.) Benth., and naming four new species as G. eximia K.A.Sheph. & A.D.Crawford, G. narcissoides K.A.Sheph. & A.D.Crawford, G. occulta K.A.Sheph. & A.D.Crawford, and G. uberiflora K.A.Sheph. & A.D.Crawford.


Introduction
The hyperdiverse Southwest Australian Floristic Region (SWAFR) harbours ~7700 native species of angiosperms (see the Western Australian Herbarium's Florabase, https:// florabase.dbca.wa.gov.au/) in an area approximately the size of Italy (see Hopper and Gioia 2004; Department of Agriculture, Water and the Environment 2020). Patterns of diversity in the south-west include both geographically restricted species and widespread species, often with strong geographic structuring of genetic diversity that has been suggested to reflect complex patterns of localised extinction and persistence in response to climatic fluctuations (Byrne 2008). The contrasting narrow and broad structuring of diversity and variable levels of genetic connectivity make taxonomic resolution of species boundaries and identification of distinct genetic groups crucial for effective conservation decisions amid large-scale human impacts in the region (Coates 2000).
Restricted to the northern sandplains of south-western Western Australia (WA), yellow-flowered shrubs in the genus Geleznowia Turcz. (Rutaceae) show indications of this varied and complicated geographic structuring of diversity. These plants often depend on disturbance for recruitment, and are frequently observed after fire, in gravel pits and along road verges. Despite often occurring in large numbers following disturbance, the plants appear to be transient in the landscape, with most senescing within a decade and only a few persisting longer. As their floral displays are stunning ( Fig. 1), wild populations are targeted for harvesting by the cut-flower industry. Geleznowia has not been brought into commercial cultivation, despite work demonstrating horticultural potential (Plummer et al. 2000), because propagating them from cuttings in sustainable numbers has been shown to be challenging (D. Growns, pers. comm.). Some taxa appear to be narrow-range endemics because they are known only from a few populations. With their scarcity and patchy distribution, these plants are potentially under threat from over-harvesting and habitat decline through climate change, weed invasion, altered fire regimes, and further development in the region, but taxonomic uncertainties prevent formal conservation assessments and development of appropriate management plans. Clarifying the taxonomy of Geleznowia is therefore essential to guide the conservation of this charismatic genus.
The taxonomic history of Geleznowia reflects some of the uncertainty about species boundaries. The genus was described by Turczaninow (1849) with a single species, G. verrucosa Turcz. Later, Bentham (1863) also recognised G. calycina (J.Drumm. ex Harv.) Benth. and G. macrocarpa Benth. Taxonomic usage of the names since then has synonymised G. calycina and G. macrocarpa under G. verrucosa, retaining Geleznowia as a monotypic genus (Wilson 1970; Council of Heads of Australian Herbaria 2007). Genetic and morphological work by Broadhurst (Broadhurst et al. 1999Broadhurst 2000) suggested the presence of at least two subspecies, with a third form hypothesised to have resulted from hybridisation between the other two; however, no taxonomic changes were published. At approximately the same time, a horticultural study (Plummer et al. 2000) recognised five distinct morphological forms that retained differences when grown together in a common garden. To recognise and provide legislative protection for some of these taxa, informal phrase names (see Barker 2005) were created in 2010 for three of the forms, namely, G. sp. Red Bluff (A. Crawford ADC 597), G. sp. Marchagee (A. Crawford ADC 1353) and G. sp. Binnu (K.A. Shepherd & J. Wege KS 1301). Although G. verrucosa and G. sp. Marchagee were reasonably widespread, G. sp. Red Bluff and G. sp. Binnu were range restricted and respectively given Priority Two and Priority Three statuses, indicating that they may warrant listing as Threatened but currently do not meet the criteria owing to insufficient data, and thus require further survey. Noting considerable variation in indumentum, leaf size and flowers in his Flora of Australia treatment, Wilson (2013) still recognised a single variable species, G. verrucosa, suggesting infraspecific variation might eventually be delimited. Although published in 2013, it is likely that his treatment pre-dated the erection of the phrase-named taxa and did not refer to them.
More recently, Geleznowia sp. Red Bluff was recognised as G. amabilis K.A.Sheph. & A.D.Crawford, retaining its Priority Two conservation status because it is currently known only from a few populations in the Kalbarri region (Shepherd and Crawford 2020). However, the resolution of G. sp. Binnu and G. sp. Marchagee, as well as a geographically restricted and morphologically divergent population designated here as 'White Peak' (Fig. 1), has faced challenges. It has proven difficult to confidently separate entities within what appears to be a species complex on the basis of the examination of herbarium specimens alone, partly because of subtle character differences that sometimes overlap because of variation during plant growth stages. Given the current and historical uncertainty about morphological taxa and putative hybridisation within this complex, genomic data represent an ideal tool to shed light on evolution in the group and help delimit species. We undertook a large field-based study to collect specimens and material for genomic analysis and applied a method based on a reduced representation of the genome through sequencing of restriction site-associated DNA (RAD: Baird et al. 2008;ddRAD: Peterson et al. 2012). The method provides thousands of single-nucleotide polymorphisms (SNPs) useful for both population genetic and phylogenetic analyses. RAD methods have been used to effectively evaluate species boundaries in other plant genera, e.g. Diospyros L. (Linan et al. 2021), Leptospermum J.R.Forst. & G.Forst. (Binks and Byrne 2022), Rhodanthemum B.H.Wilcox, K.Bremer & Humphries (Wagner et al. 2020), Seringia J.Gay (Binks et al. 2020), Triodia R.Br. (Anderson et al. 2017) and Viburnum L. (Spriggs et al. 2019).
In approaching species delimitation in a complex, we adopt a species concept following ideas developed by de Queiroz (1998de Queiroz ( , 2007 that species represent evolutionarily distinct lineages that acquire distinguishing features (e.g. reciprocal monophyly, fixed morphological differences) over time since speciation. Lineage recognition is then based on detecting one or a combination of these distinguishing features, with greater confidence for a given delimitation when more lines of evidence agree. We approach species delimitation as an integrative task (e.g. Dayrat 2005;Padial et al. 2010), comparing genomic evidence with morphological and ecological evidence to support the distinction of taxa.
This study aims to use genomic data from populations across the range of morphological forms in Geleznowia (G. amabilis, G. verrucosa, G. sp. Binnu, G. sp. Marchagee and G. 'White Peak') (Fig. 2) in combination with morphological measurements and field observations to delimit species and revise the taxonomy, highlighting patterns of evolution in the group and informing their conservation. Geleznowia forms recognised at the start of this study ( Fig. 1, 2, Supplementary Table S1). For DNA sampling, fresh leaf material was collected from 142 plants in 2020, with six individuals per form per population except in two cases (P9 and P10) where only two individuals were present. Molecular data indicated that Geleznowia is closely related to Philotheca Rudge & Drummondita Harv. (Bayly et al. 2013;Appelhans et al. 2021), so these were chosen as outgroups. We included six fresh collections each of two populations of Drummondita ericoides Harv. (Fig. 1g), two of Philotheca spicata (A.Rich.) Paul G.Wilson, one of P. pinoides (Paul G.Wilson) Paul G.Wilson, and one of P. sericea (Paul G.Wilson) Paul G. Wilson (Fig. 1f). Fresh material was collected into silica gel for immediate dehydration and freeze-dried on return to the laboratory. At least one voucher specimen per taxon per population was submitted to PERTH (Table S1).

Library preparation, sequencing, processing and assembly
Extraction of DNA followed a modified CTAB protocol (Doyle and Dickson 1987), with the addition of polyvinylpyrrolidone (PVP) to the extraction buffer. Approximately 40 mg of dried leaf material was used. Aliquots were taken twice from 13 samples to serve as technical replicates (bringing the total to 191 samples). DNA was sent to the Australian Genome Research Facility (AGRF) in Melbourne for library preparation and sequencing. Library preparation was based on a modified double-digest restriction site-associated DNA (ddRAD) protocol (Peterson et al. 2012). Modifications included slight alterations to digestion conditions, moving clean-up steps to after ligation, and including a column clean-up step. Initial screening of eight digestion enzyme combinations on three samples was followed by library generation for 191 samples by using the optimum combination (PstI/MspI) chosen based on the amount of amplified product and the absence of repeat regions. Following enzyme digestion and adapter ligation (including a barcode), samples were pooled in groups of 48 for size selection (280-375 bp) and 11 cycles of indexing polymerase chain reaction (PCR) for seven replicates (separated then re-pooled) of the group. All samples were then pooled into a single sequencing run (an entire SP200 flow cell = 2 lanes) in a 150 bp single-end configuration on an Illumina NovaSeq 6000 (ver. 1.0, see https://sapac. illumina.com/systems/sequencing-platforms/novaseq.html). Raw reads were returned in multiplexed files based on four indices.
Demultiplexing of reads was performed using ipyrad (ver. 0.9.81, see https://github.com/dereneaton/ipyrad; Eaton and Overcast 2020), including filtering for adapter contamination and low-quality base pairs, allowing a maximum barcode mismatch of one, a minimum quality score of 30 for trimming the end of reads, and a maximum number of low-quality bases per trimmed read of two. Any reads less than 50 bp long after trimming were removed. After demultiplexing and quality filtering, samples were compared for the number of retained reads. If the number of retained reads was less than two standard deviations below the average number of retained reads across samples or if it was less than 800 000, the sample was removed from the full dataset. This resulted in the removal of four samples (including two of the technical replicates).
To choose optimal assembly parameters in ipyrad, we ran assemblies for a subset of samples (two samples per population, including all original technical replicates, for a total of 73 samples) for each hundredth of the clustering threshold from 0.85 to 0.99. The resulting assemblies were then compared for a range of parameters suggested by earlier studies (Ilut et al. 2014;Mastretta-Yanes et al. 2015;Paris et al. 2017;McCartney-Melstad et al. 2019). Based on this, a clustering threshold of 0.93 was selected for assembling the full dataset.
The full assembly used ipyrad parameters corresponding to a minimum read depth for base calls of six, maximum alleles in a consensus sequence of two (diploids; n = 14; Smith-White 1954), maximum proportion of uncertain base calls and maximum proportion of heterozygous base calls in consensus sequences of 0.05, maximum proportion of SNPs in a locus of 0.2, and maximum proportion of samples sharing a heterozygous site of 0.5. We kept all SNPs that were present in at least three samples to allow flexibility in downstream filtering. The output variant call format file (VCF; containing SNPs) and loci file (full alignments) were used for downstream filtering and analyses.

Filtering
To estimate SNP error rate in relation to read depth, we used Tiger (ver. 1.0, see https://bitbucket.org/wegmannlab/tiger; Bresadola et al. 2020) and SNPs from 11 technical replicates. The heterozygous error rate was used to set minimum (17) and maximum (100) mean depth cut-offs for retained loci to avoid higher error rates at low and high depths. Following error-rate estimation, replicates were removed in downstream filters, typically keeping the sample with more retained loci in the dataset.  To assess whether there were any samples that could not be distinguished genomically (effective clones), ingroup samples were filtered for SNPs present in at least 50% of samples, and pairwise SNP comparisons were made using a custom Python script (ver. 3.8, see https://docs.python.org/3/reference/). The proportion of SNPs (called in both) differing between samples for each of the nine technical replicate pairs was determined. These proportions were averaged for all replicates and doubled to give a conservative approximate lower limit for distinguishing between samples (~0.8% difference). Pairwise comparisons with fewer differences than this limit (i.e. <0.8% different) represent effective clones, and these were randomly reduced to a single representative per clone. This resulted in the removal of 57 ingroup samples, in some cases reducing populations to a single individual or, in one case, two populations (P21 and P22) to a single individual. Consequently, downstream analysis focused on individual-based metrics of genomic relationships, rather than population-based metrics.
Custom scripts (see Data availability) were then used to generate four datasets for various downstream analyses on the basis of differing filtering and sample combinations (Table 1). For ordination and allele-frequency analyses (Dataset 1), stricter filter settings were used, requiring SNPs to be present in at least 90% of samples, to be bi-allelic and to have a minimum minor allele count of three. A single SNP was retained per locus (chosen on the basis of maximum sample coverage, and randomly in the case of ties) to reduce the potential for linkage. For analyses less sensitive to missing data (genomic distances and phylogenetic analyses; Datasets 2-4), laxer filter settings were used, requiring SNPs to be present in at least 25 or 50% of samples (see Table 1), and without restrictions on allele count.
Dataset 1 was run under an admixture model with no priors, and 40 runs per K value for K1-15. Each run included a burnin of 100 000 generations of the Markov-chain Monte Carlo (MCMC) and sampling for a further 100 000 generations. The optimal K value corresponding to the highest level of structure in the data was evaluated with the top 10 runs for each K on the basis of likelihood by using the Evanno delta K method (Evanno et al. 2005). The top 10 runs were submitted to the online web service for CLUMPAK (see http://clumpak.tau.ac. il/index.html; Kopelman et al. 2015) to obtain the ancestry proportions for major and minor modes in the runs.
Genomic distances between samples were calculated using Dataset 2. The filtered VCF was imported and converted using vcfR into a DNAbin object for the package ape (ver. 5.6.1, see https://cran.r-project.org/package=ape; Paradis and Schliep 2019). Distances between samples were calculated using the MATCHSTATES method in the R package pofadinr (ver. 1.0, see https://github.com/simjoly/ pofadinr). The method uses matches in nucleotides between two samples at each position (Joly et al. 2015) and can account for missing data. The distance matrix was visualised as a network in SplitsTree4 (ver. 4.17.1, see https://softwareab.cs.uni-tuebingen.de/download/splitstree4/welcome.html; Huson and Bryant 2006), using the Neighbour-Net method (Bryant and Moulton 2003), with splits displayed using the Equal Angle algorithm (Dress and Huson 2004).
Phylogenetic relationships among samples were estimated using Datasets 3 and 4 as input to concatenation and coalescent approaches respectively. For concatenation, loci from the filtered VCF were extracted from the full-loci file, samples were filtered from those loci, and alignments were combined and filtered for a maximum 50% missing data per site. The resulting alignment was used in a maximum-likelihood analysis in IQ-TREE (ver. 2.1.3, see http://www.iqtree.org/; Nguyen et al. 2015;Minh et al. 2020a). IQ-TREE was run with a search for the optimal model (Kalyaanamoorthy et al. 2017), 1000 ultrafast bootstrap (UFB) replicates (Hoang et al. 2018), and 1000 bootstrap replicates of the Shimodaira-Hasegawa-like (SH) approximate likelihood-ratio test (Guindon et al. 2010).
For a coalescent approach, a single SNP per locus was run in SVDquartets (Chifman and Kubatko 2014) as implemented in PAUP* (ver. 4.0a, build 168; see https://paup. phylosolutions.com/; Swofford 2003). The SNPs were concatenated into a single sequence per sample, with ambiguity codes being used for heterozygous sites. Two analyses were run, namely (1) individual lineages (i.e. each sample), and (2) samples grouped by putative species on the basis of Rectangles surrounding population labels and symbols are coloured by initial morphological identification corresponding to forms recognised at the outset of this study. Populations are numbered, increasing south to north. Smaller symbols correspond to herbarium specimens and are coloured similarly. The base map is sourced from results from PCA, genomic distances and Structure. Each SVDquartets analysis was run to sample quartets exhaustively, with ambiguous counts being distributed over compatible site patterns, with quartets assembled using QFM (Reaz et al. 2014), and with 100 bootstrap replicates. To evaluate conflicting phylogenetic signal, site concordance factors (Minh et al. 2020b) were calculated for the lineages tree by using the full alignment in IQ-TREE and sampling 100 000 random quartets around each internal branch. Site concordance factors indicate the percentage of informative sites for a given branch that support that branch in the tree. Trees were plotted using the R package ape. To visualise the conflict among bootstrap replicates, bootstrap trees were used to construct consensus networks (Holland and Moulton 2003) in SplitsTree4, including splits found in at least 10 of the 100 replicates.

Morphological measurements
Thirteen morphological characters were used for taxon comparison (see Results), with additional characters measured for taxonomic descriptions (see Taxonomic treatment) by using dried herbarium specimens (380 sheets held at PERTH including 172 vouchers collected in 2020 and 2021). Specimens were assessed to ensure that they fit the character ranges for a given taxon, but individual measurements were not recorded for all specimens. Detailed floral measurements (at anthesis) were made from field material preserved in 70% ethanol or rehydrated herbarium material by using boiling water with a drop of detergent. Assessment of the whorls of yellow petaloid bracts subtending each flower ( Fig. 1) can be challenging because of a gradual transition from green leaves to yellow bracts (with some structures being both green and yellow). Each bract was scored as such only when there was an increase in size from the leaves below and the yellow colour extended for more than half its length. Seed characters (see Supplementary Fig. S1) were measured using SmartGrain (ver. 1.2, see https://www.naro. affrc.go.jp/archive/nias/qtl/SmartGrain/; Tanabata et al. 2012) on seed collections held in the Western Australian Seed Centre (Kensington) for which the identity at the collecting locality could be verified by a specimen lodged at PERTH.
Examination of type material from MEL was facilitated through a loan to PERTH. Scanned images of type material were also viewed in Global Plants (http://plants.jstor.org/) and the Museum National d'Histoire Naturelle (P) online database (http://science.mnhn.fr/institution/mnhn/search).

Genomic analyses
The PCA (Fig. 3) generally showed continuous variation across the samples, but distinct clusters were evident (sometimes corresponding with morphology) in some dimensions. The first four principal components of the PCA explained similarly substantive variation (~40% of the total genomic variation) and allowed for the identification of six cohesive genomic clusters. The first and second components separated a group of G. verrucosa samples (Group 1) from the rest (including other G. verrucosa samples). In the third principal component, the G. sp. Binnu samples formed a distinct cluster (Group 2) from everything else. The fourth principal component more clearly demarcated a second group of G. verrucosa samples (Group 3) and separated two groups of G. sp. Marchagee samples (Groups 4 and 5) as distinct. Finally, the second and third components separated G. amabilis samples as a distinct cluster (Group 6). The morphologically divergent G. 'White Peak' (P9; two samples) and a single population of G. verrucosa (P4; one sample) consistently clustered near other groups and did not fall out as clearly distinct. For the remaining analyses, we refer to these six genomic groups and two divergent populations (P4 and P9). Structure output (Fig. 4) suggested K6 as optimal or the highest level of structure on the basis of the Evanno method. At K6, the same six genomic groups were evident as in the PCA, with the majority of samples showing no admixture among them. Some admixture was evident in a few individuals in Groups 1 and 2, with one sample of Group 2 that co-occurs with plants from Group 1 showing almost half of the assigned ancestry to Group 1, consistent with hybridisation. Other samples of Group 2 showed a smaller proportion of shared ancestry, and one individual from the Group 1 population P15 showed a small proportion of admixture with Group 2. One sample from Group 3 showed indications of possible mixed ancestry, and the divergent population P4 was unresolved with ancestry assigned to multiple groups. At K6, the divergent population P9 showed admixture with Groups 2, 4 and 5, yet, at K7, it was resolved as a distinct cluster. At K8, further population structure was shown within Group 5, consistent with more continuous variation.
Whereas major modes detected by CLUMPAK are shown in Fig. 4, minor modes were also detected (see Fig. S2), although they did not largely differ from the overall pattern recovered and showed the same six groups at K6.
The SplitsTree4 network ( Fig. 5) showed separation of the same six genomic clusters as those present in the PCA and Structure results, along with the two divergent populations (P4 and P9, both on isolated branches). Populations typically grouped distinctly within each of the six clusters, but not in all cases. Again, the G. verrucosa groups (1 and 3) were separated, and there was a deep division within the G. sp. Marchagee sampling (Groups 4 and 5). There were long splits connecting co-occurring populations P15 (in Group 1) and P16 (in Group 2), consistent with hybridisation.
The concatenated alignment was 442 161 bp, with 13 369 parsimony informative sites. Samples had an average of 24.6% (8.2-82.0%) gaps or ambiguities; outgroup samples all had >75% and ingroup samples all had <20%. The bestfit model chosen by IQ-TREE on the basis of the Bayesian information criterion was TPM3 + F + R7 (AC = CG, AG = CT, AT = GT + empirical base frequencies + seven free rate categories). In the resulting tree ( Fig. 6), all six of the genomic groups from the PCA were fully supported as monophyletic (100/100 SH/UFB). Relationships between groups were generally poorly supported for a dataset of this size, except for a sister relationship between Groups 1 and 3. The divergent population P4 was supported (100/ 100) as sister to Group 4, whereas P9 had some support (91.1/95) as sister to Group 5. Branching from the backbone was typically short compared with divergence from outgroups and reflected a lack of information about relationships between groups in the ingroup.
For the SVDquartets coalescent analyses, the lineages and species trees (Fig. 7) had 16.7 and 26.7% of quartets respectively, incompatible with the resulting trees. The individual lineages analysis recovered good support for the monophyly of most of the genomic groups from the PCA, as evident in high bootstrap values and concordance factors indicating most sites supported that grouping. Groups 2 and 3 had 99% rather than 100% bootstrap support, whereas Group 1 was not recovered as monophyletic because of the odd placement of one sample from P17 (asterisk in Fig. 7), which grouped with the outgroup, although with significant conflicting signal (see the concordance pie chart at the node subtending the remainder of the ingroup). Relationships www.publish.csiro.au/sb Australian Systematic Botany among groups had some support. Groups 2 and 6 were supported as sister (100% bootstrap support), whereas Groups 4 and 5 formed a clade together with P4 and P9 (93%). The divergent population P4 was recovered as sister to Group 4 (80%). In the species analysis, there was strong support for the sister relationships of Groups 2 and 6 (96%), and low support for a sister relationship between Groups 1 and 3 (70%). On the basis of the bootstrap tree network, support for Groups 1 and 3 as sister was potentially affected by signal linking Group 1 and the outgroup, as in the lineages tree. The relationship between Groups 4 and 5 as belonging to a single clade was poorly supported (55%), and it was likely that this was affected by the placement of P9. Groups 4 and 5 were recovered in a single clade with P4 and P9 in 74% of the bootstrap trees.

Morphological measurements
Floral and vegetative characters were compared for each genomic group identified in the genomic analyses and similar associated specimens (see Table 2). The comparisons showed unique sets of characters that could be used to demarcate each of the six genomic groups, plus the morphologically divergent population P9 (see the key in the Taxonomic treatment). Diagnostic characters included the number, shape, size and indumentum of the yellow petaloid bracts subtending each flower (see Fig. 1), the total number of flowers per inflorescence, and the size and shape of the sepals.

Discussion
Our analyses of Geleznowia have demonstrated again the potential of genomic data to distinguish separate lineages within morphologically variable and taxonomically challenging groups of plants. Population sampling of ddRAD genomic data across the genus enabled the identification of multiple distinct entities that largely maintain their distinctiveness when in sympatry. Detailed morphological review and comparison of those entities uncovered sets of characters that can be used to distinguish them. Together, our results point to unrecognised lineage diversity in the group that warrants taxonomic recognition and conservation assessment.

Species delimitation
Consistent patterns across ordination, distance, Bayesian and phylogenetic analyses, coupled with morphological measurements, point to seven morphologically recognisable and genomically distinct groups within Geleznowia. The challenge for understanding species complexes is to distinguish the fuzzy boundary where population differentiation has transitioned to effective reproductive isolation and independent evolution that we expect for species (e.g. de Queiroz 2007). One of the strongest indicators of this is co-occurrence of taxa with the maintenance of morphological, genomic or ecological differences (Harrison and Larson 2014;Rannala and Yang 2020). Within Geleznowia, this is reflected in the region near Kalbarri, where three of the groups (3, 5 and 6) occur within short distances of each other and sometimes at the same location. Field observations also include members of Group 1 co-occurring with Groups 3 and 5 in this area (the closest sampled population is P17). These groups were recovered as distinct in all analyses and represent the extremes of morphological and genomic variation in the group. Their close occurrence without admixture points to established reproductive barriers and independent evolution. Co-occurrence or nearby populations and a lack of admixture is likely to be caused by processes that result in reproductive isolation, e.g. flowering-time differences, different pollinators, and reproductive incompatibilities, among others (Grant 1971). The pattern might also result from a breeding system that reduces outcrossing (e.g. Martin and Willis 2007). Little is known about what pollinates Geleznowia, with only a few observations of moths, ants and beetles having been reported on the flowers, none of which appeared to result in pollination . The showy flowers are superficially alike and do not appear to exhibit differing specialisations for specific pollinators among the groups. Work on the breeding system and genetic diversity in morphological forms of Geleznowia (Broadhurst et al. 1999; suggested facultative selfing but also some level of inter-form compatibility (sometimes with lower seed set). Although some forms showed more selfcompatibility than others, none was completely selfing . Given that genomic groups in Geleznowia retain some level of outcrossing and genomic diversity, genetic cohesion across large distances (hundreds of kilometres in some cases; see Fig. 8) without mixing with overlapping and sometimes co-occurring groups seems unlikely to be solely due to selfing or the apomictic spread of clones. However, even if the patterns are driven or reinforced by selfing, there is an arguably effective barrier to gene flow leading to relative evolutionary independence among the groups, compatible with calling them species under the concept we are using.
Other examples of co-occurrence and retained genomic differences are found south of Kalbarri, where Groups 1 and 2 are found growing together at the same location (P15 and P16). Maintenance of morphological distinction between the two at the same location suggests distinction, even when some individuals show evidence of introgression (see below). Further south, P9 co-occurs with Group 1 (P10)    without good evidence for current admixture or morphological ambiguity. And finally, Groups 1 and 3 flower at different times at the same location (P5; Group 3 sampled the first year, both observed the second year). The field observations of co-occurrence between Groups 1 and 3 and the broad overlap in their distributions (see Fig. 8), while maintaining genomic distinctiveness point to reproductive isolation. Differentiating species boundaries from population divergence becomes more difficult when taxa are allopatric (Rannala and Yang 2020), and this was largely the case for our sampling of Groups 4 and 5. These taxa are morphologically similar and were considered one entity during sampling (G. sp. Marchagee). Although specimen distributions may partly overlap (Group 4 has a patchy and limited distribution, but P3 is ~3 km from an unsampled population of plants with morphology matching Group 5), we did not observe them co-occurring at the same location. Genomic evidence from the PCA, SplitsTree4 network and Structure analysis indicates a close relationship but with consistent differences. The phylogenetic analyses suggest that they may be closely related, but support values for this were variable (43% in the concatenation tree, 93% in the lineages tree and 74% in the species tree) and were potentially complicated by the inclusion of the divergent populations P4 and P9. Given that the two groups are not widely geographically separated on the basis of field observations and herbarium specimens and that there is no strong evidence for admixture, coupled with consistent morphological differences (see Taxonomic treatment), we recognise them as separate species. Future work may help map the distributions of the two more clearly, broaden sampling and refine understanding of their relationship.  Our analyses point to complex patterns of divergence for populations P4 and P9. They are found on isolated branches in the distance network and concatenation tree but show indications of admixture in the Structure analysis and in the conflict and concordance factors in the lineages tree that could possibly represent ancient hybrid origins. As mentioned above, P9 is morphologically and genomically distinct from co-occurring Group 1 and can be recognised as a distinct species. However, the P4 population is less morphologically distinct, with characters similar to Group 3 (see Table 2). Since it is genomically divergent but associated with Group 4, it is unclear how to treat the population. Given Groups 3, 4 and 5 occur in the region and are geographically close to P4, a hybrid origin seems possible, although perhaps not recently enough to show clear genomic admixture. We provisionally retain the population as a single genomically divergent lineage of possible hybrid origin between ancestors of Groups 3 and 4.

Taxonomic implications
On the basis of agreement between genomic and morphological differences that are typically maintained in sympatry, we recognise seven species in Geleznowia. We retain the currently circumscribed G. amabilis and provide new circumscriptions for the six remaining species. For clarity, a summary is provided in Table 3, so as to clearly link the genomic groups used within this paper to previous identifications and the new taxonomy.
It became evident through the examination of highresolution scans of type specimens recently made available via Global Plants (http://plants.jstor.org/), particularly of the holotype of G. verrucosa (KW 001001059) held in the M.G. Kholodny Institute of Botany, National Academy of Sciences of Ukraine, that the current concept of G. verrucosa had been misapplied in Western Australia. The type instead matched the morphology of Genomic group 5, originally treated in WA as the phrase-named G. sp. Marchagee. Specimens previously considered to be typical G. verrucosa are recognised here as two different and distinct species. The first, corresponding to Genomic group 1, represents one of Bentham's previously recognised species, G. calycina, which is reinstated here (and includes his G. macrocarpa as a synonym). The second, corresponding to Genomic group 3, is newly described here as G. uberiflora K.A.Sheph. & A.D.Crawford.

Evolution in the group
One pattern evident in both the SplitsTree4 network and the Structure analysis was admixture between some individuals of Groups 1 and 2 in the co-occurring populations of P15 and P16. The location of these populations matches that of a putative hybrid population from earlier work by , who found genetic and morphological evidence of hybridisation between their putative taxa. At the time,  considered Geleznowia to consist of just two subspecies, so matching their concepts to ours is imprecise. Regardless, we also found hybridisation in that population, though we interpret our results at the species rather than subspecies level (see below). Our results show strong evidence for current admixture between groups 1 and 2, with one individual of group 2 having a genomic composition consistent with an F1 or early-generation hybrid. Other individuals showed evidence of some introgressed genomic background, including the individual of group 1.  found evidence that introgression in the hybrid population was asymmetrical. Although we found some indication of introgression in both directions, more admixture was detected in individuals from Group 2 (e.g. Fig. 4), in line with asymmetric introgression. Although the admixture was evident in our genomic results, individuals were readily assigned to either morphological group without ambiguity, suggesting that the morphological signal of hybridisation was weak. This is consistent with previous observations of only a few morphological intermediates and the absence of an obvious 'hybrid zone' .
The presence of some rare admixture suggests that reproductive isolation between Groups 1 and 2 is not complete. This is not incompatible with species divergence,  Harrison and Larson 2014). Broadhurst (2000) suggested that their two putative taxa were distinct, but because there were additional 'intermediate' form(s) possibly derived from hybridisation, the two should be recognised as subspecies (however, this was not formally done). The topic of using the rank of subspecies and its evolutionary relevance has been and continues to be contentious (de Queiroz 2020(de Queiroz , 2021Hillis 2021;Burbrink et al. 2022). It has been suggested that the rank of subspecies may be used to denote species that are incompletely separated (de Queiroz 2020). In practice, evidence of genomic mixing between divergent lineages has been used to justify keeping them as single species (e.g. Georges et al. 2018), whereas lack of clear evolutionary independence has led to suggestions to sink some morphologically recognised subspecies (e.g. Prates et al. 2023). The choice remains partly subjective with regards to how much divergence or evidence of evolutionary independence is sufficient to recognise taxa. We choose to recognise Groups 1 and 2 at the level of species, despite the evidence of admixture, for a number of reasons, including the following: (1) morphological divergence between the two groups is strong, without indication that it is breaking down in the population where the two co-occur and undergo introgression; (2) observation of admixture is rare and highly localised (no admixture in the widespread Group 1 was clearly evident outside of the hybrid population); and (3) the degree of genomic divergence between them is similar to that between other genomic groups that co-occur without admixture (e.g. Groups 1 and 3) that we have interpreted as separate species. The P4 population showed indications of potential hybrid origin in the Structure and SVDquartets analyses, which may help explain why it is morphologically less distinct but genomically divergent. In the SVDquartets lineages tree (Fig. 7), there is indication that a substantial proportion of sites are in conflict and point to affinities elsewhere, although it is supported (80%) as sister to Group 4. Morphologically, these plants are closest to those in Group 3. The flowers of the P4 plants have a number of bracts similar to that of the flowers of Group 3 and a lack of bracteoles similar to some Group 3 plants (and Group 4 plants). Flower size is more similar to flowers in Group 3 than flowers in Group 4. On the basis of the SplitsTree4 network, there is no indication of recent introgression, suggesting that the ambiguity may be the result of an ancient hybridisation event. The P4 population is geographically close to the ranges of both Group 3 and Group 4, so a hybrid origin is possible. Further sampling in the region around P4 may help clarify whether there are any other similarly anomalous populations.
Our phylogenetic results indicated little support for treelike relationships among species of Geleznowia, with lower support and conflict along the backbone of the group. One possible explanation for this pattern is a rapid radiation and the fixation of ancestral variation in small populations. In line with previous work by Broadhurst et al. (1999), we detected extremely low genetic diversity within and sometimes among populations (results not shown), in some cases resulting in effective clonality. Given the ephemeral nature of populations of Geleznowia and their response to disturbance, population bottlenecks could have substantial effects on their genetic diversity (Broadhurst et al. 1999). The presence of essentially identical genomic samples is also consistent with previous breeding system work, indicating facultative selfing in the group, which varied in degree between putative taxa . Diversification within Geleznowia may be dominated by these kinds of processes, where morphological and genetic variation is rapidly fixed and conversely rapidly lost in small populations.

Conservation implications
Our proposed new taxonomy of Geleznowia highlights species of immediate conservation concern, because some of these species are highly range restricted and under threat from land clearing. For example, G. eximia (P9) and G. occulta (Group 4) are known from only one or two locations and are therefore in urgent need of further survey to determine whether there are additional populations. The situation is particularly dire for G. eximia, which is currently known from only two plants in a single population and has not been previously conservation listed given the lack of taxonomic recognition. It will be listed as a Priority One conservation species (T. Llorens, pers. comm.) but is likely to warrant consideration for assessment as Threatened. Geleznowia amabilis (Group 6) remains known from a small area near Kalbarri, and G. narcissoides (Group 2) is currently difficult to locate and in need of survey. Herbarium records suggest G. narcissoides was potentially more widespread; however, it is unclear whether all recorded populations are persisting. Given G. narcissoides was a primary target for the cut flower industry, review of the conservation management of this species is needed.
The observation of low genomic diversity (unexpected clonality) for most species of Geleznowia (all except G. amabilis) is notable in that it applies to both localised and widespread species. Conservation assessment and management would be improved by a more extensive assessment of clonality and diversity within species of Geleznowia. Given these species maintain a soil seed bank and appear to depend on disturbance such as fire to regenerate populations, populations are often transient. This complicates survey efforts to assess diversity and full geographic distributions and means that the genetic diversity represented by the soil seed bank is largely unmeasured. The provision of formal species names and their descriptions will greatly aid survey efforts and allow conservation assessments for these poorly known species. Erect subshrubs to shrubs 0.15-2 m high, branchlets terete, glandular-verrucose, glabrous or with an indumentum of simple hairs. Leaves sessile or with a short petiole, overlapping and crowded towards terminal branches, coriaceous, elliptic to obovate, glandular-verrucose adaxial surface, with or without short, simple hairs. Flowers 1-17 in terminal heads. Bracts 0-14 petaloid, yellow, surrounding flowers. Sepals 5, free, 3.5-14 mm long, imbricate and resembling bracts. Petals 5, yellow, elliptic, thicker than sepals, 4-8.8 mm long. Stamens 10, free, 2.5-5.5 mm long, glabrous. Carpels 5, free, thickened at the apex, glandular-verrucose, with two ovules per carpel; style glabrous, with a narrow to club-shaped stigma. Fruit obovoid with a single seed per locule. Seeds dark brown to black 3.6-5 mm long, with a pale aril.

Distribution
A genus of seven species endemic to Western Australia (Fig. 8).

Key to the species of Geleznowia
www.publish.csiro.au/sb Australian Systematic Botany

Distribution and habitat
Currently known only from a few populations in or near Kalbarri National Park (Fig. 8) in the Geraldton Sandplains bioregion (Department of Agriculture, Water and the Environment 2020). This species is found growing in yellow or brown sand over sandstone, or red-brown sandy loam with laterite, in coastal scrubland, dense heath, low mallee or Acacia shrubland with Calytrix, Grevillea, Calothamnus and Melaleuca.

Phenology
This species flowers from July to October, with fruits forming in October to November. The distinctive vivid goldenyellow colour of the bracts and sepals is maintained throughout flowering, although some outer bracts rarely become tinged with red towards the apex as fruits develop.

Conservation status
This species is listed as Priority Two under Conservation Codes for Western Australian Flora. Although some populations are found within a National Park, the extent of the distribution of this species remains poorly known and further survey is required.

Etymology
From the Latin amabilis, meaning worthy of love.

Distribution and habitat
Widespread through the northern sandplains from west of Gillingarra to south of Shark Bay (Fig. 8)

Phenology
Early flowering species with flowers at anthesis from May to early August, and fruits forming in late August to September. The yellow petaloid bracts usually become distinctly tinged with red post-pollination (Fig. 10c).

Conservation status
Reasonably widespread and not considered to be under threat at this time.

Typification
In April 1854, William Harvey met James Drummond in Perth and they discussed 'several new and curious genera www.publish.csiro.au/sb Australian Systematic Botany   (Harvey 1854a, p. 110). Drummond had briefly described this new species, without providing a name, as 'a stiff upright-growing shrub, about 2 ft [~0.61 m] high; the flowers are borne in corymbs from nine inches to a foot in diameter; they are not very conspicuous, but are accompanied by numerous large bracts of a golden-yellow colour, which render this one of the most showy of our native plants; it appears on sand-plains to the east and west of the southern branch of the Hill River, and in other similar situations to the south of the Irwin River' (Drummond 1853, p. 122 Of the available syntypes for Geleznowia macrocarpa, Wilson treated the Kew specimen from Hooker's Herbarium (K 000717328) as the 'holotype'. This specimen is here designated as the lectotype for Geleznowia macrocarpa. This sheet is an Oldfield collection from the Murchison River with a 'Herbarium Hookerianum 1867' stamp and was therefore available to Bentham. Similarly, the MEL 232723 specimen is an Oldfield collection from the same locality and is initialled with Bentham's characteristic 'B', so it was also seen by him. The fragment on the MEL 232721 specimen also has a note in Bentham's hand, indicating that the material is from the specimen sent to Hooker and as such is treated here as an isolectotype, even though there is no information to suggest it was an Oldfield collection. Overall, the quality of the type material is very poor, consisting of bare branchlets with a few loose leaves and flowers with the subtending bracts having fallen off, thus the total number of flowers or bracts within each inflorescence cannot be confirmed. Distinguishing features mentioned by Bentham (1863) in his protologue include narrow sepals and the '[c]occi (not yet fully ripe) more than twice as long as broad'. There are a few specimens at PERTH that have sepals narrower than typical (3.2-4.2 mm wide cf. 4.5-8.5(9.6) mm wide) and an ovary that is longer than broad (e.g. C.A. Gardner 7744 PERTH 0968617), corresponding to the flowers observed on the type. These specimens are intact and are currently identified as G. calycina because of the size and number of flowers and bracts within each inflorescence. Therefore, G. macrocarpa is currently retained herein as a synonym of G. calycina. It should be noted that within this 'narrow sepal' group there are a few specimens that have quite dense hairs and two specimens where the fruits also have simple hairs and unusual finger-like protrusions extending from the verrucose glands at the apex (e.g. A.M. Ashby 2172 PERTH 00969206; R. Bates 4052 PERTH 04274695).

Notes
This species can be distinguished from other Geleznowia by the following character combination: a low subshrub 0.3-0.75 mm high with dark green leaves; usually with 5-7 (sometimes 1-4) flowers per inflorescence, surrounded by (7)10-14 pale yellow to yellow petaloid bracts (often tinged red post-pollination), 8.2-16.9 mm long, 5-10.3 mm wide, abaxial surface glabrous or sometimes with scattered hairs 0.05-0.1 mm long; 8-12 bracteoles; sepals 9-12.3 mm long, 4.7-8 mm wide; and a broad stigma 0.4-0.5 mm wide. Geleznowia calycina most closely resembles G. uberiflora but www.publish.csiro.au/sb Australian Systematic Botany can be distinguished from it by its maximum flower number of 5-7 flowers per inflorescence (cf. usually 1-3 flowers), with more bracts ((7)10-14 cf. (5)6-8) and more bracteoles (8-12 cf. 0-6). G. calycina is also an early flowering species with flowers appearing as early as May and fruits forming in late August to September (cf. peak flowering in August-September in G. uberiflora). Moreover, the bracts and sepals in G. calycina usually turn red post-pollination (Fig. 10c), which readily distinguishes this species when it is found growing with the later-flowering G. uberiflora (Fig. 10e). Geleznowia calycina is also somewhat morphologically similar to G. narcissoides but is readily distinguished by its fewer flowers, with 5-7 flowers per inflorescence (cf. (5) 7-10 flowers), greenish-yellow to yellow petaloid bracts (cf. pale lemon yellow), and bracts, bracteoles and sepals being glabrous or with an indumentum of scattered to moderately dense hairs up to 0.1 mm long (cf. an indumentum of moderately dense to dense hairs up to 1.2 mm long). G. calycina also flowers earlier than G. narcissoides, often showing the distinctive red blush to bracts postpollination when co-occurring with G. calycina when its flowers are in bud or reaching anthesis in spring.

Distribution and habitat
Known only from one location (Fig. 8). Found growing in deep yellow sand in low shrubland with Grevillea, Acacia, Hibbertia and Conospermum.

Phenology
Flowering specimens collected in September and October, with fruits forming in late October.

Conservation status
This species is currently known only from a single population with two extant plants. As such, it is likely to meet criteria for listing as Threatened (Critically Endangered); however, further survey is urgently needed to confirm this. It is to be listed as Priority One under Conservation Codes for Western Australian Flora (T. Llorens, pers. comm.).

Etymology
From the Latin eximius, meaning exceptional or uncommon, in reference to the attractiveness and rarity of this new species.

Distribution and habitat
This species is known from a few widespread populations from north of Geraldton to north of Kalbarri (Fig. 8) in the Geraldton Sandplains bioregion (Department of Agriculture, Water and the Environment 2020). However, recent field work has failed to relocate this species in some areas where it has been previously collected, and it is unclear how many populations are currently persisting. Found growing on gentle slopes or flats in white-grey or yellow sand or brown sand over laterite in Banksia woodlands with Grevillea, Hibbertia, and small-flowered myrtaceous species or in low heath with Acacia, Banksia attenuata, Allocasuarina, Grevillea, Calothamnus, Calytrix and Stirlingia.

Etymology
From Narcissus L. and the Latin -oides (like), alluding to the showy flowers of this species being reminiscent of a doubleheaded daffodil.

Distribution and habitat
This species is currently known from only two populations west of Coorow and east of Goomalling (Fig. 8)

Conservation status
This species is currently known from only two populations, one of which is on private property. The second population, although occurring in a nature reserve, has been observed only in low numbers. It is to be listed as Priority Two under Conservation Codes for Western Australian Flora (T. Llorens, pers. comm.).

Etymology
From the Latin occultus (hidden, concealed) in reference to the cryptic nature of this species, because it was initially recognised as distinct through this molecular study, and that it is currently known from only two localities.

Distribution and habitat
This species is reasonably widespread through the northern sandplains from south-east of Goomalling and west of Dalwallinu to the Kalbarri region (Fig. 8)

Phenology
Peak flowering from August to September with fruits beginning to form in late September through to October.

Conservation status
This species is widespread and represented within the conservation estate and therefore is not currently considered to be under threat.

Etymology
From the Latin uber (fruitful, abundant) and -florus (-flowered), in reference to the numerous showy yellow flowers typical of this species.
One population from south of Eneabba (K.A. Shepherd & C.F. Wilkins KS 1712 (PERTH 09507817)) was detected as genomically divergent (Group P4) and is provisionally identified here as a hybrid G. uberiflora × occulta.

Distribution and habitat
This widespread species occurs from east of Watheroo to Dirk Hartog Island (Fig. 8) in the Geraldton Sandplains and Yalgoo bioregions (Department of Agriculture, Water and the Environment 2020). Found growing in flat or gently sloping areas, sometimes in disturbed sites such as road verges or fire breaks, in white to brown or yellow sand, sometimes over laterite or sandstone, in low, dense heath or shrublands associated with Allocasuarina, Acacia, Banksia, Grevillea and Callitris.

Phenology
Flowering commences in July and extends through to early spring with fruits forming in September to October.

Conservation status
This species is widespread and not considered to be under threat at this time.

Typification
The single specimen lodged at the M.G. Kholodny Institute of Botany, National Academy of Sciences of Ukraine (KW 001001059) was used by Turczaninow to write the protologue and is therefore recognised as the holotype (Mosyakin et al. 2019). The bottom two fragments on the specimen at the Natural History Museum London (BM 001015536) labelled 'J. Drummond 3: 8' (pencil marked 'a') represent isotype material, whereas the top three E. Pritzel 644 fragments (pencil marked 'b') do not. Five fragments mounted on the top and bottom left (K 000717329) of a sheet stamped 'Herbarium Hookerianum 1867' and two fragments (K 000717330) mounted on blue paper affixed to the bottom right and stamped 'Herbarium Benthamianum 1845' are both labelled as 'Drummond n. 8' collections and are recognised as isotypes. On the bottom right of the MEL 232716 sheet is a blue Mueller label with 'G. verrucosa Turcz, W.A. J.Dr. 8' in his handwriting, above are two '8' tags and a packet of fragmented material mounted on the top right that is a match for the protologue of G. verrucosa. As such, this material is also treated as an isotype. The mounted specimen on the left-hand side of this MEL 232716 sheet has a Drummond '83' tag and is an isolectotype of G. calycina (see typification section under that name). The MEL 232725 sheet has a blue Mueller label with 'Geleznowia verrucosa Turcz. W.A. J.Dr.' in his handwriting but a note in pencil in another hand at the bottom says 'evidently, no 8'. Paul Wilson noted that this specimen was a 'probable isotype' (in sched., 21/9/1999); however, here it is treated as a possible isotype.