Three forms of distinctiveness – geographic, evolutionary, and functional – can be used to define priority species for conservation and diagnose the ecological and evolutionary drivers of biome assembly. Geographic distinctiveness defines the restrictedness of taxa in space; evolutionary distinctiveness defines the isolation of species through time; and functional distinctiveness defines the uniqueness of species traits. Whether these axes are coupled or decoupled across taxa and across geography at large spatial and phylogenetic scales is debated. Here, we use a large continental plant radiation (family Proteaceae, subfamily Grevilleoideae) as a case study to estimate these dimensions and investigate their coupling across diversity facets: from species-level to assemblage- and biome-scale diversity and endemism. We show that species-level distinctiveness along the three axes is largely decoupled across the radiation but positively correlated at the assemblage level. For example, the Tropical Forests biome harbours evolutionarily and functionally distinct paleoendemics, and species-rich Mediterranean and Temperate Forests biomes are assembled from functionally and phylogenetically related taxa. Some biomes, however, showed stark decoupling from this general trend, with Tropical Savannas containing phylogenetically related but functionally distinct taxa, with a mix of large and small geographic ranges. Our results suggest that multidimensional measures of distinctiveness are complementary but non-interchangeable. Different species sets contribute the most to each axis of distinctiveness, and these are filtered through assembly processes to shape the spectrum of incongruent to congruent spatial patterns across biomes, underscoring the need to consider multidimensional biodiversity patterns.

Spatial, Phylogenetic, and Functional Trait Data and Imputation

We obtained distributional ranges and phylogeny for 672 Australian Grevilleoideae species from a recent study (Skeels et al. 2025). Information on plant functional traits was obtained from AusTraits (Falster et al. 2021) for 651 Grevilleoideae species. The remaining 21 species in the phylogenetic and spatial dataset had no trait information. AusTraits contained 429 candidate traits, but coverage was uneven across Grevilleoideae species: on average, each trait was recorded for only 59 species. To minimise missing data, we retained 19 traits that were recorded for more than 80% of Grevilleoideae taxa. Six of these were invariant and were dropped, and the remaining 13 traits included flowering time, plant height, plant growth form, fruit length, fruit type, fruit dehiscence, leaf shape, leaf compoundness, leaf stomatal hairs, leaf length, leaf width, sex type, and diaspore fleshiness. These traits are related to key functional aspects of Grevilleoideae ecology, phenology, dispersal ability, and life history. From these traits, we further estimated six compound variables, for a total of 18 traits: we replaced the flowering time trait with four binary variables representing the presence or absence of summer, autumn, winter, and spring flowering, as well as the total months in flower per year. In addition to leaf width and leaf length as individual traits, we also estimated the ratio of leaf length to leaf width. Leaf shape was inconsistently classified, so we simplified this category by defining two states: needle leaves versus broad leaves, with needle leaves defined as those that were filiform, acicular, terete, or strap-shaped, and broad leaves as all other forms. We then excluded the 11 species with fewer than half of the 18 traits recorded, leaving 640 species with spatial, phylogenetic, and trait data, of which 498 had complete trait data. To test the effect of data missingness on our results, we replicated analyses on a reduced dataset of 498 species with no missing observations and on an imputed dataset of 640 species (details in the Supporting Information).

