We tested whether variation in Elongate dendritic morphometric traits reflected shared evolutionary history by assessing phylogenetic signal using both multivariate and univariate approaches. We estimated Pagel’s λ (Pagel, 1999) for each individual trait, except for the CC traits, which are highly interdependent and therefore better interpreted from the multivariate results. Univariate tests were conducted with phytools (Revell, 2024) using species means and standard errors, and significance was assessed against a null model assuming no phylogenetic structure.

Because fitted λ values alone are not reliable indicators of the strength of phylogenetic signal, we used multivariate effect sizes (Z) to compare across trait sets (Collyer et al., 2022). For the multivariate analysis, we used a matrix-based generalization of Blomberg’s κ (Blomberg et al., 2003; Mitteroecker et al, 2025), implemented in the R package geomorph (Baken et al., 2021; Adams et al., 2025). Analyses were performed separately for three trait sets: (1) geometric and branching traits, (2) PH1–PH20 traits, and (3) CC1–CC20 traits.

To further evaluate how phytolith traits relate to taxonomy, we tested whether Elongate dendritic morphometric traits could classify taxa at three hierarchical levels: (1) genus (lineage; see further), (2) species, and (3) subspecies within Triticum ( Table 1 ), using Random Forest models (Breiman, 2001). Random Forest was chosen because it performs well with relatively small sample sizes and complex, correlated predictors (Qi, 2012). All analyses were conducted in R (R Core Team, 2024), with scripts archived on Zenodo.

The dataset was unevenly represented, with some groups contributing far more phytoliths than others (e.g., Triticum taxa alone accounted for 622 of the 1,850 phytoliths). In machine learning, such class imbalance is problematic because models tend to overfit larger groups while underrepresenting smaller ones, leading to biased predictions (Japkowicz and Stephen, 2002). To reduce this imbalance, we trimmed the majority groups to better match the smaller ones and removed groups with fewer than 20 phytoliths. For genus-level classification, we also merged some of the ‘potential confuser’ genera into two broader groups: those belonging to the subtribe Poeae (Alopecurus, Arrhenatherum, Helictochloa, Helictotrichon, Lolium, Parapholis) and those belonging to the supertribe Triticodae (Aegilops, Bromus, Dasypyrum, Elymus, Hordelymus, Thinopyrum), resulting in six classes: Avena, Hordeum, Secale, Triticum, Poeae, and Triticodae ( Supplementary Data Sheet 2 ). Because this six-class system combines some genera into broader phylogenetic groupings, we hereafter refer to it as lineage-level classification.

After these steps, 1,563 phytoliths remained at the lineage level, 808 at the species level, and 248 at the Triticum subspecies level. Each dataset was partitioned into training (60%), validation (20%), and test (20%) sets ( Supplementary Data File 3 ) using stratified sampling at the inflorescence sample level. Approximately 20 phytoliths per inflorescence sample were distributed across the three sets in an 60/20/20 ratio.

The trimmed phytoliths from the majority groups were retained as a second evaluation set, termed ‘test_other’ ( Supplementary Data Sheet 2 ). The taxa represented in this set were the same as those used for training and validation, but the phytoliths themselves came from different plant specimens. The ‘test_other’ set therefore provided a stricter assessment of model generalization to phytoliths from new specimens of the same taxa.

Since both PH and CC traits quantify branching architecture—albeit in different ways—only one of the two sets was included as predictors alongside geometric predictors. Preliminary tests indicated that CC traits improved classification performance relative to PH traits. Although PH is theoretically more robust (Cohen-Steiner et al., 2007), CC was more computationally efficient and better suited to large datasets, while still demonstrating strong predictive power (Hacquard and Lebovici, 2024; Li et al., 2018). Models were implemented in the R package caret (Kuhn, 2008) and evaluated using classification accuracy and Cohen’s κ (Cohen, 1960). Confusion matrices were generated and visualized as heatmaps (Wickham, 2016).

A schematic overview of the morphometric workflow is presented in Figure 4 , summarizing the main steps of the study, from sample selection and chemical processing through confocal imaging, 3D surface mesh generation, extraction of morphometric and topological traits, phylogenetic and regression analyses, and finally machine learning classification.

