Frond orientations with independent current indicators demonstrate the reclining rheotropic mode of life of several Ediacaran rangeomorph taxa

The Ediacaran biotas of Newfoundland have some of the oldest known well-dated complex soft-bodied macroorganisms (Narbonne 2005; Xiao and Laflamme 2009; Liu et al. 2015; Matthews et al. 2021) and are essential in furthering our understanding of the evolution of complex macroorganisms. The earliest macrofossils of the Ediacaran biota are widely considered to be stem eumetazoans (Dunn et al. 2021), and many had unusual fractal-like growth patterns that have no good modern analogues (Narbonne 2005; Brasier et al. 2012; Hoyal Cuthill and Conway Morris 2014). Understanding these enigmatic organisms has been a challenge since their initial discovery (Gürich 1930; Ford 1958). The Ediacaran biota has been much explored in recent years, with research debating aspects of their phylogenetic affinity (Seilacher 1989, 1992; Budd and Jensen 2017; Dunn and Donoghue 2018), ecology (Clapham et al. 2003; Darroch et al. 2018; Mitchell et al. 2020), mode of life (Glaessner 1985; Laflamme et al. 2009; Dufour and McIlroy 2017; McIlroy et al. 2021), and taphonomy (Gehling 1999; Narbonne 2005; Liu et al. 2011; Bobrovskiy et al. 2019). Inferences that most of the Ediacaran biota—particularly the fractal-branching frondose Rangeomorpha—lived erect within the water column during life (i.e., Glaessner 1985; Laflamme et al. 2007) have recently been contested (McIlroy et al. 2020, 2021, 2022; Pérez-Pinedo et al. 2022). The high fidelity of preservation of originally soft tissues such as fractal-like fronds (Fig. 1) and the lack of taphonomic evidence for an erect mode of life (e.g., swing marks imparted by a current, see Jensen et al. [2018]) challenge a priori assumptions about erect body postures. There is therefore a need to consider the alternative hypothesis that these organisms were reclined on the seafloor, as they are preserved on bedding planes (McIlroy et al. 2020, 2021, 2022; Pérez-Pinedo et al. 2022).

Sedimentological evidence for paleocurrent directions in the Ediacaran of Newfoundland is mainly in the form of current ripple trends in turbidites (i.e., Wood et al. 2003; Ichaso et al. 2007; Mason et al. 2013). At some localities, geostrophic paleocurrent directions are also inferred but are solely based on frond orientation. Preferential orientation of Ediacaran fronds has been thought to be due to current “felling” in the Mistaken Point–type biotas of Avalonia, particularly at Mistaken Point (following Seilacher 1992, 1999). The poorly supported felling model has been used to propose an erect mode of life for all the fossil organisms in the Mistaken Point assemblages, except Fractofusus spp. (e.g., Wood et al. 2003; Vixseboxse et al. 2021). An alternative model invoking rheotropic epifaunal growth of reclining taxa has been proposed (McIlroy et al. 2022), which calls into question the causal relationship between frond orientation and paleocurrents. There is therefore a clear need to assess the orientation of Ediacaran fronds from surfaces with physical sedimentary structures as independent proxies for paleocurrent direction. The clearest documented relationships between frond orientation and current structures are the fronds of the “Spaniard's Bay” assemblage (King 1982, 1988, 1990; Narbonne 2004; Narbonne et al. 2009) that are oriented subparallel to flute marks (cf. Brasier et al. 2013), but such clear relationships are uncommon in the rest of the Ediacaran assemblages of Avalonia. This paper details the orientation of Ediacaran taxa from a newly discovered fossil surface in the Fermeuse Formation at Melrose, on the southern portion of the Catalina Dome in the Discovery UNESCO Global Geopark (Fig. 2). The importance of this surface derives from its well-constrained paleocurrent indicators in the form of wavy bedding and climbing current ripples (Fig. 3).

The fossil assemblage at Melrose bears some resemblance in terms of taxonomy and community composition to that of the famous “D” and “E” Surface assemblages of the Mistaken Point Ecological Reserve (MPER) (see Clapham et al. 2003). The assemblage was discovered in 2018 by a team from Cambridge University (Emily Mitchell, Alexander Liu, and William McMahon), U.K., and is here named the “Melrose Surface,” given its location in the town of Melrose. Tectonostratigraphic reconstructions of the Conception Group suggest that deposition occurred in a continental slope environment adjacent to an island-arc setting on the northwest margin of Gondwana (Murphy et al. 2004; Pollock et al. 2009; Pisarevsky et al. 2012). The Conception Group is dominated by fine-grained turbidites (Williams and King 1979; Wood et al. 2003; Ichaso et al. 2007). The tuffite overlying the “E” Surface at MPER is dated at 565 ± 0.64 Ma (Matthews et al. 2021), slightly older than the undated Fermeuse Formation at Melrose (Fig. 2).

Ecological reconstructions of Mistaken Point–type biotas have hitherto depicted these Ediacaran organisms as living erect within the water column, influenced by and felled in the direction of paleocurrents. However, due to the lack of sedimentological proxies on such fossiliferous beds, paleocurrent direction has been inferred in some cases from frond orientation. Paleocurrent analysis of the Mistaken Point Formation has shown a predominantly southeasterly flow direction, as evinced from both the abundant fronds and the uncommon physical sedimentary structures in turbidite beds (Wood et al. 2003; Ichaso et al. 2007). It has, however, also been noted that there is a second preferential frond orientation direction that is roughly orthogonal to the downslope current, previously attributed to contour currents that felled the fronds before rapid ash burial; this aided in the high fidelity of preservation seen on these fossiliferous surfaces (Seilacher 1992, 1999; Wood et al. 2003; Narbonne 2005). There is a range of frond orientations on the MPER surfaces, and some taxa are commonly orthogonal to slope (Beothukis mistakensis; Hawco et al. 2020; McIlroy et al. 2020, 2022); some even have fronds oriented against the inferred paleocurrent direction (e.g., Bradgatia sp., orientations in Flude and Narbonne [2008]). This is more likely to be a rheotropic growth response, as it is incompatible with flow hydrodynamics (McIlroy et al. 2022; contra Vixseboxse et al. 2021).

