Australian biogeography, climate-dependent diversification and phylogenomics of the spectacular Chamelaucieae tribe (Myrtaceae)

Taxon sampling

We sampled 131 taxa representing all but one (97%, 36/37) of the genera in Chamelaucieae, and 10 of the 11 subtribes sensu Rye et al. (2020), omitting only Baeckea L. sens. str. in Baeckeinae Schauer (Table 1). The sampling strategy within Chamelaucieae attempted to maximise sampling of generic diversity, utilising available published and unpublished ribosomal and chloroplast phylogenies. Where possible, we attempted to sample at least two species per currently accepted genus, including representatives of the two most divergent clades and a representative of each formally recognised section. Representatives of all putative clades were sampled for putatively polyphyletic and paraphyletic genera (Baeckea sens. lat., Darwinia Rudge, Homoranthus A.Cunn. ex Schauer, Malleostemon J.W.Green, Verticordia DC.). We attempted to sample all morphologically orphaned groups of Baeckea that remain unplaced following the revisions of Rye and co-authors (mostly belonging to subtribe Hysterobaeckineae Rye & Peter G.Wilson; Rye et al. 2020). We therefore aimed to capture a near-complete sampling of Chamelaucieae to approximately sectional level. We also attempted to sample from the geographic extent of each genus to maximise representation for biogeographic signal; however, species-level phylogenies are unavailable for many genera, limiting assessment of whether any phylogenetically or biogeographically important clades were unsampled. We also sampled six outgroup taxa in the tribes Osbornieae (Osbornia F.Muell.), Lindsayomyrteae (Lindsayomyrtus B.Hyland & Steenis), Eucalypteae (Eucalyptus L’Hér.) and the sister tribe Leptospermeae (2× Leptospermum J.R.Forst. & G.Forst. and 1× Pericalymma (Endl.) Endl.; Wilson and Heslewood 2023). We supplemented the dataset with 35 samples from the Maurin et al. (2021) Myrtales dataset generated as part of the Plant and Fungal Trees of Life initiative (PAFTOL; www.paftol.org).