Ancestral state reconstruction indicated that low occurrence of Elongate dendritic was likely ancestral within Pooideae (>85% posterior probability; Figure 5 ). Independent losses occurred in a few clades, such as Meliceae, while high production evolved multiple times independently, particularly within the Triticeae and Aveninae. The common ancestor of Triticeae likely had high Elongate dendritic occurrence (posterior probability = 0.79), as did the broader Triticodae crown group (posterior probability = 0.64). Similarly, Avena and several close relatives within Aveninae showed ancestrally high occurrence (posterior probability = 0.94). These results suggest that Elongate dendritic occurrence predates domestication, since it was already present in common ancestors, including the Triticodae crown group, which originated ~25 million years ago (Orton et al., 2021).

Out of 51 Elongate dendritic morphometric traits, 29 traits showed significant phylogenetic signal (λ = 0.98–0.46; mean = 0.75), indicating that closely related taxa tended to produce phytoliths with similar 3D morphometric traits. The strongest signals were observed in average dendritic length (λ = 0.98), solidity of the ‘core body’ (λ = 0.88), and maximum core length (λ = 0.86), while only four persistent homology (PH) traits (PH1, PH2, PH3, PH8) showed significant signal (λ = 0.85–0.68) ( Supplementary Data Sheet 3 ).

The multivariate analyses showed that the geometric and branching traits carried the strongest phylogenetic signal overall, followed by the PH traits, while the connect components (CC) traits carried the weakest signal ( Figure 6a ). Both the geometric and branching set and the PH traits produced significant results for two summary statistics quantifying the degree of phylogenetic signal (traceK and detK), whereas the CC traits were significant only for detK. Taxa also separated clearly along the first two K-components, which represent trait combinations capturing the strongest phylogenetic structure ( Figure 6b ; Mitteroecker et al., 2025).

Within the geometric and branching trait set, surface area, volume, and solidity contributed most strongly to K1, while average dendritic length, number of dendritic branches, and surface area were most influential along K2. For the PH traits, PH19, PH9, and PH10 loaded highest on K1, whereas PH18, PH20, and PH12 dominated K2. The traits most strongly associated with the K-components did not always align with those highlighted by the univariate λ results, especially within the PH trait set. This likely reflects differences in how λ and κ capture phylogenetic signal, as well as the influence of trait variance in univariate analyses (Collyer et al., 2022). Despite these differences, both univariate and multivariate results consistently indicated that morphometric traits captured meaningful evolutionary relationships.

These results demonstrate that Elongate dendritic carry a clear phylogenetic imprint, both in their occurrence patterns and in their 3D morphometric features, with closely related taxa tending to resemble each other more strongly.

Domestication was associated with differences in Elongate dendritic morphometric traits. Binary logistic regression pointed to a link between domestication status (wild vs. domesticated) and morphometric traits. Seven rotated components (RCs), selected from a scree plot, were used as input factors ( Supplementary Data File 3 ). The logistic regression model fit the data significantly better than the null model (χ²₇ = 22.93, p = 0.0018). Four RCs predicted domestication status (p < 0.05; Table 2 ): RC4 increased the odds (OR = 3.37, 95% CI [1.18, 13.20]), while RC6 (OR = 0.58, 95% CI [0.34, 0.93]), RC2 (OR = 0.53, 95% CI [0.30, 0.87]), and RC7 (OR = 0.54, 95% CI [0.31, 0.89]) were associated with reduced odds. The remaining RCs (RC1, RC3, RC5) were not significant ( Table 2 ).

This model suggested that the morphometric variation captured by specific rotated components (RC4, RC2, RC6, and RC7) may help discriminate between domesticated and wild species. RC2, highly loaded (loading > 0.60) by solidity of the ‘core body’, solidity, sphericity of the ‘core body’, and occupancy of the ‘core body’, indicated that domesticated taxa had a denser overall structure. RC7, highly loaded by flatness of the ‘core body’, indicated that wild taxa had a flatter inner structure. RC6, highly loaded by maximum height, maximum width, and maximum width of the ‘core body’, indicated that domesticated taxa generally had larger phytoliths. RC4 values, highly loaded by several CC traits, were higher in wild taxa (especially Aegilops), reflecting more complex dendritic branching in some wild taxa ( Figure 7 ; Supplementary Table S3 ). However, some domesticated taxa, particularly within Triticum, also exhibited high CC values, suggesting that certain cultivated taxa retained extensive branching traits.