The “Melrose Surface” assemblage, while similar in taxonomic character to the older Mistaken Point assemblages, is part of the later Fermeuse Formation of the St. John's Group (Hofmann et al. 2008) (Fig. 2). The lithofacies of the Fermeuse Formation suggest deposition in an upper slope to distal prodelta setting (King et al. 1988; Mason et al. 2013), and thus in shallower waters than the ecologically similar Mistaken Point biota, which is from a deep basinal setting. The presence of what is normally considered a deep basinal biota in the shallower Fermeuse Formation requires detailed sedimentological consideration beyond the scope of this work, but it may be possible to invoke a regional increase in relative sea level. The fossiliferous surfaces of the Fermeuse Formation are typically dominated by matgrounds with Ediacaria and Cyclomedusa morphs of Aspidella sensu lato (Gehling et al. 2000; Brasier et al. 2010) and the protist-like Palaeopascichnus delicatus (Hawco et al. 2021). Beds deposited in similar paleoenvironments with discoidal features comparable to Aspidella have been used to suggest that these discoidal features may be pseudofossils, possibly formed by fluid escape or matground-associated sediment injection rather than by organisms (Menon et al. 2016). The Melrose assemblage considerably extends the upper stratigraphic ranges of several taxa (Fractofusus misrai, Pectinifrons abyssalis, Primocandelabrum sp., Bradgatia sp., and Charniodiscus sp.).

Traditional paleobiological reconstructions of stemmed frondose Ediacaran fossils have followed Glaessner (1985), interpreting these fossils as pennatulacean-like organisms that lived erect in the water column. The descriptive terminology used (i.e., holdfast, stem, frondlet) is unfortunate, in that it evokes affinities to extant taxa and infers biomechanical properties and life habit. The most significant divergence from this paradigm of ubiquitously erect Ediacaran macroorganisms was the work of Seilacher (1992, 1999), which suggested a more protistan-like mode of life for some taxa and inferred a reclining mode of life for others, such as Beothukis. This alternative hypothesis—in which life position should be the same as that in which fronds are preserved (McIlroy et al. 2020, 2021) unless there is evidence to the contrary—has been successfully applied to demonstrate that Charniodiscus concentricus (and possibly the frondose portion of Charniodiscus procerus), Arborea arborea, Arborea spinosa, and Arborea longa were likely held erect or recumbent in the water column (Laflamme et al. 2004, 2018; Pérez-Pinedo et al. 2022).

The “Melrose Surface” assemblage is taxonomically similar to the “D” and “E” Surface assemblages at Mistaken Point (shared occurrence of all taxa in the “Melrose Surface”) and has a comparable community composition but shows a much lower generic diversity (5 vs. 8 [“D”] and 12 [“E”]), a lower Shannon diversity (0.375 vs. 0.7 [“D”] and 1.52 [“E”]), and a lower Shannon evenness (0.23 vs. 0.33 [“D”] and 0.61 [“E”]) (cf. Clapham et al. 2003).

Fractofusus misrai, the most abundant species in the Mistaken Point Ecological Reserve (i.e., Mitchell et al. 2015; Mitchell and Butterfield 2018), dominates the “Melrose Surface” assemblage (91% of all specimens) but is otherwise very rare in the Catalina Dome (cf. Hofmann et al. 2008). The second most abundant species on the surface is P. abyssalis (6% of the total specimens), which is occasionally found on surfaces with F. misrai (i.e., MPER “D” Surface [Clapham and Narbonne 2002; Clapham et al. 2003; Bamforth et al. 2008], although we note the lack of statistically significant co-occurrence indices in Eden et al. [2022]) and is otherwise unreported from the Catalina Dome. Pectinifrons abyssalis has been interpreted as having frondlets elevated into the water column (Bamforth et al. 2008). Similarly, “multifoliate” rangeomorphs such as Bradgatia sp. and Primocandelabrum sp. are poorly represented in the assemblage (2% and 0.5%, respectively). These have typically been interpreted as living erect in the water column (e.g., Flude and Narbonne 2008; Hofmann et al. 2008) and preserved by felling by tuff-laden turbidity currents (though this is at odds with the hydrodynamics of the inferred turbidity currents, as noted by McIlroy et al. [2022]). A single Charniodiscus sp. is the only arboreomorph fossil on the “Melrose Surface.” Charniodiscus has also been occasionally found in association with F. misrai, P. abyssalis, and Bradgatia (cf. “D” and “E” Surfaces at MPER; Clapham and Narbonne 2002; Pérez-Pinedo et al. 2022; but see Eden et al. 2022).

Fractofusus misrai was discovered by Anderson and Misra (1968) in the Mistaken Point Ecological Reserve, where it is restricted to the Mistaken Point Formation and was named by Gehling and Narbonne (2007) along with Fractofusus andersoni. In the Catalina Dome, F. misrai is otherwise only known from three ambiguous specimens from the Murphy's Cove Member, which has been correlated with the Mistaken Point Formation of the Avalon Peninsula as well as a possible specimen in the Fermeuse Formation (Hofmann et al. 2008).