For the 640 species dataset, 170 trait values were missing out of 8,320 (~2%). To impute this missing data, we followed a similar protocol to Helmstetter et al. (2026). First, we log-transformed, centred, and scaled all quantitative values to unit variance. To account for phylogenetic non-independence, we decomposed the phylogenetic distance matrix containing all 640 species into a set of eigenvectors using PVRdecomp in the PVR package (Santos 2012) and used the first 10 eigenvectors as additional traits to train the imputation model. As we had mixed data types of continuous and categorical, we used the missForest function in the missForest package (Stekhoven and Bühlmann 2012) to impute traits based on a random forest algorithm. We repeated the procedure excluding the phylogenetic eigenvectors and compared the imputation accuracy with out-of-bag errors (the prediction error estimated from observations withheld from each bootstrap sample during random forest training) based on normalised root mean squared error (RMSE) for continuous traits and the proportion of falsely classified (PFC) for categorical traits. We found a marginal improvement in both NMRRE and PFC when including the phylogenetic vectors (phylogenetic NRMSE = 0.59, non-phylogenetic NRMSE = 0.63; phylogenetic PFC = 0.08, non-phylogenetic PFC = 0.09), suggesting phylogenetic signal in the traits. We then asked how congruent the imputed and non-imputed functional dendrograms were, as well as how congruent each was to the phylogeny, by estimating Generalized Robinson-Foulds distances (Smith 2021) using the SharedPhylogeneticInfo function in the TreeDist package (Smith 2020). We repeated all downstream analyses on both imputed and non-imputed datasets and found them to be qualitatively similar, driven by a very high correlation in functional distinctiveness metrics (see below) between the two approaches (r = 0.93; Supporting Information), and a high degree of shared phylogenetic information (77.9%). As such, we focus on the imputed dataset in the main text and present the non-imputed results in the Supporting Information.

Geographic, Functional, and Evolutionary Distinctiveness of Species