Our LDA for domestication showed that LD1 provided some separation between wild and domesticated taxa. However, substantial overlap between individual phytoliths remained ( Figure 8 ), suggesting that although domestication influenced morphometric traits, high within-group variation overlapping with between-group variation limited the ability to reliably distinguish individual phytoliths. The model achieved a classification accuracy of 62.3%, which was statistically higher than the no-information rate of 54% (binomial test, p < 0.001), but the improvement was modest.

To investigate the relationship between morphometric traits and ploidy level (diploid < tetraploid < hexaploid), we applied ordinal logistic regression. Seven RCs, selected from the scree plot, were used as predictors in a cumulative logit model ( Supplementary Table S3 ). The ordinal logistic regression model did not significantly improve over the null model (χ²(7) = 5.30, p = 0.624). None of the RCs demonstrated a statistically significant effect, as their 95% confidence intervals all included 1 (the null hypothesis value) ( Table 3 ). For example, RC1 had an odds ratio of 1.27 (95% CI [0.84, 1.95]), indicating that increases in RC1 did not significantly alter the odds of transitioning from one ploidy category to a higher one. Although the threshold between tetraploid and hexaploid was substantially greater than zero (log odds = 1.167, SE = 0.260, t = 4.493), this shift was not explained by any of the retained RCs.

The LDA for ploidy level revealed minimal separation among individual phytoliths from diploid, tetraploid, and hexaploid taxa, with substantial overlap between groups ( Figure 9 ). A binomial test indicated that the observed accuracy (48.6%) was not significantly greater (p = 0.243) than the no-information rate (NIR = 47.7%), suggesting that ploidy level did not strongly influence the morphometric traits of individual Elongate dendritic.

All Random Forest models performed significantly better than chance (p < 0.05 across all tests; Table 4 ). The best classification performance was achieved at the lineage level (six classes; 844 entries), where the model reached a test accuracy of 49.5%, a strong improvement over the no-information rate (21.7%). At the species level (39 classes; 808 entries), classification accuracy decreased to 34.0% but remained far above the no-information rate (6.4%). At the Triticum subspecies level, the test accuracy was 31.9%, compared to a no-information rate of 14.9%.

Following the evaluation framework of Landis and Koch (1977), the test kappa values for lineage (0.38), species (0.30), and Triticum subspecies (0.25) indicated a fair level of agreement (κ ranging from 0.21 to 0.40) between predicted and true classes. Although classification accuracy declined with increasing taxonomic resolution, morphometric traits retained taxonomic signal across all levels. The reduced accuracy likely reflects either limited taxonomic resolution of phytoliths at finer scales or insufficient training data at lower taxonomic levels, where many classes had few samples.

The confusion matrices showed that, at the lineage level, Avena was sometimes misclassified as Secale, followed by Triticodae confusers, Hordeum, and Triticum. Poeae confusers, by contrast, were well separated. Hordeum was primarily misclassified as Avena or Triticodae confusers. Secale was commonly misclassified as Triticum or Triticodae confusers. Triticum was misclassified mainly as Triticodae confusers, but also as Secale and, to a lesser extent, as Hordeum or Avena. Finally, Triticodae confusers were most frequently misclassified as Triticum, but also as Hordeum and Secale ( Supplementary Table 4 ). At the species and Triticum subspecies levels, misclassifications were more scattered, likely reflecting greater trait overlap across taxa and the increased classification complexity due to the higher number of classes ( Supplementary Table S4 ).

To identify which traits the models relied on most, we examined the top five importance scores for the classification of lineages, species, and Triticum subspecies ( Figure 10 ). At the lineage level, the most important traits were ‘core body’ occupancy (relative to the entire phytolith), flatness, average dendritic branch length, and the solidity of both the whole phytolith and the ‘core body’. At the species level, the top-ranked traits included occupancy of the ‘core body’, dendritic branch density, and the standard deviations along the second and third principal components—related to phytolith width and height variation, respectively. For Triticum subspecies, the most influential traits were solidity, sphericity, total dendritic length, and the standard deviations along the second (~width) and third (~height) principal components of the ‘core body’.