DNA extractions were performed using QIAGEN DNEasy Plant Minikits (Hilden, Germany). All other laboratory work was completed at the Australian Genome Research Facility, Melbourne. Genomic DNA extracts were fragmented enzymatically to ~350 base pairs, except in samples that were already sheared through natural degradation. DNA libraries were prepared with NEBNextUltra II FS Library Prep Kit (New England Biolabs, Ipswich, MA, USA). Hybridisation capture was carried out on pools of 12–16 libraries, enriched using the Angiosperms353 v1 bait kit (Johnson et al. 2019) with V5 chemistry; probes were produced by myBaits target capture kit (Daicel Arbor Biosciences, Ann Arbor, MI, USA; catalogue number #308108.v5). After hybridisation capture, enriched libraries were sequenced on Illumina NovaSeq. 6000 SP (Illumina Inc., San Diego, CA, USA) with v1.5 chemistry and generating 150-bp paired-end reads.

Bioinformatic data processing

High-throughput 150-bp paired-end and combined unpaired reads were imported into CLC Genomics Workbench (ver. 20.0.2, see https://digitalinsights.qiagen.com/) for trimming using a quality score limit of 0.05 (Phred score ~13) and discarding reads with length <50 nucleotides. The cleaned reads for each individual were randomly sampled to 4 million reads, to reduce computational burden for downstream processing steps.

The cleaned reads were used to recover the Angiosperms353 target sequences using HybPiper (ver. 2.1.1, see https://github.com/mossmatters/HybPiper; Johnson et al. 2016) with a target file comprising sequences of Leptospermum erubescens Schauer, Pericalymma spongiocaule Cranfield and Eucalyptus grandis W.Hill, sourced from the Kew Tree of Life Explorer (Baker et al. 2021). The HybPiper ‘assemble’ command was run using BLASTx (ver. 2.14.0; Camacho et al. 2009) to map reads to the translated target sequences and the mapped reads were subsequently assembled de novo using SPAdes (ver. 4.0.0, see https://github.com/ablab/spades; Bankevich et al. 2012). Coding sequences (hereafter exons) were extracted using exonerate (Slater and Birney 2005), and the regions flanking the exons (i.e. introns and untranslated regions, UTRs) were recovered using the ‘--run_intronerate’ flag and combined with the exons to form supercontig sequences. The ‘paralog_retriever’ command was used to recover putatively paralogous exon sequences.

We used the Paragone (ver. 1.0.0, see https://github.com/chrisjackson-pellicle/ParaGone) pipeline to process the HybPiper output (Jackson et al. 2023) that provides an implementation of the orthology inference approach developed by Yang and Smith (2014). Following multiple sequence alignment, Paragone uses HmmCleaner (ver. 0.180750, 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, see http://trimal.cgenomics.org/; Capella-Gutiérrez et al. 2009) to remove sequencing and alignment errors, and rogue taxa from multiple sequence alignments. We used the monophyletic outgroups (MO) algorithm to recover putatively orthologous gene alignments for the exon data using an internal outgroup file comprising Leptospermum erubescens, Leptospermum polygalifolium Salisb. and Pericalymma spongiocaule. Paragone was used to generate cleaned sequence alignments (i.e. all the pre-paralogy resolution subcommands were run) for the supercontig data. We used the MAFFT (ver. 7.471, see https://mafft.cbrc.jp/alignment/software/; Katoh and Standley 2013) ‘mafft_auto’ algorithm, gene tree inference was performed using FastTree 2 (ver. 2.1.1, see http://www.microbesonline.org/fasttree; Price et al. 2010) and TreeShrink was run with a q value of 0.2 for the alignment steps. Alignments with fewer than 20 sequences were removed from downstream analyses.

Phylogenetic analyses

The exon and supercontig gene alignments output from Paragone were each concatenated using Geneious Prime (ver. 2022.0.1, see https://www.geneious.com) and we used IQ-TREE (ver. 2.2.0, http://www.iqtree.org/; Nguyen et al. 2015; Chernomor et al. 2016) to generate maximum likelihood (ML) phylogenies. ModelFinder (ver. 1.0, see http://www.iqtree.org/ModelFinder/; Kalyaanamoorthy et al. 2017) was used to estimate the best partitioning scheme and models (–MFP + MERGE) with the relaxed hierarchical clustering algorithm (Lanfear et al. 2014) restricting the search to the top 10% of partition merging schemes (–rclust 10) for the concatenated data. Branch support was estimated using the best partitioning scheme with 1000 ultrafast bootstrap replicates (Minh et al. 2013) and the IQ-TREE implementation of the approximate likelihood ratio test (aLRT; Guindon et al. 2010), with interpretation of support following Minh et al. (2013) and recommendations in the IQ-TREE manual (see http://www.iqtree.org/doc/Frequently-Asked-Questions).

In addition to the concatenated analyses, we estimated coalescent-based species trees for the exon and supercontig data sets using weighted ASTRAL (wASTRAL; Zhang and Mirarab 2022) as implemented in Accurate Species Tree Estimator (ASTER, ver. 1.19, see https://github.com/chaoszhang/ASTER). We inferred gene trees for each locus using IQ-TREE with support estimated using the approximate Bayes test (Anisimova et al. 2011). The wASTRAL algorithm attempts to account for gene tree uncertainty in the species tree estimates by weighting to reduce the impact of long terminal branches and low branch support. We followed the authors’ recommendations (see https://github.com/chaoszhang/ASTER/blob/master/tutorial/wastral.md) for upper and lower thresholds for local Bayesian support (1.0 and 0.33 respectively) and used local posterior probability (LPP) to measure internal branch support for the species tree estimate.

Both the concatenated and coalescent based species trees were further interrogated using quartet sampling (QS; Pease et al. 2018). The QS approach subsamples quartets from the target tree and alignment to produce four metrics: Quartet Concordance (QC), Quartet Differential (QD), Quartet Informativeness (QI) and Quartet Fidelity (QF). Taken together, QS scores describe the degree of topological variation in the data (Pease et al. 2018), providing a valuable alternative to measures such as bootstrap that can be a poor measure of variation in the data, particularly among large data sets (Lanfear and Hahn 2024). We ran QS analyses with 500 replicates and a log-likelihood cut-off of 2.

We focused on the supercontig concatenated dataset for our phylogenetic dating, biogeographic and diversification analyses as this is more phylogenetically informative (double the total alignment length) than the exon dataset. We used the concatenated as opposed to the coalescent topology as our QS analyses indicated that low QC support along the backbone is due to low phylogenetic informativeness of individual genes rather than significant conflict or skewed discordance (Pease et al. 2018). Inaccurate or poorly resolved gene trees may limit the reconstruction of a more accurate species tree than those from concatenation approaches (Gatesy and Springer 2014). Furthermore, studies have shown that concatenation approaches can perform equally well or even better than coalescent approaches that attempt to account for incomplete lineage sorting (Tonini et al. 2015).

Phylogenetic dating

We used BEAST (ver. 2.4.7, see https://www.beast2.org/; Bouckaert et al. 2014) for divergence age estimates. Although there are pollen microfossils of Chamelaucieae (Thornhill et al. 2012, 2015), these cannot be confidently assigned to genera or groups within the tribe due to homoplasy and a lack of diagnostic characters. We therefore used the Chamelaucieae + Leptospermeae crown age (uniform prior: 47.8–103 Ma) from Thornhill et al. (2015), based on an Eocene Myrtaceidites Cookson & K.M.Pike fossil (Thornhill and Macphail 2012), as a secondary calibration in the BEAST analyses, following Nge et al. (2022). The maximum age constraint for our analyses was set to 103 Ma for Myrtaceae based on the stem age of the family from the BEAST analyses of an angiosperm-wide study with 238 fossil calibrations (Ramírez-Barahona et al. 2020).

We used a gene-shopping approach using the SortaDate pipeline (see https://github.com/FePhyFoFum/sortadate; Smith et al. 2018) to account for computational limitations and achieve adequate convergence (>200 effective sample size; ESS) from the BEAST analyses. We selected a subset of the most clock-like loci of the 300+ nuclear markers using SortaDate (root-to-tip variance filtering option). We tested for the effects of different subset sizes (10, 20 and 30 loci) on divergence age estimates. As the effects of subset sizes are minimal, we used the 30 most clock-like loci as the main dataset for this study, as a balance between including more loci while achieving computationally feasible convergence in BEAST (see Supplementary Table S1). We constrained the starting topology using a chronogram obtained from treePL (ver. 1.0, see https://github.com/blackrim/treePL; Smith and O’Meara 2012). The input topology for treePL was obtained using IQ-TREE. The treePL topology was calibrated using the same secondary calibration points for our BEAST analysis, allowing the runs to initialise successfully, as the specified calibration nodes are required to overlap with the respective nodes in the starting topology. We also conducted sensitivity tests for the dating analyses. We tested the effects of genetic locus sampling on our dating analyses by conducting BEAST analyses on subsets of the 10, 20 and 30 most clock-like loci from the gene-shopping approach. As divergence age estimates were not significantly different for these subsets (Supplementary Fig. S1 and Table S1), we used the 30-gene subset for subsequent biogeographic and diversification analyses.

The configuration XML file for BEAST was created using BEAUti (ver. 2.4.7, see https://beast.community/beauti), with the 30 most clock-like loci concatenated as the input alignment and a constrained treePL starting topology provided as a Newick tree file. The GTR + G substitution model was applied for the entire alignment (Abadi et al. 2019). The relaxed log-normal clock model and birth–death tree priors were chosen in BEAUti. Secondary calibrations were implemented as uniform priors, with the 95% confidence intervals (CI) from Thornhill et al. (2015) specified as minimum and maximum bounds. Three independent BEAST runs were performed with 50 million Markov Chain Monte Carlo (MCMC) chains each, to ensure convergence towards a global instead of a local optimum, sampling every 1000th generation. Convergence (>200 ESS) for each run was assessed using Tracer (ver. 1.6, A. Rambaut, M. A. Suchard, D. Xie and A. J. Drummond, see https://github.com/beast-dev/tracer/releases/tag/v1.6) with the log files as input. The first 20% MCMC chains for each run were discarded as burn-in. The output tree files from the three runs were subsequently combined using LogCombiner (ver. 2.4.7, see https://beast.community/logcombiner) and annotated in TreeAnnotator (ver. 2.4.7, see https://www.beast2.org/treeannotator/). The output dated phylogenies were visualised in FigTree (ver. 1.4.3, A. Rambaut see https://github.com/rambaut/figtree/releases) and final figures edited in Inkscape (ver. 1.2, see https://inkscape.org, accessed April 2024). All outgroups were pruned from the dated BEAST phylogeny using the ‘drop.tip’ function in the R package ape (ver. 5.3, see http://ape-package.ird.fr/; Paradis et al. 2004) prior to biogeography and diversification analyses, as we focused on the crown diversification of the tribe in this study. Investigating the stem diversification of Chamelaucieae would require a more extensive and complete sampling of sister outgroups than is in scope for this paper.

Biogeographic analyses

The geographic distribution of each sequenced taxon was determined using vouchered herbarium records obtained from the Australasian Virtual Herbarium (see http://avh.chah.org.au, accessed March 2023) and cross-validated with Florabase for recent taxonomic changes of Western Australian taxa (see https://florabase.dbca.wa.gov.au, accessed April 2023). Cultivated and erroneous records (e.g. outliers, records in the ocean) were excluded. We delimited six geographic regions that largely correspond to present-day Australian biomes (Table 2): south-west Western Australia (W; Mediterranean climate), arid Eremaean interior (A), eastern Australia (E; temperate to subtropical), northern Australia including the Kimberley and Top End regions (N; monsoonal), Cape York peninsula and north-eastern Australia (C; monsoonal, wet tropical and subtropical) and the semi-arid transition zone between W and A (T). Region boundaries were based on the Interim Biogeographic Regionalisation of Australia (IBRA; Environment Australia 2000), following previous studies with some modifications (see Supplementary Table S3; Cook et al. 2015; Nge et al. 2022, 2023). As many Chamelaucieae endemic or near-endemic taxa occur around the semi-arid periphery of south-west Western Australia (i.e. part of the Eremaean but not strictly in the arid interior), we delimited the two adjacent IBRA regions (Yalgoo and Coolgardie) as a separate region (T). We separated regions N and C across the Gulf of Carpentaria for northern Australia, as this has been demonstrated to be a significant biogeographic barrier for many taxa (Bowman et al. 2010). The New Caledonian Sannantha virgata (J.R.Forst. & G.Forst.) Peter G.Wilson was coded as occurring in eastern Australia, to limit the number of biogeographic regions due to computational constraints and because the study focuses on biogeographic disjunctions within Australia.

Biogeographic analyses conducted using BioGeoBEARS (ver. 1.1, N. J. Matzke, see https://github.com/nmatzke/Matzke_R_binaries; Matzke 2013) in R (ver. 3.5.1, R Foundation for Statistical Computing, Vienna, Austria, see https://www.r-project.org/) inferred relative probabilities of ancestral areas for each node in the dated phylogeny. We included variations of the dispersal–extinction–cladogenesis (DEC; Ree and Smith 2008) and dispersal–vicariance (DIVA; Ronquist 1997) models. We excluded BAYAREA models (Landis et al. 2013) as these do not account for vicariance (Pontes-Nogueira et al. 2021). As in previous studies, we included additional model parameters such as jump dispersal (j), dispersal as a function of geographic distance (x) and time stratified parameters that included information on how respective regions across Australia have shifted through time (Nge et al. 2022, 2023). We implemented a scenario of progressive aridification and climate cooling since the Oligocene (32 Ma) by progressively larger geographic dispersal distances between regions (W+T, N+C, E), with an expanding arid centre as a barrier for dispersal (see Supplementary Table S4). Results from the biogeographic analyses were compared using Akaike Information Criterion scores (AIC; Akaike 1974) in the BioGeoBEARS package to determine the best model for our dataset.

Homalocalyx F.Muell., Ochrosperma Trudgen, and Calytrix were underrepresented geographically compared with other genera in our study. To account for this, we conducted two separate BioGeoBEARS analyses with different geographic codings for these genera. The first analysis used actual geographic coding of the sampled species. The second took into account results from a previously published study by Nge et al. (2022) with near-complete taxonomic sampling and complete geographic representation for Calytrix, and unpublished data for Homalocalyx and Ochrosperma (with complete sampling of all species for both) that inferred widespread ancestral areas for these genera. Accordingly, the second coding gave Calytrix, Homalocalyx and Ochrosperma widespread distributions. Results from both coding regimes were similar and we focus on the second analysis as we believe this is a more accurate representation of the biogeographic history of these genera (see the Supplementary material for details, Fig. S2). We also conducted a sensitivity analysis for BioGeoBEARS with the wASTRAL (supercontig) topology instead of the IQ-TREE topology, dated using treePL instead of BEAST or starBEAST due to computational constraints. This analysis also used the widespread geographic distribution coding for Homalocalyx, Ochrosperma and Calytrix.

We estimated the number and directions of dispersal–vicariant patterns to test H2 (asymmetrical dispersal patterns) in BioGeoBEARS using biogeographic stochastic mapping (BSM; Matzke 2015). Event frequencies were taken from the mean of event counts across 50 BSMs.

Diversification analyses

We applied the paleoenvironmental birth–death models of Condamine et al. (2013) using RPANDA (ver. 1.8, see https://cran.r-project.org/package=RPANDA; Morlon et al. 2016) in R to test for climate-dependent diversification patterns (H3). We compared pure birth (no extinction), constant birth–death, and time-dependent diversification models with other birth–death models that incorporated global temperature and sea-level changes. Variations in speciation and extinction rates (constant v. exponential) were also included as part of the model testing (Table 3), resulting in 14 different models. We chose the best-fitting model as the one with the highest AIC weight (out of a value of 1).

We used Bayesian analyses of macroevolutionary mixtures (BAMM; Rabosky 2014) on our dated BEAST phylogeny to test for diversification rate-shifts (H4). Three sampling fractions were specified to test the effects of incomplete taxonomic sampling: (1) the sampling fraction for each subtribe was specified for BAMM (Supplementary Table S5), (2) a global sampling fraction for the entire Chamelaucieae tribe was specified based on the total number of described species from Rye et al. (2020) (18.8%, 124/660 spp.) and (3) a more conservative global sampling fraction for the tribe based on the total number of described and phased-named taxa from Rye et al. (2020) (15%, 124/810 spp.). The BAMM analyses were run for 5 million and 20 million MCMC generations for the global sampling and subtribal sampling fractions respectively, to achieve convergence (>200 ESS). Convergence of each run was assessed using BAMMtools (ver. 2.1.10, see https://cran.r-project.org/package=BAMMtools; Rabosky et al. 2014) in R. The first 10% of BAMM runs was discarded as burn-in. Credibility diversification rate shifts were also determined using BAMMtools, with two different expected rate shifts as our priors: (1) one expected rate shift event as recommended for smaller phylogenies with <200 tips, and (2) three expected rate shifts to account for the wider unsampled diversity of the tribe (660+ species).

We also used GeoHiSSE (ver. 1.9.6; Caetano et al. 2018) in R to compare diversification rates between south-west Western Australia and other regions (i.e. geographic-dependent diversification). GeoHiSSE implements 36 models to test geographic-dependent speciation or extinction and hidden states (Supplementary Table S6). As GeoHiSSE only analyses binary traits, the geographic coding used in our biogeographic analyses were collapsed to two regions (W and non-W). All taxa with occurrence records in W (even non-endemics) were coded as ‘W’. AIC scores were used to determine the best-fitting model.

Data set construction

The average number of reads that were mapped to target sequences using HybPiper was ~340,000 comprising ~10% of our subsampled sequence reads (Supplementary Table S2). In general, gene recovery was high across samples averaging ~310 genes with sequences and ~255 genes recovered with at least 50% of the target length. After paralog filtering, 295 loci were retrieved from the MO option in Paragone. The final exon and supercontig alignments are 155,922 bp (292 loci) and 367,371 (341 loci) in length respectively, including 131 taxa.

Phylogenomic relationships and divergence time of Chamelaucieae

We use the following terms to discuss support values in our phylogenetic trees: (1) nodes with aLRT = 100, BS = 100, LPP = 1, QC/QD/QI = 1/NA/1 are fully supported, (2) nodes with aLRT > 80, BS > 95 LPP > 0.85, QC > 0.5 are strongly supported, (3) nodes with aLRT 80, BS 95, LPP 0.85–0.6, QC 0.3–0.5 are moderately supported, and (4) nodes with aLRT < 80, BS < 95, LPP < 0.6, QC < 0.3 are weakly supported. Support values for nodes across Chamelaucieae are generally higher for the supercontigs than the exons dataset. In general, we recovered backbone topologies with strongly supported values from aLRT/BS/LPP and weaker support from QS.

Several notable differences in topology and support values were detected across different analyses (concatenated v. coalescent) and from different datasets (exon v. supercontig) (Fig. 2, 3). Despite these differences, there are several consistent findings from our results. Chamelaucieae was recovered as monophyletic, sister to Leptospermeae with full support in all analyses across both datasets. Ochrospermatinae is recovered as the sister taxon to all other Chamelaucieae subtribes, albeit with weak support across analyses in both datasets. Rinziinae sens. lat. is non-monophyletic with strong to moderate support across all analyses and both datasets.

Fig. 2.

Phylogenetic relationships of subtribes in Chamelaucieae based on concatenated (IQ-TREE) and coalescent (ASTRAL) approaches for supercontig and exon datasets. (a) IQ-TREE and (b) ASTRAL topologies based on the supercontig dataset. (c) IQ-TREE and (d) ASTRAL topologies based on the exons dataset. Support values are based on approximate likelihood ratio tests (aLRT)/ultrafast bootstrap (BS) for IQ-TREE and local posterior probability (LPP) for ASTRAL. Nodes have full support (100/100/1; aLRT/BS/LPP) unless otherwise highlighted with red circles, and support values shown. Rinziinae* includes Enekbatus and Rinzia, and Rinziinae** includes Euryomyrtus.

Fig. 3.

Maximum likelihood IQ-TREE phylogeny of Chamelaucieae based on the supercontig dataset. Support values are approximate likelihood ratio tests (aLRT) and ultrafast bootstrap (BS). Nodes with full support (100/100, aLRT/BS) are indicated with black circles. Nodes without full support are highlighted as red circles, with actual support values shown. Cladogram depicts phylogenetic relationships across subtribes, and Clade A, Clade B.

All subtribes except for Rinziinae sens. lat. are recovered as monophyletic with strong support based on concatenated and coalescent support values. Support for these relationships was lower based on the quartet sampling approach, with most subtribes (6/11) recovered as monophyletic with strong support from the exon dataset. This increased to eight subtribes with strong support based on the supercontig dataset (with monophyly of Micromyrtinae and Thryptomeninae respectively now strongly supported). The remaining subtribes are recovered as monophyletic with moderate to weak support from quartet sampling.

Clade A, consisting of five subtribes (Chamelauciinae, Calytricinae, Alutinae, Micromyrtinae and Thryptomeninae), had moderate to strong support across analyses and the two datasets (Fig. 2). Clade B (Hysterobaeckeinae and Rinziinae sens. str.) was recovered with strong support from the supercontig dataset but moderate to weak support from the exon dataset, across both concatenated and coalescent approaches. The phylogenetic placements of Stenostegiinae and Astarteinae are uncertain. The placement of Stenostegiinae alternates from being either sister to the remainder of Clade A (supercontig IQ-TREE, exon ASTRAL, exon IQ-TREE) or sister to Astarteinae (supercontig ASTRAL), with high to full support from aLRT/BS/LPP and weak support from QC. Astarteinae either forms a grade with Clade A (supercontig dataset) or is sister to Clade B (exon dataset), with strong to moderate support from aLRT/BS/LPP and weak support from QC. Different topologies were recovered from the exon dataset depending on the analysis method. The concatenated approach showed Rinziinae sens. str. to be non-monophyletic, with Euryomyrtus Schauer forming a grade with Astarteinae, the rest of Rinziinae sens. str. (Enekbatus Trudgen & Rye, Rinzia Schauer) and Hysterobaeckeinae. By contrast, the coalescent approach showed Astarteinae forming a grade with Clade A (Rinziinae sens. str. as monophyletic and Hysterobaeckeinae).

The topological placement of subtribe Micromyrtinae also alternates within Clade A, with differing topologies from concatenated and coalescent approaches across both datasets. Resulting IQ-TREEs showed that the subtribe is recovered as sister to Alutinae + Calytricinae from the supercontig dataset and sister to Calytricinae only from the exon dataset, with weak support (aLRT/BS/LPP/QC). However, in coalescent ASTRAL trees, Micromyrtinae is sister to Chamelauciinae with weak to strong (LPP = 0.98 from supercontig dataset) support.

Several genera in Chamelaucieae including Homoranthus, Verticordia, Darwinia, Baeckea sens. lat., Malleostemon and Hysterobaeckea (Nied.) Rye were recovered as non-monophyletic with full to strong support based on aLRT/BS/LPP but with low QC values.

The stem age of Chamelaucieae was estimated at 63.5 Ma (95% confidence interval; CI: 47.8–88.8 Ma) and the crown radiation at 42.4 Ma (95% CI: 31.2–61.1 Ma; Fig. 4). All subtribes diverged and radiated in the Oligocene to early Miocene (20–33 Ma; Fig. 4). The New Caledonian Sannantha virgata diverged from sister Australian relatives at c. 6.7 Ma (95% CI: 4.1–9.7 Ma).

Fig. 4.

Maximum likelihood ancestral range estimation of Chamelaucieae (Myrtaceae) based on the best model DIVA + j + x with time stratification (dotted lines) from BioGeoBEARS. Phylogenetic relationships shown are based on the IQ-TREE supercontig topology. Dated phylogeny was constructed with secondary calibrations implemented in BEAST, with 30 top clock-like loci subset from the Angio353 bait kit. Pie diagrams at nodes show relative probability of ancestral states, with common area combinations highlighted in colours and the remainder in grey. See Supplementary Table S7 for all possible area combinations included. The map shows region delimitation across the Australian continent based on IBRA regions: south-west Western Australia (W), semi-arid adjacent transition zone (T), arid centre (A), eastern Australia (E), Cape York and north-eastern Australia (C), and northern Australia (N) that includes the Kimberley and Top End. The red line highlights the Eocene–Oligocene boundary (33 Ma). Subtribes are highlighted in different colours. Sannantha virgata is a New Caledonian endemic.

Biogeography

The best-fitting biogeographic model with the lowest AIC score was DIVA + j + x, with the addition of aridification scenario time stratified distance matrices (Supplementary Table S7). The ancestral area for Chamelaucieae was unresolved, with near equal probability of ancestral area occurring across multiple regions of the continent, supporting H1 (Fig. 1, Supplementary Table S8, Fig. S13). However, the ancestral area of Chamelaucieae excluding subtribe Ochrospermatinae (sister to remainder of the tribe) was estimated to be south-west Western Australia (region W; Probability, P = 0.99, node126). Similarly, the ancestral area of the deeper backbone divergences was predominantly estimated to be W (Fig. 2). Furthermore, 5 of the 11 sampled subtribes had the ancestral area estimated to be W: Micromyrtinae, Chamelauciinae, Astarteinae, Rinziinae and Hysterobaeckeinae. The ancestral area for Alutinae + Calytricinae + Micromyrtinae was estimated to be in central Australia (region A; P = 0.58, node 234). The biogeographic results from our wASTRAL topology yielded largely similar results, of widespread ancestral distribution for Chamelaucieae and W as the ancestral area for the deeper backbone divergences (Supplementary Fig. S3). A notable difference between our main results (IQ-TREE topology) and the sensitivity analysis (wASTRAL topology) is the ancestral area of Stenostegiinae and Astarteinae. This was recovered as region N from the wASTRAL topology but W from the IQ-TREE topology due to different placements of the two subtribes in the respective phylogenies (Fig. 4, S3).

Directionality and type of biogeographic events

Most biogeographic events comprised dispersals (49%, 104/213) followed by within-area speciation (42%, 89/213) (Table 4). Very few vicariance events (10%, 20/213) were inferred. Range expansion dispersals were more common (42%, 90/213) than founder events (7%, 14/213).

The direction of dispersal events was found to be strongly asymmetrical, in support of H2 and dominated by dispersals from south-west Western Australia to the semi-arid transition zone (T; 22% of all dispersal events) and arid interior (A; 12%) (Fig. 5, Table 5). Multiple independent dispersal events from W to T occurred in Chamelaucieae including Verticordia sens. lat., Balaustion Hook., Cheyniana Rye and Malleostemon (Fig. 4). Overall, south-west Western Australia was the largest source area of dispersals (49% of all dispersal events) and the semi-arid transition zone the largest sink (33%). The central arid region served as an important connecting area for other parts of Australia, with dispersals recorded outwards to all other regions (Table 5). Dispersals out of northern (N) and eastern (E) Australia were limited, confined to only dispersals into Cape York and north-eastern Australia (C). Dispersal modes also differed between regions: dispersals from W to T were largely range expansions (21%), with very few long-distance founder events (1%). By contrast, dispersal rates from W to A were more similar among range expansions (8%) and founder events (3%). Not surprisingly, range expansions occurred more frequently between adjacent areas (e.g. W and T).

Fig. 5.

Summary of dispersal events (range expansion and founder events) and associated directionality for Chamelaucieae across Australia, estimated with biogeographic stochastic mapping in BioGeoBEARS (Table 3). Size of arrows represents frequency of dispersal events, red being the most common (W → T). Dashed arrows are the least common dispersal events and events with a frequency of less than 1 are not shown. See Table 5 for frequencies of dispersal events for each area. Letters represent different regions across Australia based on IBRA regions: south-west Western Australia (W), semi-arid adjacent transition zone (T), arid centre (A), eastern Australia (E), Cape York and north-eastern Australia (C) and northern Australia (N) that includes the Kimberley and Top End.

Temperature-dependent speciation and diversification rate shifts in Chamelaucieae

Temperature-dependent diversification models were overwhelmingly the best-fitting models in the RPANDA analyses, supporting H3 (Table 3). The model with the highest AIC weight (0.66) has speciation varying exponentially with temperature without extinction. In addition, the second- and third-best models were also temperature-dependent, with speciation varying exponentially with temperature and extinction constant (AIC weight 0.23), and both speciation and extinction varying exponentially with temperature (AIC weight 0.08). Speciation rates declined from the initial crown radiation in the Eocene (c. 40 Ma) until after the E-O boundary at 34 Ma (Fig. 6). The speciation rate remained relatively constant thereafter until an increase during the Middle Miocene Climatic Optimum (MMCO; 14–17.5 Ma) and subsequent decline towards the present (Fig. 6). A declining speciation trend towards the present was also found in our BAMM analyses across all different sampling fraction and rate shift regimes (Fig. 6, Supplementary Fig. S4–S7).

Fig. 6.

Diversification dynamics of Chamelaucieae. (a) Temperature-dependent paleoenvironmental model showing speciation decline towards the present. The Eocene–Oligocene boundary is marked with a red column. (b) Declining speciation rates shown from BAMM analyses. (c) The most credible diversification rate shift configuration from BAMM showing a rate shift at the crown of Chamelaucieae, with sampling fractions specified for each subtribe and an expected rate shift prior of 1.

One diversification rate increase was detected, at the crown radiation of Chamelauciinae in the late Oligocene to Miocene (c. 26.85 Ma, 95% CI: 19.1–37.98 Ma), by the BAMM analyses with subsampling fractions specified for each subtribe (Fig. 6, Supplementary Fig. S8–S12). The increased radiation of this subtribe initially occurred in south-west Western Australia based on our BioGeoBEARS results, with subsequent dispersal and smaller radiations in eastern Australia and into the arid interior (Fig. 4). By contrast, no rate shifts were detected when we used a global sampling fraction (both actual and conservative regimes; Fig. S8–S12 and Table S9). Our GeoHiSSE results indicate that model 11 (weighed AIC = 0.5/1) was the best-fitting model, followed by model 17 (weighed AIC = 0.38/1; Table S5). Taken together, these results support H4 (lack of diversification rate shifts in south-west Western Australia).

Phylogenetic relationships

Our study presents different topologies for Chamelaucieae depending on the dataset used (exon v. supercontig) and phylogenetic method applied (coalescent or concatenated), with differing levels of support depending on the metric used (see results above). Coalescent and concatenated approaches have increasingly been shown to potentially yield differing topologies even from the same phylogenomic dataset (e.g. Antonelli et al. 2021; Helmstetter et al. 2024, 2025; Nge et al. 2024b; Zuntini et al. 2024). Similarly, differing support values may be obtained from different metrics. Our topologies generally have high aLRT/BS (concatenated) and LPP (coalescent) support values, in contrast to low QC values (indicative of weak majority but not conflicting support for alternate topologies). Although having high LPP support in contrast to low QC values might appear counterintuitive, this is not unusual as these are calculated differently (Sayyari and Mirarab 2016; Pease et al. 2018). The LPP is recommended as the general measure of support for the species topology given that this indicates the probability of a selected node being the true node given the set of gene trees available, assuming the presence of incomplete lineage sorting (ILS; Sayyari and Mirarab 2016). Low QC values along the backbone topologies of Chamelaucieae along with high QD and QI values are likely due to a rapid radiation in the early evolutionary history of the tribe – due to poor phylogenetic informativeness (low QC) likely due to ILS, rather than discordant skew or conflict (low QD) or poor genetic sampling (low QI) (Pease et al. 2018). Ultimately, each different methodological approach has inherent strengths and limitations, and should be assessed critically for each individual case (Joyce et al. 2024, in press). Nevertheless, we argue that these discordances in our topologies should not significantly affect our biogeographic results for several reasons. All subtribes apart from Ochrospermatinae have higher species richness in W, resulting in all backbone nodes apart from the crown node having W as ancestral. The unresolved ancestral area of the Chamelaucieae crown group would likely remain true as Ochrosperma is consistently recovered as the earliest diverging lineage in all phylogenomic analyses, albeit with lower support from the coalescent analyses. All analyses recovered similar topologies within subtribes that were strongly supported as monophyletic (except Rinziinae sens. lat.).

Our phylogenomic trees have some notable differences from previous results based largely on Sanger plastid and rDNA ETS sequence data. All subtribes were recovered as monophyletic except for Rinziinae sens. lat., in which Astus Trudgen & Rye and Triplarina Raf. are separated from Rinziinae sens. str. (Rinzia, Enekbatus and Euryomyrtus), a result foreshadowed by Rye et al. (2020) based on weak support for the combined clade, compared to strong support for each of the two subclades. Ochrosperma was recovered as sister to the remaining Chamelaucieae in all cases. Ochrosperma has the ancestral anther type (dorsifixed, versatile) and what is possibly the ancestral seed type (reniform) that may justify the position as sister to the remainder of the tribe. However, Calytricinae is not sister to Ochrosperma in our study, in contrast to Rye et al. (2020). We inferred Calytricinae as sister to Alutinae in all analyses. In Rye et al. (2020), Alutinae and Micromyrtinae were sister taxa, in contrast to our study in which both tribes were not recovered as sister in every analysis. These incongruences between nuclear (our study) and plastid + ribosomal ETS topologies from previous studies suggest the presence of non-bifurcating evolution in the early evolutionary history of the tribe, likely attributable to incomplete lineage sorting, hybridisation or a combination of both (Folk et al. 2017; Nge et al. 2021a; McLay et al. 2023). Low or negative QC and low QD values noted for Homoranthus and Darwinia also indicate the presence of reticulation among shallower nodes, supporting previous claims of recent hybridisation between extant species of Chamelaucieae (Briggs 1964; Rye et al. 2013; Rye 2022).

Baeckea sens. lat. is non-monophyletic in our topology, confirming the findings of others based on plastid data (Lam et al. 2002; Rye et al. 2020). Baeckea sens. str. is the only genus missing from the tribal sampling. Previous phylogenies did not resolve the position, e.g. within a large polytomy of subtribes in Rye et al. (2020). All of the orphaned ‘Baeckea’ species belong to Hysterobaeckeinae. Some of these orphaned Baeckea species require evaluation as to whether inclusion in existing sister genera is possible or the description of new genera is required: B. elderiana E.Pritz. (sister to Kardomia Peter G.Wilson), B. muricata C.A.Gardner (sister to Hysterobaeckea on a long branch), B. pentagonantha F.Muell. + B. sp. Dudawa (sister to the Kardomia–B. elderiana clade) and B. sp. London Bridge (sister to Cheyniana + Oxymyrrhine Schauer). Other orphaned Baeckea species are discussed below.

Micromyrtus flava (J.W.Green) Rye & Peter G.Wilson is recovered nested within other species of the genus, supporting the recent recombination from the previous monotypic genus Corynanthera J.W.Green (Rye and Wilson 2022).

Scholtzia Schauer is paraphyletic with respect to Babingtonia Lindl. and several orphaned Baeckea species due to the presence of S. sp. Geraldton amongst the grade of Baeckea orphans leading to Babingtonia. A number of orphaned Baeckea species (B. robusta F.Muell., B. staminosa E.Pritz., B. subcuneata F.Muell., B. sp. Walkaway and B. sp. Walyahmoning) are phylogenetically placed in the subtribe Hysterobaeckeinae (rather than Baeckeinae) in the Scholtzia clade, with Babingtonia apparently nested within this, and require further investigation to establish whether the latter two genera require recircumscription or whether additional genera should be erected for the atypical species. Baeckea staminosa in particular may be best placed in Scholtzia as this was resolved sister to two core Scholtzia samples, and has 2-lobed anthers (a common state in Scholtzia) and 2 collateral ovules per loculus, (1 or 2 superposed ovules per loculus in Scholtzia) in comparison with 4 or more ovules per loculus in other members of the Hysterobaeckeinae.

Malleostemon is monophyletic with the removal of two species that are morphologically atypical. Malleostemon sp. Adelong falls sister to Kardomia + Baeckea elderiana, requiring further assessment of generic limits in this clade, potentially requiring a new genus. Malleostemon sp. Officer Basin falls within Hysterobaeckea, consistent with its morphology (Rye 2016) and the species should be described in this genus (B. L. Rye, in prep.).

The genera Darwinia, Homoranthus and Verticordia were also recovered as non-monophyletic in our study, and this is supported by densely sampled ribosomal and plastid data (M. D. Barrett, unpubl. data). In particular, the three Western Australian Darwinia species (D. diosmoides (DC.) Benth., D. meeboldii C.A.Gardner and D. speciosa (Meisn.) Benth.) were paraphyletic with Actinodium Schauer nested within, and remote from the eastern Australian species (D. fascicularis Rudge and D. micropetala (F.Muell.) Benth.) that grouped with Homoranthus, a relationship previously noted by Lam et al. (2002). The core of Homoranthus is monophyletic but the genus is rendered polyphyletic by the somewhat divergent position of the aberrant H. homoranthoides (F.Muell.) Craven & S.R.Jones that has previously been placed in a monotypic genus Schuermannia F.Muell. (Mueller 1853). Our study is the first to include multiple Verticordia species and the genus is resolved as grossly polyphyletic with the genera Actinodium, Chamelaucium, Darwinia, Homoranthus and Pileanthus Labill. nested within Verticordia. Detailed taxonomic revision and new circumscriptions for subtribe Chamelauciinae are in progress and will be published in due course.

The polymorphic genus Rinzia was comprehensively sampled for this study with seven species representing all five sections, including type species of four of the sections (R. carnosa (S.Moore) Trudgen, type of sect. Mesostemon Rye; R. fumana Schauer, type of the autonymic sect. Rinzia; R. longifolia Turcz., type of sect. Discolora Rye; and R. polystemonea, type of sect. Polyandra Rye). We were unable to sample R. orientalis Rye, the type of sect. Semasperma Rye (Rye 2017) but the section is represented by two other species, R. ericaea (F.Muell. ex Benth.) Rye and R. icosandra (F.Muell. ex Benth.) Rye. Despite the morphological diversity within Rinzia, the genus is confirmed as monophyletic, an important conclusion given previous suggestions that more than one genus should be recognised (Rye 2017). However, a uniting morphological apomorphy for Rinzia is not known and the genus is defined within subtribe Rinziinae by exclusion from the more uniform Euryomyrtus and Enekbatus.

The widespread, variable genus Hysterobaeckea is resolved as monophyletic, with the inclusion of Malleostemon sp. Officer Basin in that genus (B. L. Rye, in prep.). The most-divergent species in the genus is Baeckea aff. behrii, restricted to a small patch of granite on Mount Ziel in central Australia, the highest peak west of the Great Dividing Range and disjunct from the remainder of Hysterobaeckea to the south. This disjunction could possibly have resulted from an early vicariance event as central Australia dried, prior to the genus Hysterobaeckea reinvading the arid zone from the south. Hysterobaeckea is an unusual example of a primarily arid zone genus, occurring from south-eastern Australia across the north of the Nullarbor Plain and in arid margins of south-western Australia.

Species with 1-locular, indehiscent fruits, i.e. nuts, with a terminal style (the old tribe Chamelaucieae sensu Johnson and Briggs 1984, i.e. the clade comprising subtribes Alutinae, Calytricinae, Chamelauciinae, Micromyrtinae and Thryptomeninae here) are resolved as monophyletic with molecular data for the first time. Atypical nuts with the style inserted in a depression occur in Malleostemon (all species), Scholtzia (in 2 species only, most species having a multi-locular, indehiscent fruit), 1–2 species of Astartea DC., Hysterobaeckea and Babingtonia (the 3 latter genera usually with a multi-locular dehiscent fruit, i.e. a capsule). In all these taxa the presence of styles inserted in pits at the nut apex hints at a more recent ancestry from capsules.

All subtribes are supported by morphological characters, as provided below for the new subtribe Triplarininae and more narrowly circumscribed subtribe Rinziinae, and as previously provided (Rye et al. 2020) for the other 10 subtribes.

Cheyniana, Oxymyrrhine and Baeckea sp. London Bridge all have very reduced ±globular anthers, lacking an enlarged connective gland, in agreement with the positioning within a single clade. Baeckea sp. London Bridge may require recognition in a distinct genus as this orphaned species has unusual dark reddish flowers that may help support such a conclusion.

Bird pollination arose independently at least five times within Chamelaucieae (Keighery 1982; Brundrett et al. 2024), with a radiation of bird-pollinated bell inflorescences (at least 35 species in Darwinia; subtribe Chamelauciinae), best known for the Mountain Bells in the Stirling Range of south-west Western Australia. The stem divergence of the Mountain Bell clade (represented by Darwinia meeboldii in our study) from Actinodium occurred in the Miocene c. 13 Ma (95% CI: 8.9–18.6 Ma) but dozens of putatively bird-pollinated Darwinia species were not sampled, preventing reconstruction of bird pollinating syndromes in this genus. Other presumably bird-pollinated species occur sporadically in Verticordia sensu lato (13 species suggested to be bird-pollinated by Brundrett et al. 2024, e.g. the sampled V. humilis Benth. and several other species not sampled e.g. V. grandis J.Drumm.) and Chamelaucium (at least four unsampled species; Brundrett et al. 2024). In other subtribes of Chamelaucieae, only two independent origins of bird-pollinated lineages are noted in the subtribe Hysterobaeckeinae, with only one extant species per lineage – Balaustion pulcherrimum Hook. and Cheyniana microphylla (C.A.Gardner) Rye (Rye 2009). In view of the striking radiation of the Mountain Bells, the lack of radiations in some other bird-pollinated Chamelaucieae is puzzling and warrants further investigation.

Widespread ancestral area and continental vicariance

The ancestral area for Chamelaucieae is unresolved in our analyses, with equal probabilities of different area combinations. Most of these (80%+) have more than one region, suggesting that the ancestral area for the tribe was widespread across Australia historically, thus supporting the Peripheral Vicariance Hypothesis (PVH; H1). Early divergences within the tribe were vicariant events across the northern, eastern and south-western parts of Australia: Ochrosperma (regions E/N), Stenostegia A.R.Bean (N) and the remaining Chamelaucieae (W as ancestral). Our results therefore support an emerging pattern of peripheral vicariance recovered from multiple Australian mesic plant groups (Nge et al. 2022), suggesting that PVH is a dominant biogeographic signal in the Australian flora. Suggestions or explicit recovery of this pattern have been found for Nicotiana L. (Ladiges et al. 2012), Boronia Sm. (Duretto and Ladiges 1998) and Rhamnaceae (Ladiges et al. 2005; Nge et al. 2021b, 2023).

Interestingly, except for the ancestral node of Chamelaucieae, all nodes in the backbone had south-west Western Australia (W) as ancestral. Other biogeographic studies of Australian plants have also shown W to be ancestral (Mast and Givnish 2002; Jabaily et al. 2014; Lamont et al. 2016; Nauheimer et al. 2018). However, Ladiges et al. (2012) and McLay et al. (2016) cautioned that these findings may be artefactual and this was supported by biogeographic sensitivity tests of Nge et al. (2022). W was estimated to be ancestral for Calytrix (Myrtaceae) when default biogeographic models were applied in BioGeoBEARS, without the addition of dispersal distance parameters and time stratification but when these were added the ancestral area was widespread (Nge et al. 2022). However, our study has W as ancestral for most Chamelaucieae clades even with the addition of distance and time stratified parameters in the models. A large proportion of the tribe could well have originated in W, congruent with an earlier presence and subsequent dominance of sclerophyllous elements replacing rainforest in that region since the Eocene (Hill 1998, 2017; Hill and Brodribb 2001; Crisp and Cook 2013). Alternatively, severe extinction of Chamelaucieae lineages in other parts of Australia may have erased vicariant signals for this group (Sanmartín and Meseguer 2016). Indeed, our paleoenvironmental diversification results indicate that there have likely been significant extinctions (or at the very least, a decrease in speciation) in Chamelaucieae during the Eocene–Oligocene. Testing for the effects of extinction as a parameter in biogeographic models is worth further investigation, e.g. with the application of PhyBears (N. J. Matzke, unpubl. data).

We found that most clades in Chamelaucieae diversified within the respective geographic regions (within-area speciation), with few dispersal events between regions. This has also been shown from previous studies on Calytrix (Nge et al. 2022), Cryptandra (Nge et al. 2023), Pomaderris (Nge et al. 2021b), Banksia L.f. (Cardillo and Pratt 2013) and Goodeniaceae (Jabaily et al. 2014). Clades that correspond to different geographic regions of Australia with little biogeographic exchange is a common pattern across the Australian flora – a pattern that can even be seen from phylogenetic studies without an explicit biogeographic analytical component (e.g. de Lange et al. 2010; Clowes et al. 2022; Duretto et al. 2023; Orel et al. 2023; Wilson and Heslewood 2023). Vicariance in southern Australia has been attributed to the uplift of the calcareous Nullarbor Plain in the Miocene (14 Ma; Crisp and Cook 2007). This edaphic barrier limits dispersal between the mesic regions of southern Australia, resulting in more within-area speciation rather than dispersal events across the barrier. However, the Nullarbor uplift is clearly not the only driver for these geographic patterns, as vicariances were not restricted to the Miocene. For example, we infer a divergence between eastern Australian Homoranthus and south-west Western Australia Verticordia sens. str. at c. 22.8 Ma (95% CI: 16.4–32.3 Ma), preceding the uplift of the Nullarbor at 14 Ma. We believe that peripheral vicariance, argued above and in previous studies, is a more convincing explanation.

Asymmetrical dispersal into arid Australia and speciation decline

Chamelaucieae has a surprisingly high proportion of dispersals amongst all biogeographic events (47%; Table 4). In a separate study of Calytrix within the tribe, dispersals accounted for only 21% of events (Nge et al. 2022), as was found in a global study of Solanaceae (20% in Dupin et al. 2017). The high proportion of dispersals in Chamelaucieae is largely driven by range expansions rather than founder events, mainly from W to the adjacent T (Table 5). We show that transitions from W to the semi-arid (T) and arid (A) interior are common in Chamelaucieae, supporting H2. Similar patterns have also been noted in Hakea Schrad. & J.C.Wendl. (Cardillo et al. 2017), Triodia R.Br. (Toon et al. 2015), Lomandra Labill. (Gunn et al. 2020, 2024), Schoenus L. (Elliott et al. 2021), Goodeniaceae (Jabaily et al. 2014) and Cryptandra (Nge et al. 2023) but not for Pomaderris (Nge et al. 2021b). There are three likely reasons for the dispersal bias. Firstly, many lineages in the mediterranean W are dry adapted and may therefore be pre-adapted to colonise the drier habitats of T and A (Cardillo et al. 2017), whereas many lineages in eastern Australia inhabit wetter temperate and subtropical areas. However, few transitions from mediterranean southern South Australia (as part of eastern Australia) into the arid zone have been recovered. Hopper (1979) suggested that the transition zone of W is a fertile area for recent allopatric speciation due to greater fragmentation of habitats and populations from past climatic change than in the arid zone and wetter areas of W. Secondly, high species and lineage diversity in W may simply mean that the region has more potential colonising lineages to disperse into adjacent arid regions than regions with more depauperate floras. Thirdly, West Wind Drift may bias dispersals from west to east in continental Australia and the wider Southern Hemisphere (Sanmartín et al. 2007; Nauheimer et al. 2018). Dispersal towards the east is also evident in our phylogenetic placement of Sannantha Peter G.Wilson, with the New Caledonian endemic (S. virgata) nested in an eastern Australian clade.

Our study provides counterevidence to the hypothesis of Australian arid flora being primarily derived from dispersal of coastal floras that were pre-adapted to xeric conditions (Burbidge 1960; Schodde 1989; Crisp et al. 1999; Ladiges et al. 2011). Although this may be true for some taxa, as hypothesised for Amaranthaceae and Chenopodiaceae (Kadereit et al. 2005; Cabrera et al. 2011), we show that arid zone lineages in Chamelaucieae are mainly derived from adjacent regions, in Hysterobaeckeinae, Chamelauciinae, Rinziinae sens. str. and Calytricinae (Nge et al. 2022) but appear to be absent in Astarteinae. Species in Hypocalymma (Endl.) Endl. and Astartea (Astarteinae) are largely restricted to the mesic fringes of south-west Western Australia, where most inhabit swampy habitats (Rye 2013; Rye et al. 2013; Keighery et al. 2023). Another genus from this subtribe, Cyathostemon Turcz., does reach inland well into T and almost reaches to the border with A. This subtribe is likely limited by the requirement for wetter niches, either from transitions of more dry-adapted ancestors similar to most other present-day Chamelaucieae or persistence in a wetter niche from the past.

Despite frequent transitions from south-west Western Australia to more arid areas, these did not often result in radiations in the arid zone. Only one clade (Hysterobaeckea) shows convincing evidence of an arid (semi-arid) radiation in the Miocene (c. 20 Ma). We acknowledge that more extensive sampling may be required to further test this pattern for other potential clades, e.g. in Rinziinae and Pileanthus–Chamelaucium–Verticordia subg. Eperephes A.S.George. However, this is an emerging pattern in multiple sclerophyllous mesic Australian plant groups (Cardillo and Pratt 2013; Cardillo et al. 2017; Elliott et al. 2021; Habib et al. 2022; Nge et al. 2022; Gunn et al. 2024). The previous dominant narrative of sclerophyllous plant radiations coinciding with aridification episodes in the Miocene (Crisp et al. 2004; Crayn et al. 2006) should be revisited. As a result of denser species-level sampling, we know that multiple transitions from mesic regions into the arid zone have occurred for plants in Australia but these lineages often fail to radiate, with only a few examples of arid radiations having been documented (Kadereit et al. 2005; Cabrera et al. 2011; Dodsworth et al. 2020; Renner et al. 2020; Čalasan et al. 2022). We also show a declining trend of speciation with increasing aridification towards the present, similar to Cryptandra (Nge et al. 2023), Acacia Mill. (Renner et al. 2020) and a wide range of other sclerophyllous flora in Australia (F. J. Nge, unpubl. data). Even arid elements of the Australian biota (Hancock et al. 2018; Hammer et al. 2021; Tiatragul et al. 2023) show decreases in speciation towards the present, suggesting that the consequences for speciation during intensifying Pliocene and Pleistocene aridification (Martin 2006; Byrne et al. 2008; Byrne and Murphy 2020) have been severe. We detect an increase in speciation for Chamelaucieae during the relatively mild MMCO (15–17 Ma), followed by further speciation decline, accelerated during the Pliocene and coinciding with increasing aridification during that time (Martin 2006; Fujioka et al. 2005, 2009). Whether synchronous radiations of sclerophyllous elements in the Miocene across Australia are a response to increased seasonality, warming in the MMCO, gradual aridification or other factors, these warrant further testing using diversification models. However, we also note that the strongly asymmetrical sampling of the Chamelaucieae, weighted toward geographic representation and deeper divergences (genus or subtribal level), will strongly under-represent recent species-level divergence events in this study. A more complete species-level phylogeny may significantly alter conclusions about the effect of increasing aridity for Chamelaucieae.

Our study did not focus explicitly on biome evolution of Chamelaucieae, as incorporating the evolution of specific biomes into biogeographic models (e.g. Landis et al. 2021) requires additional simulations and information beyond the scope of this study. The initiation of the formation of Mediterranean climates in temperate Australia is unclear, although Beard (1977); Lamont and He (2017) suggested that Mediterranean climates first appeared in north-western Australia during the Oligocene (30 Ma) and subsequently moved southwards as the continent drifted north. The precise timing of the commencement of Australia’s northern monsoonal climate has also been contentious (Bowman et al. 2010).

Many clades in Chamelaucieae were already widespread prior to the likely initiation of mediterranean, monsoonal and arid biomes in Australia. Hence surmising that lineages would have transitioned and adapted to these newer biomes from older ones is reasonable. Present-day biomes largely correspond to the geographic regions in our study, with some minor differences. The mesic Murray Basin region of South Australia (SA) and parts of western Victoria (Vic.) have a mediterranean climate similar to south-west Western Australia. Despite this, southern Australia (SA + Vic.) is relatively depauperate in Chamelaucieae, with only a few species, e.g. Homoranthus homoranthoides, Rinzia ericaea, Micromyrtus ciliata (Sm.) Druce, Thryptomene micrantha J.D.Hook. and T. ericaea F.Muell present in the region. One of the two south-eastern taxa sequenced in our study (Homoranthus homoranthoides) shared affinities with other eastern Australian taxa that occupy non-Mediterranean biomes whereas Rinzia ericaea belonged in a clade with western species. Hence, some of the biogeographic connections of South Australian Mediterranean lineages are driven by geographic proximity with other mesic regions instead of biome fidelity with south-west Western Australia, a similar pattern found in other plant groups (Nge et al. 2021b, 2022, 2023; Gunn et al. 2024). Micromyrtus delicata A.R.Bean, found in the Wet Tropics region of Queensland, did not make a biome transition into rainforests, as this species is found in xeric open woodlands on sparse volcanic substrates (Bean 1997). In general, many species in Chamelaucieae are confined to skeletal oligotrophic substrates across a wide range of biomes. These include outlying taxa found in New Caledonia (Sannantha virgata) and Asia (species of Seorsus Rye & Trudgen; Rye and Trudgen 2008; Wilson and Heslewood 2014). Incorporation of biotic traits with biome evolution modelling (Bouchenak-Khelladi et al. 2015; Onstein et al. 2015; Onstein and Linder 2016; Landis et al. 2021) is a promising future avenue to test for interactions between these factors in the evolution of Chamelaucieae.

Eocene–Oligocene pulse event and lineage persistence in south-west Western Australia

We show that the crown radiation of Chamelaucieae pre-dates the E-O boundary and that there is a significant decline in speciation rate after the initial crown radiation until shortly after the E-O boundary pulse event, supporting H3 (Nge et al. 2020). Paleoenvironmental diversification models show that diversification of Chamelaucieae has been most affected by temperature changes, as opposed to other factors (e.g. time, sea-level). The tribe diversified more rapidly during a historically warmer climate in the Eocene and declined as the climate cooled towards the present. The rapid declining trend between the Eocene and the E-O boundary may be due to the effects of extinction during the E-O boundary pulse event pushing back age estimates of node divergence (Crisp and Cook 2009) or this may represent actual decline during that period. Untangling these two explanations will require further investigation. Our study differs from the findings of Nge et al. (2020) that showed an increase in diversification for Chamelaucieae during the E-O boundary. We believe that our current study, with more comprehensive sampling and a better-resolved dated phylogeny, offers a more accurate estimate of diversification in the tribe. Therefore we show that Chamelaucieae exhibits a similar diversification decline at the E-O boundary to other Australian plant groups (Nge et al. 2020). Crayn et al. (2006) showed that crown radiation of the temperate clade of Elaeocarpaceae (formerly Tremandraceae) occurred immediately after the E-O boundary, as has been documented for several Australian animal groups (Schweizer et al. 2011; Almeida et al. 2012; Kayaalp et al. 2013; Brennan and Oliver 2017; Marki et al. 2019). These studies reinforce the importance of the E-O boundary extinction pulse event as a significant driver for the evolution and diversification of our biota.

Our BAMM analyses detected only one diversification rate shift in south-west Western Australia, despite most extant species and lineages occurring in that region (supporting H4). This finding, along with our non-significant GeoHiSSE results, strongly indicates that high species diversity in south-west Western Australia is best attributed to more time for species accumulation rather than elevated diversification rates. The lack of diversification rate shifts is common in the south-west Western Australian flora (F. J. Nge, unpubl. data), although additional sampling of species-rich subtribes for Chamelaucieae may very likely reveal other rate shifts not detected in this study (Cusimano and Renner 2010). Diversification rate shifts in Australian plants have mostly been linked to trait adaptations to oligotrophic environments (Jordan et al. 2008; Carpenter et al. 2014; Onstein et al. 2016; Hayes et al. 2021), aridity (Renner et al. 2020, 2021; Hua et al. 2022) or fire (Crisp et al. 2011; Thornhill et al. 2019; Lamont et al. 2020) rather than being strictly geographically dependent.

Our dated phylogeny of Chamelaucieae also reveals several long-persisting lineages (isolated, long-branched, species-poor lineages; Nge et al. 2023). These include Ericomyrtus serpyllifolia (Turcz.) Rye (30.5 Ma; 95% CI: 23.0–45.1 Ma; three species of Ericomyrtus Turcz. not sampled), Aluta maisonneuvei (F.Muell.) Rye & Trudgen (33.0 Ma; 95% CI: 23.4–46.3 Ma; four other species of Aluta Rye & Trudgen not sampled) and Stenostegia congesta A.R.Bean (40.5 Ma, 95% CI: 29.6–57.8 Ma; monotypic, apparently ‘relictual’ and highly localised), although we note that, other than the monotypic Stenostegia A.R.Bean, the divergence age of these lineages may change significantly with additional ingroup sampling. Long-persisting lineages have also been detected in south-west Western Australian Cryptandra (Nge et al. 2023), Calytrix (Nge et al. 2022), Tremandra R.Br. ex DC. (Crayn et al. 2006), Franklandia R.Br. (Sauquet et al. 2009) and others (Hopper and Gioia 2004), and have been linked to the tectonic and climatic stability of the region (OCBILs; Hopper 2009; Hopper et al. 2021).

9. Myrtaceae tribe Chamelaucieae subtribe Rinziinae

Myrtaceae tribe Chamelaucieae subtribe Rinziinae Rye & Peter G.Wilson, Australian Systematic Botany 33: 203 (2020).

Type: Rinzia Schauer

Mostly prostrate or small shrubs to 1 m, rarely to 2.5 m tall, glabrous, sometimes with adventitious roots. Leaves opposite, decussate, small, usually shortly petiolate. Peduncles 1- or 2-flowered. Bracteoles usually persistent. Pedicels absent to distinctly longer than peduncles. Flowers 5-merous; hypanthium circular in TS. Sepals much shorter than petals. Petals white or pink. Staminodes usually absent, rarely numerous. Stamens 3–30, all antipetalous or opposite both the petals and sepals. Anthers dorsifixed, versatile or attached on front of broad filament; thecae parallel, longitudinally dehiscent; connective gland free. Ovary 2–5-locular; placentas ±sessile; ovules 2–10(–12) per loculus. Style base deeply inset in depression. Fruits inferior to largely superior, indehiscent (in Enekbatus) or dehiscent by 3 valves. Seeds 1–2.3 mm long, reniform or rarely obovoid, arillate or with a large flat hilum in most species (except when in an indehiscent fruit); testa crustaceous, usually thick (thin when fruit is indehiscent). x = 11.

Subtribe Rinziinae has 3 genera and 44 species, with most occurring in south-western Australia. This subtribe also occurs in south-eastern Australia, including Tasmania and has a single representative in the central mainland.