To understand the contribution of unusually distinct species to macroecological patterns, we estimated evolutionary distinctiveness (ED) and functional distinctiveness (FuD) for each Grevilleoideae species, and we used the negative log-transformed geographic range size as a measure of each species’ geographic range restriction (RR), such that high values correspond with greater range restriction (smaller range sizes). ED was estimated using the equal splits measure (Redding and Mooers 2006) in the picante R package (Kembel et al. 2010) on the species tree while retaining all sampled species. FuD was estimated following Skeels and Yaxley (2023), by first estimating a functional trait distance matrix using Gower distances with the daisy function in the cluster package (Maechler 2018), and then performing hierarchical clustering on this distance matrix using six different algorithms with the hclust function from the stats package (R Core Team 2021): UPGMA, WPGMA, WPGMC, UPGMC, complete linkage, and Ward’s minimum variance.  Hierarchical clusters were summarised as dendrograms representing a hierarchy of functional distinctiveness among species. WPGMC and UPGMC algorithms had negative branch lengths and were discarded. The remaining four dendrograms were quality assessed following the method of Maire et al. (2015), which estimates the mean squared deviation (mSD) of the original Gower distances and the distances estimated with the hierarchical clustering method. We found the UPGMA algorithm had the lowest mSD (0.02), followed by WPGMA (0.07), complete linkage (0.1), and Ward’s minimum variance (0.12). Using the UPGMA dendrogram, we estimated FuD similarly to ED, by using the equal splits measure (Redding and Mooers 2006) with the evol.distinct function in the picante package (Kembel et al. 2010). We quantified phylogenetic signal in each log-transformed metric using Pagel's λ and Blomberg's K with the phylosig function in the phytools package (Revell 2012). We estimated the linear relationships between log-transformed FuD, ED, and RR using phylogenetic generalised least squares models (PGLS) with the gls function in the nlme package (Pinheiro et al. 2020). We fitted PGLS for three different sets of relationships (FuD ~ ED,  RR ~ FuD; RR ~ ED) across 640 species with complete data, and for each compared three models for the phylogenetic correlation structure of the residuals (Brownian motion, Pagel's λ, and Ornstein–Uhlenbeck) and selected the best fit with model selection on AIC. Finally, to formally test whether the biota of each biome differs from the rest of the radiation, we fitted a phylogenetic ANOVA using the phylANOVA function in phytools (Revell 2012) for each biome and metric combination,~ w~here we contrasted biome member species against all non-members, and P-values were Benjamini–Hochberg-adjusted. Because many species occur in more than one biome, a species was treated as a member of every biome in which it occupies at least one ecoregion (Olson et al. 2001).

Geographic, Functional, and Evolutionary Distinctiveness of Assemblages

For calculations of assemblage distinctiveness measures, we defined assemblages as the species whose distributions overlap each 50 km × 50 km equal-area grid cell (hereafter “sites”) across Australia. We group assemblages into their constituent biomes (Olson et al. 2001), with the seven key biomes (with a simplified naming convention) in Australia being: Tropical Forests, Tropical Savannas, Temperate Forests, Mediterranean, Temperate Grasslands, Montane, and Deserts.

Species ranges were obtained from a recent study (Skeels et al. 2025) in which Grevilleoideae species distributions were estimated using alpha hulls derived from cleaned observation records from the Atlas of Living Australia (ala.org). Next, for each site, we estimated assemblage-level patterns of distinctiveness by taking the mean of log-transformed species-level metric values of ED, FuD, and RR for all species present. To estimate whether distinctiveness metrics were higher or lower than expected by chance, for each site, we computed the standardised effect size (SES) by comparing the observed distribution of distinctiveness metrics to a distribution generated under a null model (SES = (observed − mean null) / SD null). We used the curveball algorithm (Strona et al. 2014) implemented in the vegan package (Oksanen et al. 2019) to generate 999 null assemblages for every site. The curveball algorithm holds both cell richness and species occupancy frequencies fixed while randomising species identities across sites. Sites were classified as lower than expected (SES < −1.96), not significant (−1.96 ≤ SES ≤ 1.96), or higher than expected (SES > 1.96) for each metric. We then examined the pairwise associations among the three SES metrics (FuD vs ED; RR vs FuD; RR vs ED) by fitting, for each contrast, an ordinary least-squares (OLS) regression. To test whether assemblage-level distinctiveness metrics were associated with another commonly used metric of relatedness and divergence within assemblages, we estimated mean pairwise phylogenetic distances (MPPD) and mean pairwise functional distances (MPFD), and their SES values using the PhyloMeasures package (Constantinos and Brody 2015), which computes them analytically under a uniform null model rather than by randomisation. We fitted an OLS regression of ED_SES and MPPD_SES and of FuD_SES and MPFD_SES.

PE and FE each capture two of the three distinctiveness dimensions, by quantifying the geographic rarity of species within an assemblage, weighted by evolutionary history and functional trait distinctiveness, respectively. For each 50 km × 50 km grid cell, we estimated PE following Rosauer et al. 2009) and FE following Skeels and Yaxley 2023) using the phylo_endemism function in the phyloregion package (Daru et al. 2020), with the phylogeny and functional dendrogram, respectively. To explore the distribution of sites with narrowly distributed species and long branches or short branches of the phylogeny or functional dendrogram, as well as regions with a mix of both short and long branches, we used the categorical analysis of neo- and paleo-endemism (CANAPE) approach (Mishler et al. 2014). For phylogenetic endemism (PE), we simulated the distribution of 999 randomised communities using the curveball algorithm (Strona et al. 2014) in the canapeR package (Nitta et al. 2023), from which we estimated PE and relative PE (RPE). RPE is the ratio of the scaled PE values on the actual phylogenetic tree and the scaled values of PE on a tree with equal branch lengths (comparison tree). Sites were first identified as centres of endemism where PE on the observed or comparison tree was significantly high (one-tailed, α = 0.05); among these, paleo-endemic sites had significantly high RPE and neo-endemic sites significantly low RPE (two-tailed, α = 0.05); sites significant in PE but not RPE were classed as mixed endemism, and sites significantly high in PE on both trees at α = 0.01 were classed as super-endemic. We repeated the procedure using the functional dendrogram to estimate hotspots of functional endemism (FE) and regions of significantly long functional dendrogram branches, where narrow-ranged species are more functionally distinct than expected against the null; we refer to these as divergent-endemism hotspots. On the other hand, regions of significantly short functional dendrogram branches, where narrow-ranged species are less functionally distinct than expected against the null, we refer to as redundant-endemism hotspots.