To evaluate model generalizability, we tested classification accuracy on phytoliths from completely unseen inflorescence samples—i.e., those held out during the initial split due to being majority groups. While these particular inflorescence samples were not seen during training, they still belonged to taxa represented in the training set. At the lineage level, 287 Triticum phytoliths were included in the test_other set. The model correctly classified 53% of them as Triticum, while the remainder were misclassified as Triticodae confusers and, to a lesser extent, as Secale, Avena, or Hordeum ( Supplementary Table S4 ).

At the species level, we evaluated 991 phytoliths from 18 unseen species. Accuracy on this test_other set was 9.7%, considerably lower than the 33.9% accuracy observed on the standard test set. The observed accuracy did not exceed the no-information rate of 21.0% (p = 1.00), indicating that the model failed to generalize to these unseen taxa. This was likely the result of limited training data, where only a single inflorescence sample per species was included. Our results suggest that comparing phytoliths from an unseen specimen of a given species to just one specimen in the training set is insufficient for robust prediction. For Triticum subspecies, the test_other set included phytoliths from eight previously unseen samples. Accuracy was 15.8%, which did not differ significantly from the no-information rate (14.7%; p = 0.36), indicating that the model performed only marginally better than chance on these unseen individuals. As with the species-level model, the limited training data (only one inflorescence sample per subspecies) likely restricted model performance ( Supplementary Table S4 ).

In this study, we evaluated the distribution of Elongate dendritic across Pooideae inflorescence samples. The occurrence categories applied here (absent, low, high) should be regarded as estimates, with the primary aim of identifying where Elongate dendritic occurred across the Pooideae to ensure representative sampling for morphometric analysis. Nonetheless, these estimates revealed that high occurrence of Elongate dendritic appears to have evolved independently in multiple cereal-rich clades, including Triticeae and Avena. In contrast, complete absence was rare, observed only in a few lineages such as Meliceae, while many lineages exhibited low occurrence. This evolutionary pattern illustrates that the capacity to produce Elongate dendritic was phylogenetically widespread across Pooideae, but high occurrence was concentrated in select lineages—many of which include important domesticated cereals such as wheat, barley, rye, and oats, as well as non-cereal taxa ( Figure 5 ). These findings align with Novello and Barboni (2015), who reported that Elongate dendritic occurred frequently in wild African grasses. In archaeological contexts, Elongate dendritic should therefore not be automatically interpreted as evidence for cereals. This perspective is consistent with ICPN 2.0 guidelines, which recommend cereal attribution only in secure contexts, such as storage or cooking vessels, or specific thin sections (ICPT, 2019).

Albert et al. (2008) reported that an Elongate dendritic relative abundance of around 7–8% within a phytolith assemblage may indicate domesticated cereals. In future work, we aim to revisit the inflorescence samples and quantify relative abundances of Elongate dendritic more precisely. This will allow us to test whether relative abundances are consistently linked to domestication status or if they reflect deeper phylogenetic structuring.

Our results supported the first part of hypothesis H1, that domestication influences the morphometric traits of Elongate dendritic. Domesticated taxa generally produced larger, denser, and flatter phytoliths, while some wild taxa exhibited more complex dendritic branching ( Figure 7 ). Nonetheless, LDA revealed substantial morphometric overlap at the individual level, highlighting high intra-sample variation and the limited reliability of single phytoliths for distinguishing wild from domesticated forms ( Figure 8 ). These results are consistent with earlier studies linking domestication to increased phytolith size (Piperno, 1988; Piperno and Stothert, 2003; Pearsall, 1989; Vrydaghs et al., 2009), although more recent work on rice has shown considerable size overlap between wild and domesticated populations (Wang et al., 2019).

Our results did not support the second part of hypothesis H1, that ploidy level influences the morphometric traits of Elongate dendritic. No consistent morphometric patterns were detected across ploidy levels, suggesting that polyploidy had little effect on Elongate dendritic traits ( Table 3 ; Figure 9 ). As ploidy information was obtained from the literature, future work should include cytogenetic verification to refine these findings (Diet et al., 2022).