Orientation data of F. misrai from the “E” Surface at Mistaken Point (Seilacher 1992; Gehling and Narbonne 2007) suggest no preferential orientation of the population. The appropriateness of log transformations to study population structure is much debated, and two approaches have been conducted on morphometric traits. Clustering algorithms and Bayesian information criterion (BIC) analysis of F. misrai on transformed data from both the “D” and “E” Surfaces suggest that both assemblages were composed of specimens of different ages and with a large size range, suggesting a slow-growth model with aseasonal (continuous) reproduction (Darroch et al. 2013). Similarly, alternative studies on non-transformed data suggest that the F. andersoni assemblage from the “H14” Surface (Bonavista Peninsula) represents a multigenerational population, possibly the result of continuous asexual reproduction via stolons or propagules (Mitchell et al. 2015).

Pectinifrons abyssalis was described from the Mistaken Point Ecological Reserve (“MPD,” “MPN,” “SH,” and “Watern Cove” Surfaces) and Green Head (Bamforth et al. 2008) and is considered to be a rangeomorph with a curved pedicle rod and two series of parallel first-order rangeomorph branches that are typically found on the concave side of the rod. The only known exceptions to this branch positioning are two specimens from the “Mistaken Point North” Surface that may have been twisted or folded before burial (Bamforth et al. 2008). The first-order rangeomorph branches of Pectinifrons are interpreted as being held above the seafloor, emerging from the upper surface of the pedicle rod. No preferential orientation was found at the Mistaken Point Ecological Reserve, whereas an apparent bimodal orientation was reported from Green Head (Bamforth et al. 2008).

Bradgatia was originally described by Boynton and Ford (1995) based on material from the Charnwood Forest (U.K.), while the Newfoundland material was later reviewed by Flude and Narbonne (2008). Bradgatia linfordensis is a multifoliate rangeomorph, showing several first-order branches stemming from a central region. The fronds have typically been interpreted as being held in clusters subvertically in the water column and felled by depositional events. Taphomorphs of Bradgatia include elongate specimens (“I-shaped”) to specimens showing first-order branches oriented in all directions (“O-shaped”), while in others the fronds are all on one side of the point of origin (Flude and Narbonne 2008).

It has been previously noted that Bradgatia sp. are commonly bimodally oriented in both the inferred up-current and down-current directions on the “G” Surface at MPER (Flude and Narbonne 2008), while at other localities (“E” and “D” Surfaces at MPER) there are also abundant specimens oriented perpendicular to the paleocurrent (Flude and Narbonne 2008; Vixseboxse et al. 2021; McIlroy et al. 2022). Orientation data from Bishop's Cove potentially show unimodal orientation, with several specimens oriented orthogonally and against the paleocurrent direction (Flude and Narbonne 2008). These orientation data have not, however, been incorporated into paleobiological models of the genus.