Our results partially supported hypothesis H2, that 3D morphometric traits of Elongate dendritic can distinguish taxa at different taxonomic levels within subfamily Pooideae. We detected significant phylogenetic signal in several traits, including solidity, core body solidity, maximum core body length, volume, surface area, average dendritic length, and number of dendritic branches ( Figure 6 ; Supplementary Data Sheet 2 ). Many of these traits also emerged as important for classification across taxonomic groups ( Figure 10 ). Among the classification approaches tested, Random Forest models performed best at the lineage level (accuracy 49.5%, κ = 0.38), with performance decreasing at the species (accuracy 34.0%, κ = 0.32) and Triticum subspecies levels (accuracy 31.9%, κ = 0.25) ( Table 4 ). These results support earlier findings that phytolith morphometry contains phylogenetically informative variation (Brightly et al., 2024; Chen et al., 2020; Gallaher et al., 2020; Hošková et al., 2020, 2021; Vrydaghs et al., 2009). The strong phylogenetic structuring of Elongate dendritic morphometric traits may also explain why it was difficult to separate wild from domesticated taxa, given their close taxonomic relatedness.

The results of this study are promising, indicating that classification based on isolated Elongate dendritic is feasible, particularly at the lineage level. However, the accuracy values achieved by our models were not necessarily high enough to be reliable for archaeological inference ( Table 4 ).

A likely explanation for the moderate classification accuracy, particularly at the species and subspecies levels, is the limited number of phytoliths analyzed per sample (approximately 20 per inflorescence). Such small sample sizes may not adequately capture the considerable intra-specimen morphometric variation. Increasing the number of phytoliths per specimen would likely mitigate this limitation. However, the use of 3D confocal microscopy remains both time-intensive and resource-demanding, which constrains the scalability of this approach. As a practical alternative, a 2D framework based on standard light microscopy could be developed, incorporating either manual or automated outlining of individual Elongate dendritic (see, e.g., Lloyd et al., 2024), followed by geometric and topological data analysis. Although this approach sacrifices three-dimensional detail, it would enable the acquisition of substantially larger and more representative datasets. Moreover, it would differ from earlier 2D analyses by Ball et al. (2017), which examined dendritic lobes in articulated phytoliths, by focusing instead on isolated Elongate dendritic, which broadens the number and type of samples that can be analyzed.

Future work could explore deep learning approaches which learn directly from 3D meshes such as MeshCNN (Hanocka et al., 2019), rather than relying on predefined features. Deep learning methods have shown promise for Elongate dendritic in 2D, both for morphotype recognition of isolated phytoliths (Andriopoulou et al., 2023) and taxonomic attribution of silica skeletons containing Elongate dendritic (Berganzo-Besga et al., 2022). However, applying deep learning approaches to our 3D dataset may also present challenges, including overfitting due to small sample sizes, preprocessing demands (e.g., mesh standardization), high computational requirements, and the black-box nature of deep learning algorithms in general.

Regarding model evaluation, the ‘test_other’ set—composed of phytoliths from the same taxa as those used for training but collected from different individual plants—provides a useful check on how well the model generalizes to new specimens of familiar taxa ( Supplementary Table S4 ). Currently, it represents only a subset of the taxa used for training and does not include any entirely new taxa outside the training set, so it cannot evaluate performance on species outside the training set. Expanding this set to include all taxa present in the training set, as well as taxa not used for training, would allow a more comprehensive assessment of model performance across both familiar and entirely new taxa. It is also important to evaluate how the model perform on Elongate dendritic originating from different biological, environmental or depositional contexts. These may influence morphometric traits in ways that affect classification performance. Several of these sources of variation are discussed below.

A range of factors may influence morphometric variation in Elongate dendritic, many of which we could not control in this study. For example, previous work has shown that dendritic morphometric traits vary within the inflorescence across bract types (glumes, lemmas, paleae), across spikelet or panicle positions (upper, middle, lower), and even between accessions of the same species (Ball et al., 2017). Other factors may influence Elongate dendritic morphology but are largely unstudied.