The data presented in this study were collected from a single fossiliferous surface of the Fermeuse Formation that crops out on the southeastern margin of the Catalina Dome, Bonavista Peninsula (see O'Brien and King 2005) (Fig. 2). The surface is part of a succession of thin- to medium-bedded gray and green fine sandstones and siltstones, with some coarse-grained sandstones (Fig. 4). The “Melrose Surface” fossils are preserved on a 2-cm-thick gray to dark gray wavy bedded siltstone, which is overlain by a 5- to 10-cm-thick sandstone with soft-sediment deformation. Under the siltstone bed is a 0.1- to 0.6-mm-thick ripple cross-laminated sandstone bed that indicates a paleocurrent orientation of 102° SE (78° NE with respect to magnetic north); this feature is consistent with the cleavage-enhanced wavy bedding. The fossiliferous horizon is ferruginous, probably reflecting the weathering/oxidation of pyrite.

Data collection was grid based, with 41 squares of approximately 1.5 × 1.5 m completely documented with respect to their paleontology. Grids were photographed and orientation data collected (n_t = 208 specimens of five different taxa measured: Fractofusus n_f = 190; Pectinifrons n_p = 12; Bradgatia n_b = 4; Primocandelabrum n_pr = 1; and Charniodiscus n_c = 1). Orientation was recorded for four additional Bradgatia specimens from an adjacent slab where access is more difficult, which belongs to the same surface.

Fragmental specimens of all studied taxa were excluded from the analysis due to analytical difficulties and potential to introduce sampling bias. Only Fractofusus was included in the statistical analysis due to the low numbers of well-preserved specimens of all other taxa. Given the low numbers of kinked Fractofusus specimens (seven specimens) and given the bidirectional straight nature of the remaining specimens, orientation was projected under and over 180° separately, opposed to projections on 360° (following Gehling and Narbonne 2007; Mitchell et al. 2015; Vixseboxse et al. 2021). In this way, we avoided, plotting two images of virtually identical vectors, which could have resulted in misleading interpretation of the data. Climbing ripple cross-lamination was used to infer a southeasterly paleocurrent orientation of 102° (78° to magnetic north) (cf. Mason et al. 2013) (Fig. 3), providing an independent and unequivocal current direction against which to compare frond orientations.

In this study, two quantitative morphometric traits (length and width) of the taxa were recorded and corresponding size-frequency distributions were analyzed (i.e., Darroch et al. 2013). Morphometric traits were retrieved from field photography and analyzed with the software ImageJ (Schneider et al. 2012). Study of size-frequency distributions is a commonly used method in biology to explore population structure (Darroch et al. 2013). It can provide insights regarding the coefficient of variation across different morphometric variables, the skewness or relative contribution of different age/size classes, and the mode as a proxy for mortality rate and reproductive behavior (e.g., Bak and Meesters 1999; Meesters et al. 2001; Darroch et al. 2013). First, normality of the data was assessed by Shapiro-Wilk test, and logarithmic transformation of the data was conducted as required, as it normally generates a more precise representation of population structure (following Darroch et al. 2013) (Fig. 5). However, alternative approaches with non-transformed data have also been conducted, revealing multigenerational populations (e.g., Mitchell et al. 2015). The data were analyzed by the Gaussian finite mixture model-based clustering algorithms of the package mclust (Scrucca et al. 2016). This method resolves the most likely number of modes of the morphometric traits that represent different age/size classes within the population. The best model was chosen by a criterion for likelihood-based model selection (BIC) (see Darroch et al. [2013] for full methodology). Because univariate size-frequency distribution analyses could result in misleading conclusions, bivariate analyses were also conducted (Fig. 5).

Variables that are measured in scales that are circular or cyclical rather than linear or continuous commonly arise when recording orientation or time, both of which can be accurately represented as a circumference. Using linearized scales for circular variables can lead to misinterpretation of the natural behavior of the circular variables and therefore alternative statistical treatment is required (Landler et al. 2020). Herein we apply an integrated approach to analyze multivariate datasets including quantitative and circular variables together, without treating orientation as a continuous variable. This approach allows us to interpret information derived from all the recorded variables integrated.

The circular variable of species-specific orientation was initially visualized by plotting nonparametric density curves with angular histograms, where bandwidths were chosen by visual examination (Will 2016) (Fig. 6). To test potential departures from uniformity in the distribution of the circular variable, several statistical tests were considered (see Vixseboxse et al. 2021). A Rayleigh test, designed to detect unimodal von Mises (the circular analogue of the normal distribution) departure from uniformity, was conducted (Batschelet 1981); however, the Rayleigh test has very low potential to detect multimodal distributions. Rao's spacing test was conducted to test for deviations from uniformity without the assumption of a von Mises distribution. Both visualization of angular histograms and Rayleigh and Rao's spacing circular uniformity tests were performed using the package circular in R (Agostinelli and Lund 2017; Vixseboxse et al. 2021). However, the relative performance of these tests against multimodal departures from uniformity bears uncertainty (Landler et al. 2018). As multimodal deviations from uniformity are suspected based on the summary plots (Fig. 6, see Supplementary R script), and following Landler et al. (2018, 2019, 2020), the Hermans-Rasson test was chosen as the most suitable test to explore multimodal departures from normality. The test was performed using the package CircMLE (Fitak and Johnsen 2017; Landler et al. 2019; Vixseboxse et al. 2021), reducing type I error (Landler et al. 2020). The hypothesis of a von Mises distribution was also tested using Watson's goodness of fit (Agostinelli and Agostinelli 2018; Vixseboxse et al. 2021) (see all results in Table 1). As a preliminary exploration of orientation, Gaussian finite mixture model-based clustering algorithms were employed to identify the most likely number of modes in orientations from under and over 180° (Fig. 6). This analysis was deemed appropriate, as the maximum degree of dissimilarity of the data was 180° (in the case of Fractofusus specimens due to their bipolar growth).

Traditional parallel coordinate plots (PCPs) are sometimes employed to explore the relationships between quantitative variables, treating circular variables as regular continuous variables. These variables can be scaled, with corresponding values in each variable correlated by connecting lines. Although this approach allows the observation of the relationships in multivariate datasets (Hardle and Simar 2019), it is unsuitable for circular variables, because extreme values (0° and 360°) would bear the greatest degree of dissimilarity, when in fact they are essentially equal. Therefore, modified PCPs were used to visualize multivariate data with a circular variable. A circle of radius 0.5 was chosen to represent orientation drawn as an ellipse (sensu Will 2016). An aspect ratio of 1:1 would generate a circle but would affect the readability of the figure; the mechanics of the PCP do, however, remain unaltered. This approach allows the visualization of relationships across circular and continuous variables. Because the projection of an ellipse could visually alter the spatial interpretation of the figure, density curves used in univariate summary plots were also projected (Figs. 7, 8, Supplementary R script).

To further explore these patterns in the dataset, we conducted cluster analysis on a dissimilarity matrix. Gower's distance was used to account for dissimilarity that introduced a circular variable (Will 2016). In the case of circular variables, distance is measured in either direction around the circumference, with 180° being the largest possible difference. The function for implementing a dissimilarity matrix with a circular variable was retrieved from Will (2016). Gower's general coefficient has applications in the package ade4 (Dray et al. 2007).

After generating a distance matrix including the circular variable, we performed cluster analysis. Ward's hierarchical bottom-up agglomerative polythetic clustering method was conducted using the hclust function in R (Will 2016). All groups of clusters were processed, and the cluster that resulted in the smallest increase in error sums of squares was selected (Everitt and Hothorn 2011). We employed circular variables and a Gower's dissimilarity matrix; therefore, the distance is non-Euclidean. It accounts for a pseudo-error sum of squares (Will 2016). Clusters were visualized in a dendrogram in which height is a nondimensional indication of the dissimilarity between clusters (Figs. 7, 8). Potential solutions were indicated by crosscutting horizontal lines. When the vertical dissimilarity between splitting clusters is of a small magnitude, the real existence of subsequent clusters is ambiguous. To resolve the most parsimonious solution, the Caliński-Harabasz pseudo F-statistic (Caliński and Harabasz 1974; Will 2016) was employed (Figs. 7, 8), available in the clusterSim package (Walesiak and Dudek 2020). This test explores the relationship of between-clusters sums of squares compared with within-cluster sums of squares across several cluster solutions (k). However, pseudo sums of squares are considered in this case (Will 2016). The largest value of G1 (ratio of between-cluster variability over within-cluster variability) indicates the best solution for the number of clusters. The potential one-cluster solution has no mathematical value, because the calculation is not allowed (Will 2016). The cluster solution was visualized by modified PCPs wherein individual observations belonging to clusters were colored differently. To further facilitate the readability of the figure, medoids of the clusters were projected. A medoid of a cluster is “an observation that has the smallest average dissimilarity with all the other observations in the cluster” (Will 2016: p. 19) (Figs. 7, 8) following Kaufman and Rousseeuw (2009) (see Will [2016] and Tran [2019] for full methodology of circular statistics). The viability of the one-cluster solution was inspected by uniformity and von Mises distribution tests (Vixseboxse et al. 2021) conducted before the cluster analysis described earlier and resolved by modified monothetic clustering analysis (Tran 2019).

Monothetic cluster analysis was additionally conducted to support the chosen solution for number of clusters on the variables suspected to drive clustering patterns and test the viability of the one-cluster solution (see Tran 2019). This test is “a divisive clustering method that uses a hierarchical, recursive partitioning of multivariate responses based on binary decision rules that are built from individual response variables” (Tran 2019: p. xiv). Estimation of the optimal number of clusters solution was based on a combination of M-fold cross-validation (Hastie et al. 2016) and permutation-based hypothesis test at each split (Hothorn et al. 2006) adapted to work with monothetic cluster analysis and circular variables with non-Euclidean distances (see Tran 2019).

Because MSE can be calculated for the one-cluster solution, the cross-validation method can compare the one-cluster solution with multi-cluster solutions (Sneath and Sokal 1973; Chavent 1998; Piccarreta and Billari 2007). This test was combined with permutation-based hypotheses, which involved testing the null hypothesis of independence (H₀ = splitting two clusters are identical) at each node of the tree retaining statistically significant p-values (Hothorn et al. 2006; Tran 2019). This method was conducted using the monoClust package in R (Tran et al. 2021).

Untransformed length- and width-frequency distributions were right-skewed (Fig. 5) and moderately depart from normality (Shapiro-Wilk test, α = 0.01, p-value [length] = 0.0066, p-value [width] << 10⁻⁴). Log-transformed quantitative variables (following Darroch et al. 2013) were normally distributed (Shapiro-Wilk test, α = 0.01, p-value [length] = 0.0170, p-value [width] = 0.6448). Univariate analysis conducted using Gaussian finite mixture model-based clustering algorithms resulted in a single mode or age/size group as the most likely cluster solution assuming both equal and unequal variances (Fig. 5). Subsequent bivariate analysis generated best-fitting models assuming ellipsoidal, diagonal, and spherical distributions. Only the models assuming ellipsoidal distributions allowing for unequal variances are biologically realistic (Fraley and Raftery 2007; Darroch et al. 2013). BIC results from the bivariate models derived from the clustering algorithms are congruent with the univariate models in selecting one cohort as the optimal solution (Fig. 5; see Supplementary R script for BIC values,).

Regarding circular variables, the angular histogram with density curves suggested multimodal (likely bimodal) distributions of essentially identical recorded orientations both under and over 180°. These distributions were supported by the results from the Gaussian finite mixture model-based clustering algorithms, which resolved two modes as the most likely number of clusters in the orientation-frequency distribution (under 180°: 43° [NE] and 140° [SE]; over 180°: 222° [SW] and 321° [NW]) (Fig. 6). All the tests employed (Rayleigh, Rao's spacing, and Hermans-Rasson) rejected the null hypothesis of uniformity (p-value < 0.05) as well as the existence of von Mises distributions (Watson's test), indicating nonnormal and nonrandom species-specific orientations (see all results in Table 1).

Modified PCPs were congruent with angular histograms and the results from the Gaussian finite mixture model-based clustering algorithms (Fig. 6) in suggesting observable subgroups in recorded orientation: NE,SE and SW,NW under and over 180° respectively (Figs. 7, 8). On the other hand, quantitative variables showed a more unified trend with not easily recognizable subgroups (Fig. 5). The dendrogram cluster obtained by Ward's hierarchical bottom-up agglomerative clustering method on the dissimilarity matrix resulted in two groups, with a three-cluster solution bearing little magnitude between splitting clusters as indicated by cross-cutting horizontal green and red lines (Figs. 7, 8). Caliński-Harabasz pseudo F-statistic resolved the two-cluster solution as the most likely scenario, with the largest value of G1 for both orientations under and over 180° (Figs. 7, 8). Finally, the medoids of the clusters were projected (under 180°: 43° [NE] and 130° [SE]; over 180°: 223° [SW] and 322° [NW]) giving essentially identical orientation results to those obtained from the orientation-frequency distribution analysis (Tables 2, 3).

Monothetic clustering can operate on non-Euclidean distances, and it was employed to compare the selected two-cluster solution to a potential one-cluster solution. M-fold cross-validation (CV) showed a great difference between the one-cluster and the multi-cluster solutions. The solution made by 10-fold CV shows multi-cluster outcomes showing a great decrease in CV_k (CV_k; mean squared error [MSE]) and reduced error bars (MSE ± 1 SE) for both under and over 180° (Fig. 9). The cluster solution derived from the permutation-based hypothesis test confirmed two clusters for under and over 180° driven by orientation (Fig. 9). The two confirmed clusters driven by orientation, and illustrated by the medoids, reveal two groups of Fractofusus, with one group being oriented NE-SW and the other oriented SE-NW.

Only 11 specimens of P. abyssalis are reported from the “Melrose Surface” and were therefore excluded from the clustering analyses. The specimens do not show a clear orientation pattern, ranging from being almost parallel (n_p = 2) to the inferred paleocurrent direction to being almost perpendicular (n_p = 7), and even oriented against it (n_p = 2) (Table 4, Fig. 10).

One specimen appears to be twisted (Fig. 1B), showing first-order branches on the right-hand side of the pedicle rod in the south half of the fossil and on the left-hand side on the north half. Preservation of the first-order branches appears to be similar between the two halves, but no higher-order branches can be observed.

Only four specimens of Bradgatia are reported from the gridded section of the “Melrose Surface” and four additional specimens are located on an adjacent slab. Due to the low number of specimens, they were therefore excluded from the clustering analyses. The four specimens from the gridded section have the “U-shaped” morphology (sensu Flude and Narbonne 2008) and are characterized by first-order branches arising from a central position, and therefore determine an apical–basal polarity. Three specimens are oriented in a close to orthogonal position to the paleocurrent direction, and another specimen is oriented with the frond tips in the up-current direction. The remaining four specimens show parallel and oblique orientations to the inferred paleocurrent (Table 4, Fig. 10).

The Ediacaran-aged Avalon assemblage (Waggoner 2003) is characterized by immotile epibenthic communities of macroorganisms that are preserved in situ. In the absence of evidence for motile scavengers, and given the rarity of motile macroorganisms in general (Brasier et al. 2010; Liu et al. 2010; although see horizontal burrows and surface trace fossils in Liu et al. [2014] and Liu and McIlroy [2015]), ecological mechanisms such as competition and dispersal are thought to have led to gradational epifaunal tiering (vertical subdivisions of the water column) and evolutionary innovations in Ediacaran communities (Clapham and Narbonne 2002). A tiered epifaunal structure is congruent with suspension feeding or organic matter absorption as a means of resource partitioning (Bottjer and Ausich 1986), whereby nutrients are provided by bottom currents. The evolution of stems in demonstrably erect taxa increased height and allowed for exploitation of higher tiers, resulting in a competitive advantage and increased feeding capacity (Clapham and Narbonne 2002; but see Mitchell and Kenchington 2018). Reclining organisms such as the spindle-shaped Fractofusus and Beothukis, which lay flat on or partly within matgrounds, occupied the basal tiers.

Previous inferences on Ediacaran feeding strategies such as suspension feeding or osmotrophy in the Rangeomorpha are based upon analogies with extant Pennatulacea and are inferred based on comparable epifaunal tiering (Laflamme et al. 2009). However, the morphological analogy has been demonstrated to be imperfect (Antcliffe and Brasier 2007), particularly given the absence of preserved pores for fluid flow/organic matter uptake (Liu et al. 2015; but see Butterfield 2022). The organic matter reservoirs of dissolved organic matter (DOM) (Johnston et al. 2012; Liu et al. 2015) are unlikely to have been sufficiently concentrated to support osmotrophy (Fakhraee et al. 2021; McIlroy et al. 2021). In the absence of good evidence for either suspension feeding or osmotrophy, it has been proposed that some Ediacaran organisms could have relied on chemosymbiotic strategies (i.e., Seilacher 1984; McMenamin 1998; Dufour and McIlroy 2017; McIlroy et al. 2021, 2022; Taylor et al. 2021; Pérez-Pinedo et al. 2022).

As Ediacaran seafloors lacked significant amounts of macrobioturbation before ~550 Ma (Jensen 2003), the abundant matground-associated seafloors are likely to have experienced buildup of toxic sulfides in their porewaters (Dufour and McIlroy 2017). This toxicity would have posed a threat to sessile epifaunal sediment–reclining organisms like Fractofusus and Beothukis, which maintained constant association with the Ediacaran sulfide-rich seafloor. The basal surface of Fractofusus is characterized by a large surface area, which has been considered an adaptation for maximizing a chemosymbiotic epithelium, and to allow transport of oxygen into porewaters to mitigate the impact of reduced forms of sulfur on the lower surface. These organisms are thought to have increased oxygen content at their lower surfaces by diffusion through the body and/or ciliary irrigation. Increasing oxygen supply to sulfidic porewaters increases productivity of chemolithoautotrophic sulfur-oxidizing bacteria and effectively detoxifies porewaters. Fractofusus could thus have obtained nutrients from chemolithoautotrophic microbes/symbionts via phagocytotic processes across their fractal-like lower epithelia (Dufour and McIlroy 2017).

Fractofusus is the numerically dominant species in the lower tiers at Mistaken Point and the “Melrose Surface” assemblage (n_f = 190; 91%). The Gaussian finite mixture model-based clustering algorithms resolved one mode, interpreted as one age/size class, as the most likely number of modes both for the univariate and bivariate size-frequency distributions. Different population dynamics can generate unimodal distributions such as pioneer colonizing populations, populations with slow growth rates but high recruitment and reproduction rates, and populations showing aseasonal or continuous reproduction (Darroch et al. 2013).

As shown by the medoids projected in the PCPs illustrating the results of the polythetic cluster analysis and the optimal cluster solution (Figs. 7, 8, Tables 2, 3), there are two orientation groups. Groups under and over 180° are essentially symmetrical, as they differ by ~180° from each other. This pattern validates our initial assessment of the almost null effect of the low number of kinked specimens in the analysis on both sides of the circumference both under and over 180°. These two orientation groups show different positions with respect to the single southeasterly paleocurrent: the first orientation group is closer to a paleocurrent-orthogonal orientation (diverging at a >50° angle), whereas the second orientation group is oblique to the paleocurrent, diverging at only an ~35° angle. A near-orthogonal orientation would have maximized the aspect ratio of the organism and caused the flow of relatively unexploited water over all rangeomorph elements, leading to an increased capacity to absorb DOM from the water column. However, the feasibility of the sole reliance on osmotrophy is questionable for such large organisms (Dufour and McIlroy 2017). The high aspect ratio of current perpendicular Fractofusus would also increase the capacity for oxygen uptake from the upper surface, thereby potentially increasing diffusion through the organism into the sediment and concomitantly increasing productivity of sulfur-oxidizing bacteria below the organism. On the other hand, this orientation would have led to increased potential for current-related damage due to shear stress. Fluid separation caused by the interaction between bottom currents and Fractofusus would have led to reversed velocity profiles and eddying along the entire length of the organisms, increasing the eventual possibility of being detached and washed away by the paleocurrent (Southard 2019). This same trade-off applies to an oblique orientation relative to the paleocurrent, which would minimize aspect ratio and absorption, as distal down-current rangeomorph elements would be exposed to more exploited water. However, this position would also reduce current-related damage by confining current eddying to the distal end of the organisms. These trade-offs and orientations are suggestive of a reclining position in life (Seilacher 1992; Gehling and Narbonne 2007; Taylor et al. 2022). The two orientation groups show very small differences in morphometric traits, as shown in Tables 2 and 3 and the PCPs (Figs. 7, 8), agreeing with the Gaussian finite mixture model-based clustering algorithms that resolved one age/size class (Fig. 5).

Making a well-grounded decision on the number of clusters resulting from clustering algorithms is critical. Even though the Caliński and Harabasz pseudo F-statistic (Caliński and Harabasz 1974) traditionally shows good performance, it bears limitations, as it cannot select the one-cluster solution (Will 2016). This obstacle was solved in this approach by monothetic clustering (see Tran 2019). M-fold cross-validation showed that the one-cluster solution is statistically unlikely, as it shows a very high mean squared error (CV_k). Higher numbers of cluster solutions showed more viable CV_k and MSE values and reduced error bars compared with the one-cluster solution. However, picking the smallest CV_k (minCV rule) is an ingenuous approach, as it can result in an unrealistically high number of clusters if the error rates constantly decrease (Tran 2019) (Fig. 9). The two-cluster solution driven by orientation selected by permutation-based hypothesis test corroborated previous polythetic clustering results both under and over 180°, illustrated in the PCPs (Figs. 7–9). A three-cluster solution was also highlighted in polythetic clustering as the second most likely solution, as illustrated by the Caliński and Harabasz pseudo F-statistic, Ward's hierarchical bottom-up agglomerative clustering dendrogram (Figs. 7, 8), and the clustering algorithms on size-frequency distributions (Fig. 6). Both cross-validation and the permutation-based hypothesis test bear certain stochasticity in their algorithms. However, the variability in the outcomes of cross-validation is more relevant and may have a greater effect on its performance (Tran 2019).

Given the absence of appropriate extant analogues for comparison, clustering analyses may offer very ambiguous solutions where no correct answer can be guaranteed. The outcomes derived from clustering analysis require a certain degree of interpretation to identify structural patterns in the dataset (Tran 2019). Overall, the selected two-cluster solution regarding orientation maintains both an optimal cohesion within clusters and external isolation between clusters and is congruent with all polythetic and monothetic clustering analyses. It is also consistent with the trade-off between the nutrient/oxygen gathering strategies and the current shear stress experienced by a linear benthic organism living fixed to the seafloor while exposed to a current. Both orientation clusters offer adaptative advantages to a reclining lifestyle.

Pectinifrons abyssalis has been reported from the Avalon Peninsula (Bamforth et al. 2008), but its presence in the Catalina Dome assemblages has not been documented hitherto (Fig. 1). The species is found in association with F. misrai, although the relative abundance of both species is highly variable between fossil assemblages, and the co-occurrence is not statistically significant (Eden et al. 2022). The species has been reconstructed by Bamforth et al. (2008) as being tethered to the seafloor by a curved “pedicle rod,” with two rows of first-order branches held erect in the water column. Several specimens up to 74 cm in length have been reported. The first-order branches are typically visible on the concave side of the pedicle rod; the specimens that are subject to a change in the curvature of the rod also show a change on the side on which the first-order branches are visible (Bamforth et al. 2008).

The P. abyssalis specimens at Green Head have bimodal orientations that are incompatible with felling by a paleocurrent (see McIlroy et al. 2022). There are several specimens that have their frondlets oriented perpendicular and against both the turbidity flow and contour current direction. Moreover, there are no specimens oriented in the direction of the inferred contour current at Green Head (Bamforth et al. 2008). Near-substrate upstream-directed velocity pulses have never been reported in turbidity flows (McIlroy et al. 2022), as such felling cannot happen in the up-current direction. Variations in down-current velocity within the net downstream motion of the fluids have been described (e.g., Baas et al. 2011; Kostaschuk et al. 2018). Kelvin-Helmholtz waves at the interface with the ambient fluid at high Froude numbers include up-current directions of particle movement but not net near-bed upstream velocity pulses within the rapidly moving density current (McIlroy et al. 2022). Therefore, vortices within turbidity currents cannot be used to explain bimodal felling directions; instead, the orientation of the organisms must reflect the presence of positive and negative rheotropic growth controlled by a clear-water background current, not a density flow event (see McIlroy et al. 2022).

This rationale is further supported by P. abyssalis specimens on both the Mistaken Point Ecological Reserve (Bamforth et al. 2008) and the “Melrose Surface” that do not show any preferential orientation, suggesting that the mode of life of the organism was not significantly impacted by paleocurrent. Only one row of first-order branches is discernible in all observed specimens. As the species had been described as a filter feeder erected in the water column, capable of reaching considerable sizes, it seems unlikely that it would not have been affected by currents. Therefore, we suggest the possibility that Pectinifrons might have lain its first-order branches on the substrate to feed via a chemosymbiotic relationship with sulfur-oxidizing bacteria, as suggested for F. misrai and other Rangeomorpha (Dufour and McIlroy 2017; McIlroy et al. 2020, 2021), or perhaps even like the modern Zoothamnium (Moriyama et al. 1999; McIlroy et al. 2009).

One specimen from the “Melrose Surface” shows a change in the pedicle curvature approximately coincident with a change in the side of the pedicle with first-order branches (twisted specimen in Fig. 1B). In this specimen, first-order branches are in opposite and almost parallel orientations (N-S) on the concave sides of the pedicle rod, and the whole organism has a rough W-E orientation. It is possible that the west side had been flipped via mechanical damage and folded onto the opposite side, but could equally be due to growth. This “twisting” differs from non-idiomorphic “kinked” Fractofusus specimens that were bent and repositioned with the original side still facing down (Taylor et al. 2022).

Bradgatia is a highly variable taxon, interpreted to have been tethered to the seafloor with first-order branches erected in the water column. In the Newfoundland successions, it is considered to have been susceptible to “felling” by inferred paleocurrents (Flude and Narbonne 2008). The specimens described here do not show evidence of felling, as several are instead oriented in a close to orthogonal position and against the only measurable paleocurrent, as is the case for several specimens in all the surfaces sampled by Flude and Narbonne (2008) for both turbidity flows and contour currents. Even in surfaces where Bradgatia specimens show relatively unimodal distributions (e.g., Bishop's Cove and arguably the “E” Surface) and potential alignment to a paleocurrent (contour and turbidity currents respectively), several specimens are oriented perpendicularly and against these currents despite the high aspect ratio that they would have had should they have lived with an erect mode of life. This trend becomes more evident in the “G” and “D” surfaces (Flude and Narbonne 2008). Moreover, there are specimens oriented close to perpendicular to the paleocurrent that are close to being oriented 180° from each other at the “Melrose Surface.” These orientations make the felling precept unlikely, even if an orthogonal bottom current to the one observed was inferred, as they would have been felled into and against the current (Fig. 11) (see McIlroy et al. 2022).

The specimens show a tightly constricted arrangement of first- and higher-order branches, typical of U-shaped Bradgatia. The most basal first-order branches on the down-current side of some specimens appear to be inflated relative to those on the up-current side. The lack of current reorientation suggests that Bradgatia was most likely growing along the seafloor, rather than living erect in the water column as typically reconstructed (cf. McIlroy et al. 2021, 2022). The tightly organized branching hierarchy coupled with the basal down-current inflation suggests a functional morphology that simultaneously maximizes the surface area in contact with the seafloor and exchange with the water column. Similar lifestyles have been suggested for the rangeomorphs Beothukis mistakensis (McIlroy et al. 2020) and Fractofusus misrai (Dufour and McIlroy 2017).

Paleocurrent direction in Ediacaran assemblages of Avalonia has occasionally been inferred from the preferential orientation of purportedly felled fronds. There is thus an assumption of erect modes of life that calls into question the causal relationship between inferred paleocurrent direction and the preserved orientation of fronds. In this study, we present an integrated approach conducted on a newly described fossiliferous surface combining: (1) unequivocal physical sedimentological evidence for paleocurrent direction in the form of climbing ripple cross-lamination; (2) a novel approach based on polythetic and monothetic clustering analyses reflecting the circular nature of the recorded orientation of fronds; and (3) an intuitive representation of the data based on modified PCPs including the projection of a circular variable.

1. The Fractofusus population presents two orientation groups with respect to the southeasterly paleocurrent (102° SE). The first group is oriented >50° to the paleocurrent, whereas the second orientation group diverges only at an ~35° angle from the paleocurrent. An almost current-transverse orientation could have increased the aspect ratio of the specimen, maximizing the ability to absorb DOM as well as the capacity to diffuse oxygen from the dorsal to the ventral surface. This dorsoventrally transported oxygen could have increased productivity of chemolithoautotrophic bacteria and reduced buildup of sulfides. On the other hand, a more oblique orientation would have reduced specimen aspect ratio as well as reducing the effect of lee-side current eddying and erosion, and therefore the eventual possibility of being washed away by the paleocurrent.

2. Pectinifrons specimens do not show any preferential orientation, suggestive of a mode of life not heavily influenced by paleocurrent direction. This notion challenges the previous reconstruction of Pectinifrons as a filter feeder but is compatible with a reclining mode of life with first-order branches in permanent contact with the seafloor. Alternative feeding strategies via a chemosymbiotic relationship with sulfur-oxidizing bacteria should be explored.

3. Bradgatia specimens show oblique to current-transverse and upstream orientations to the only measurable paleocurrent. The most basal first-order branches on the down-current side appear to be inflated compared with the respective first-order branches on the up-current side. These lines of evidence suggest that Bradgatia grew along the seafloor, responding to the current in a rheotropic manner, rather than being erect in the water column. The suggested reclining mode of life is indicative of a functional morphology that simultaneously maximizes surface area in contact with the seafloor and exchange with the water column suggestive of chemosymbiosis.

The analytical method herein presented offers great potential to explore datasets consisting of fronds with unidirectional directions with respect to the inferred paleocurrents, particularly at sites with independent current indicators. The findings of this work support recent interpretations of some rangeomorphs as being reclining (McIlroy et al. 2021, 2022) showing a rheotropic response to the paleocurrent (cf. rheotropic bryozoan growth in Ryland et al. [1977]), which will require a closer look at paleoecological datasets that have generally assumed an erect mode of life for all the Rangeomorpha excluding Fractofusus (i.e., Clapham and Narbonne 2002; Clapham et al. 2003; Liu 2011; Hoyal Cuthill and Conway Morris 2014; Antcliffe et al. 2015; Mitchell and Kenchington 2018; Vixseboxse et al. 2021).

We acknowledge the constructive criticism by Emily Mitchell and an anonymous reviewer and associate editor A. Tomasovych, as well as the support from E. Samson and Discovery UNESCO Global Geopark. Fieldwork and fossil casting in Newfoundland were conducted under permit from the Parks and Natural Areas Division of the Government of Newfoundland and Labrador. This project was supported by a Discovery Grant and Discovery Accelerator Supplement from the Natural Sciences and Engineering Research Council to D.M.

The authors declare no competing interests.

Data available from the Dryad Digital Repository: https://doi.org/10.5061/dryad.dbrv15f4x.