First, environmental conditions are known to influence phytolith morphology, with size often varying in response to factors such as soil moisture, temperature, and light regime, while shape appears comparatively stable across these gradients (Ball and Brotherson, 1992; Dunn et al., 2015; Hošková et al., 2021; Liu et al., 2016). For Elongate dendritic, however, no studies have specifically addressed how environmental or biogeochemical factors affect their morphometric traits. In addition, factors such as silica availability, nutrient balance, and lignification may be important, but remain to be investigated for Elongate dendritic.

Second, plant developmental stage may contribute to variation. For example, bulliform phytoliths in rice vary in size and ornamentation across growth stages (He et al., 2024), although Hošková et al. (2020) found minimal ontogenetic variation in Bilobate and Polylobate forms. Controlled growth experiments across developmental stages and inflorescence bract types (glumes, lemmas, paleae) will be essential to clarify how such variation might affect taxonomically informative traits.

Third, post-depositional and laboratory-induced alterations warrant consideration. Natural and anthropogenic post-depositional processes, including physical mechanisms such as breakage and erosion, chemical effects such as dissolution, and exposure to heat, could alter the morphology of Elongate dendritic (Cabanes et al., 2011; Cabanes and Shahack-Gross, 2015; Madella and Lancelotti, 2012; Wang and Shang, 2023). These alterations may affect features such as the dendritic branches; the psilate or nodulate surface texture; and the overall length of the phytolith, which may be shorter due to breakage. Laboratory preparation methods may also introduce morphological changes. Some samples were processed using wet oxidation alone, while others underwent a combination of wet and dry oxidation. Dry oxidation may cause deformation (shrinkage) in some phytoliths (Pearsall, 1989), likely due to dehydration during heating, although Parr et al. (2001) found no difference in Bilobate size between dry and wet oxidation. In addition, in modern inflorescence samples, isolated Elongate dendritic exhibit broken or missing dendritic processes or can be fractured in length, meaning their morphology is not always perfectly intact, even in reference material. Understanding how these processes affect the morphometric traits of Elongate dendritic, and whether such changes influence taxonomic attribution by classification models, is crucial.

Finally, broader taxonomic comparisons could help contextualize observed variation. While Elongate dendritic are most commonly produced in the inflorescence bracts of Pooideae, they have also been documented in other Poaceae subfamilies and, more sporadically, in unrelated monocot families such as Marantaceae and Arecaceae (Novello and Barboni, 2015; Ge et al., 2020; ICPT, 2019). Light microscopy could provide a first assessment of whether these forms are morphologically comparable to those in Pooideae. Regardless, biogeographic and ecological context remains essential: for example, Elongate dendritic found in tropical families such as Marantaceae are unlikely to be relevant in temperate archaeological contexts, where such plants would not have been present. Likewise, Elongate dendritic found in phytolith assemblages lacking Spheroid echinate or Spheroid ornate but containing grass silica short cells, are most likely derived from grasses.

As the Elongate dendritic dataset continues to grow and morphometric variation becomes better characterized, the next step should be to apply classification models to archaeological assemblages. Crucially, this does not require the inclusion of all species analyzed in this study. Our broad sampling strategy, designed to test generalizability across Pooideae, inevitably introduced high morphometric variation and thereby lowered classification accuracy. In contrast, archaeological contexts typically involve a narrower set of expected taxa, defined by region, chronology, and other contextual evidence. This restricted scope opens the possibility of building focused, high-accuracy models tailored to specific case studies. For example, in a medieval Northwest European context, there would be no need to include non-local progenitor taxa such as Triticum urartu or T. turgidum subsp. carthlicum, and excluding non-applicable taxa could significantly improve model performance.

More broadly, we hope that future archaeological research can integrate cereal taxonomic questions more explicitly into phytolith analysis. An ideal future outcome from this project would be the development of a tool with a user-friendly interface for attributing isolated Elongate dendritic, one that reduces reliance on extensive human training or subjective visual comparison. By leveraging models trained on morphometric patterns, such tools could standardize identifications, reduce inter-observer bias, and make phytolith analysis more accessible to a wider range of researchers. In turn, this would support more consistent and reproducible interpretations in archaeological contexts.