‘Conga lines’ of Ediacaran fronds: insights into the reproductive biology of early metazoans

Late Ediacaran strata from Newfoundland, Canada (~574–560 Ma) document near-census palaeocommunities of some of the earliest metazoans. Such preservation enables reproductive strategies to be inferred from the spatial distribution of populations of fossilized benthic organisms, previously revealing the existence of both propagule and stoloniferous reproductive modes among Ediacaran frondose taxa. Here, we describe ‘conga lines’: linear arrangements of more than three closely spaced fossil specimens. We calculate probabilistic models of point maps of 13 fossil-bearing bedding surfaces and show that four surfaces contain conga lines that are not the result of chance alignments. We then test whether these features could result from passive pelagic propagules settling in the lee of an existing frond, using computational fluid dynamics and discrete phase modelling. Under Ediacaran palaeoenvironmental conditions, preferential leeside settlement at the spatial scale of the conga lines is unlikely. We therefore conclude that these features are novel and do not reflect previously described reproductive strategies employed by Ediacaran organisms, suggesting the use of mixed reproductive strategies in the earliest animals. Such strategies enabled Ediacaran frondose taxa to act as reproductive generalists and may be an important facet of early metazoan evolution.


Introduction
Metazoans exhibit a remarkable range of reproductive strategies, each of which has varied costs and benefits [1].Tradeoffs between these costs and benefits can result in mixed reproductive strategies, as exhibited by extant members of some of the earliest diverging metazoan phyla, including poriferans and cnidarians [1,2].Mixed reproductive strategies may manifest in a number of ways.For example, a single organism may employ multiple reproductive modes either simultaneously, or throughout its lifespan [3,4].Alternatively (and often in addition), reproductive mode may be controlled by phenotypic plasticity [1,4]: the production of varied phenotypes from a genotype in response to different environmental cues.The abundant presence of mixed reproductive strategies in the earliest diverging phyla may suggest that the first animals also exhibited mixed reproductive strategies.
Rangeomorphs have been suggested to be capable of reproducing by stolons, fragments, buds and pelagic propagules (figure 1) on the basis of spatial distributions across multiple scales, as well as direct evidence from fossil impressions [11,[26][27][28][29][30].Different reproductive strategies manifest as different spatial patterns of sessile organisms.The use of methods such as spatial point pattern analysis (SPPA) can disentangle reproductive processes from environmental and interactive influences.Such analyses have been used to describe the patterns produced by known reproductive processes in extant marine sessile invertebrates (e.g.[31,32]).The near census preservation of Avalonian palaeocommunities permits the application of these techniques to elucidate rangeomorph reproductive biology [11,25,26,33,34].For example, the rangeomorph Fractofusus has been statistically demonstrated to form small, isotropic (i.e.non-directional) and hierarchical clusters, comprising smaller specimens clustered around larger specimens, which were inferred to arise from stoloniferous reproduction [11].
More recently, observed filamentous connections between rangeomorphs were interpreted to be stoloniferous, on the basis of the filaments terminating at the bases of individual specimens, and their large size (sometimes extending over a metre) being incompatible with fungal or bacterial affinities [27].Evidence for pelagic propagules derives from the wide palaeogeographic distribution of taxa such as Charnia [26,28], and current-influenced spatial distributions [11], but the nature of the propagules remains unresolved.Larvae, gametes, embryos and microscopic buds or fragments are all possible forms of pelagic propagules given existing data and rangeomorph phylogenetic placement [19,35].Putative macroscopic buds and fragments have been suggested to occur in Culmofrons and Avalofractus  royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 231601 [29,30], although the approximately fractal branching [36] of rangeomorphs may make it hard to distinguish such fragments from whole organisms (e.g.[29]).

'Conga lines' and explanations for aligned specimens
On the H5 fossil-bearing bedding plane in the Discovery Geopark, Newfoundland, Canada (figure 2), we observed examples of frondose specimens forming closely spaced, linear arrangements of 3-6 specimens, which are here termed 'conga lines' (figure 3).In two observed cases, conga lines appear to have a frond at one end of the line that is substantially larger than the others (figure 3).Such aligned arrangements of specimens have not previously been reported from Ediacaran fossil assemblages.The spatial arrangement of sessile marine invertebrates reflects the interplay between abiotic factors, reproductive processes and interactions such as competition [38].Linear arrangements of sessile organisms, such as those that form the conga lines, could only feasibly arise from a much narrower set of such factors.First, they could result from random chance, particularly if specimens occur in high densities or are aggregated due to other processes, due to the increased probability that they will seemingly align.Second, the interaction of a directional abiotic factor, such as a current, with reproductive processes could result in a linear pattern of recruitment, which would be influenced by the nature of the reproductive propagules.We consider two end members of reproductive propagules based on their dispersal capacity: pelagic and philopatric (sensu [39]).Here, pelagic is used to describe passive propagules that are neutrally or positively buoyant, and those capable of active swimming, with long-distance (teleplanic or anchiplanic) dispersal capacity.Pelagic propagules may include larvae, gametes, embryos, and microscopic buds or fragments, with larvae typically developing planktonically or demersally [40].Conversely, philopatric (remaining within centimetres to a couple of metres of their parent) propagules are negatively buoyant (e.g.[41]) or capable of active swimming, have an aplanic dispersal capacity (immediate settling), and may demonstrate benthic (including parental) development [40,42].Philopatric propagules include macroscopic buds or fragments, some embryos and larvae (often brooded or crawling [40,43]).Flow patterns around an obstacle may lead to increased settling of passive pelagic propagules in a linear fashion, or alternatively, the sequential production of philopatric propagules may interact with a current to create linear features in the lee of the parent organism.Third, linear arrangements of sessile organisms may be purely biological, e.g. a consequence of a non-branching, reproductive runner-like stolon connecting multiple organisms in sequence.
In this study, we seek to determine whether the newly documented Ediacaran conga lines are the consequence of chance alignments by calculating probabilistic models of point maps of 13 fossil-bearing surfaces from Newfoundland, Canada.We then test whether conga lines could be formed from passive pelagic propagule settlement in the lee of an obstacle, using computational fluid dynamics and particle tracking.Following these tests, we then consider the most likely mode of formation of the conga lines and the consequences for early metazoan evolution.

Data collection
We analysed point maps of portions of 13 late Ediacaran fossiliferous surfaces from Newfoundland, Canada, comprising two surfaces newly mapped in this study (H5 and MUN from the Discovery Geopark), supplemented by 11 further surfaces compiled by [44] incorporating data from [25] (Bishop's Cove in the Spaniard's Bay area; Goldmine, H14 and H26 in the Discovery Geopark; Brasier (also known as BR5), Clapham's Pigeon Cove (CPC), 'D' Surface, 'E' Surface, 'G' Surface, Pizzeria and Lower Mistaken Point (LMP) in Mistaken Point Ecological Reserve; figure 2).Together, these point maps document 16 002 specimens over 425 m 2 (electronic supplementary material, table S1).This dataset excluded available published data from surfaces that are unsuitable for spatial analyses due to tectonic deformation (Shingle Head), non-typical sedimentology (UIC), small areal extent (Bristy Cove) or ambiguity surrounding the contemporary nature of discoidal fossils (Bed B).
Well-preserved areas of the H5 and MUN surfaces [37,45] were moulded in the field using silicone rubber.A total of eight moulds encompassing 4.11 m 2 were created for H5, and five moulds encompassing 3.91 m 2 were created for MUN.Jesmonite casts were subsequently made from these moulds (housed in the Sedgwick Museum of Earth Sciences, Cambridge, UK; H5: CAMSM X 50342.1-9,MUN: CAMSM X 50340.1-6).Photogrammetry of casts under controlled lighting conditions was used to produce high-resolution two-dimensional orthomosaics, which were then scaled and aligned using photogrammetric maps of the entire surface taken in the field.The orthomosaics were subsequently checked against LiDAR scans for accuracy (using Agisoft Metashape Professional 1.8.4., cf.[25]).Vector maps were then produced from these orthomosaics in Inkscape 1.0.1, documenting specimen spatial position, dimensions (frond length, frond width, stem length and stem width), orientation and taxonomic identification (cf.[25]), and processed using a custom R script (https://github.com/kmdelahooke/dex) to generate point maps of the two surfaces (electronic supplementary material, figures S1 and S2).To correct for tectonic deformation, the point coordinates were retrodeformed by returning ovate discoidal fossils to circles (cf.[11,12]).To investigate the influence of retrodeformation and the presence of large cracks, small holes and eroded patches on our analyses, additional point maps of the H5 surface were created for retrodeformed data, non-retrodeformed data, and a case where all  small holes and eroded areas of the surface were masked (i.e.removed), to assess the sensitivity of our methods to small-scale, patchy erosion.

Probabilistic modelling
To investigate whether the conga lines result from chance alignments, we used both analytical calculations and numerical probabilistic modelling.First, we calculated the probability of finding the observed number of conga lines on each surface analytically (i.e. by calculating the exact solution using derived equations, §2.2.3).Second, we verified the analytical calculations using numerical Monte-Carlo simulations of the H5 surface ( §2.2.4).Probability analyses were performed using custom code in R 4.2.2 (R Core Team 2022), using the spatstat package [46].All codes are available on Github: https:// github.com/kmdelahooke/Conga-Lines.

Spatial point process analysis
The spatial distribution of points (specimens) on each surface was characterized prior to probabilistic calculations, as the background distribution of points naturally affects the likelihood that chance alignments will be found.At the simplest level, points can be randomly distributed, aggregated (points closer together than random) or segregated (points further apart than random; figure 4a-c).Intuitively, if the density of points for any of these patterns increases, the probability of chance alignments also increases.If the pattern is aggregated, it would be more likely that one point is proximal to another, but the inverse is true if points are highly segregated.Second-order summary functions, such as the K-function or pair correlation function (PCF), can be used to quantify these patterns (aggregation, segregation or a random distribution) across a range of spatial scales [47].
Point process models are statistical models that generate point patterns within an observation window.In order to allow the point pattern to be characterized by a series of parameters, these models can be fitted to a summary function of a given point pattern, such as density (first-order intensity function) or the K-function or PCF (second-order functions [38,47]).Most simply, the points could be randomly distributed: complete spatial randomness (CSR, figure 4a) is modelled by a homogeneous Poisson process parameterized by density (λ), equivalent to the number of points (n) (here specimens) over a given area (A) (here mapped area).Alternatively, points could be aggregated, which can be modelled by processes including a Thomas cluster process that can describe reproductive aggregations (figure 4b), or heterogeneous Poisson processes, which best describe aggregation due to shared habitat associations [47].Thomas cluster processes are parameterized by the density of parent points (κ) (themselves following a homogenous Poisson distribution), the mean number of offspring points in the cluster (μ) and the standard deviation of the distance of the offspring points to their parent point (σ).Conversely, points could be segregated, which can be modelled by a suite of processes including the Gibbs process (figure 4c).Many Avalonian populations show clustered distributions, of which Thomas cluster (TC) and double Thomas cluster (DTC) processes dominate [25].

Background distribution model fitting and evaluation
Both homogenous Poisson and Thomas cluster processes were fitted to point patterns for each surface using the package spatstat [46] as follows: (1) Homogeneous Poisson process (CSR) models were fitted to the density (intensity) function of the point pattern using the maximum likelihood method [48], generating parameter λ. (2) Single TC models were fitted to the PCF of the point pattern using the minimum contrast method, generating parameters κ, μ, σ [49].(3) The best process to describe a point pattern was determined using Diggle's goodness-of-fit test [47,50].( 4) The point process model with the highest goodness-of-fit (p d ) value was used as the background distribution when calculating the probability of chance alignments (e.g.[25,47]).

Predicted number of congas-analytical method
The predicted number of chance alignments was determined analytically assuming both CSR and Thomas clustered background distributions.The predicted number of chance three-point alignments, of maximum spacing r mm and maximum angle θ °, was derived from Ripley's K-function (electronic supplementary material, §2.1), which describes the likelihood of finding another point within a certain radius of a given point of the point pattern [51].
Given a CSR background distribution, where λ is the density of points and A is the mapped area, the predicted number of three-point alignments on a surface (E pois ) is, (2.1) Using a Thomas clustered distribution, where μ, κ, σ are the fitted TC parameters (see §2.2.1), and A is the mapped area, the predicted number of three-point alignments on a surface (E TC ) is, (2.2)

Predicted number of congas-numerical method
The analytic models were validated numerically on the H5 surface as follows: (1) CSR and TC models were fitted, as outlined in §2.2.2 (steps 1 and 2), to the entire surface point pattern.
(2) Ten thousand Monte-Carlo simulations of both random and clustered point processes were generated using spatstat, using the parameters from fitted homogeneous Poisson (λ = 140 m −2 ) and TC models (κ = 18 m −2 , μ = 6, σ = 72.63mm), and within a window defined by the H5 surface outline.(3) For each simulation, we ran an algorithm (https://github.com/kmdelahooke/Conga-Lines) to detect linear features consisting of at least three aligned points (parameterized by r and θ), and the number of three-point alignments found at that iteration was recorded.This algorithm visits each point in the point pattern (p 1 ) and finds any points that are within a radius r mm (p 2 ).For each of these points, any points within r mm (p 3 ) are detected (figure 4d).All prior points in the sequence are excluded from the search.The bearing (angle) between p 1 and p 2 (b) as well as p 2 and p 3 (b′) is calculated (figure 4d).If the bearings are within a certain angle θ, such that | b′ − b |< θ, then the coordinates of each point were recorded.The number of alignments found for each simulated point pattern was returned.(4) Step 3 was repeated, while holding all other parameters at a control on the surface (electronic supplementary material, tables S2 and S3), for the following range of values: (i) r : 10, 12, 15 and 20 mm (control: 15 mm) encompassing the variation in observed spacing of conga lines on the H5 surface. (ii) θ : 20°, 25°, 30° and 45° (control: 30°), encompassing strong to weak alignments.The relationship between the predicted number of three-point alignments and each parameter could then be compared between both analytical and numerical methods (see §2.2.6).

Probability distribution
A probability distribution for each set of parameters was drawn from the Monte-Carlo simulations by plotting the number of times a simulation returned a given number of three-point alignments.This probability distribution follows a Poisson distribution, , where x is the observed number and E D is the predicted number (E pois or E TC ) of three-point alignments.

Comparison of observed versus predicted numbers of conga lines
To calculate the probability of seeing the observed number of conga lines found on each surface by random chance, the following steps were taken: (1) The conga-line-finding algorithm was applied to the point maps of each fossil surface for a range of radius r between 5 and 50 mm, at 1 mm increments, to find the observed number of conga lines on the surface.(2) CSR and TC models were fitted to the surface point pattern, following §2.2.2 (Steps 1 and 2).
(3) Model fit was evaluated using Diggle's goodness-of-fit test and the best-fit model was found.(4) Using the best-fit background distribution, the predicted number of congas for that surface was calculated using equations (2.1) and (2.2) for the range of r. (5) The probability of the observed number of congas found (Step 1) given the predicted number of congas (Step 4), for each value of r, was calculated using the Poisson distribution ( §2.2.4). (6)The range of r in which the probability was less than 0.05 was noted, in addition to whether more or fewer congas were found than predicted.

Computational fluid dynamics
Computational fluid dynamics (CFD) is a tool for simulating fluid flow under varying conditions by solving a suite of governing equations that describe the flow.CFD analyses are increasingly used in Ediacaran research to describe how fluids interact with the morphological characteristics of taxa, such as Tribrachidium [52], Parvancorina [53], Ernietta [54,55], Arkarua [56], Pteridinium [57] and Pectinifrons [58], permitting inference of feeding mode and mobility.
Here, we use established CFD methods alongside discrete phase modelling (DPM).DPM is a particle tracking method that enables the dispersal and settlement of propagules to be modelled [59,60], to test whether the conga lines on the H5 surface could be formed from passive propagule settling in the lee of an existing frond.

Steady-state flow fields
CFD was implemented within Ansys Fluent (2021 R2).The frond was modelled as a rectangle in the centre of a 2 m × 1 m, two-dimensional computational domain, which was large enough to accommodate the extent of downstream eddies and permit the development of a boundary layer [61].The floor and obstacle were defined by no-slip boundaries, the top of the domain was defined by a free-slip boundary, and the inlet and outlet were defined at opposing ends of the domain.The inlet had a uniform flow profile and the outlet had a pressure equal to the inlet.A quadrilateral dominant mesh royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 231601 with a maximum element size of 5 mm (89 711 elements) was used.Inflation layers, consisting of 24 maximum layers and a growth rate of 1.2, were used to resolve the boundary layer close to the wall, generating a velocity profile consistent with the law of the wall (electronic supplementary material, figure S3).A mesh sensitivity analysis, where the element size was varied between 2 and 10 mm (524 080 to 24 938 elements), was carried out to determine the impact of the mesh parameters on eddy length.A characteristic length of 1 m, at an inlet velocity of 10 cm s −1 yielded a Reynold's number of 65 478, indicative of a turbulent flow regime.Steady-state RANS simulations were run primarily using the k-ω SST turbulence model.The Reynold's stress model is more accurate but does not converge for many target velocities, thus the k-ω SST was used to balance the greatest accuracy with computational tractability (e.g.[55]).Additional parameters used to describe the flow and explanations can be found in electronic supplementary material, §2.3 and tables S4, S5.
The following simulations were carried out: (1) Two-dimensional CFD simulations using the k-ω SST turbulence model were replicated over: (i) A range of frond heights (1, 2, 4, 6, 8 and 10 cm), corresponding to the size distribution of frondose specimens on the H5 surface.(ii) A range of velocities (1, 5, 10 and 15 cm s −1 ), that encompass and exceed the range of velocities typical in modern deep-water settings (1-10 cm s −1 [62]).
(2) A three-dimensional simulation (10 cm s −1 , 5 × 0.5 × 2 cm frond), which allows features such as horseshoe vortices to be modelled, was run to verify the two-dimensional result.(3) Simulations using the RSM turbulence model were run using a frond height of 5 cm and velocities of 1, 5 and 15 cm s −1 to compare to the k-ω SST turbulence model.

Discrete phase modelling
To model the distributions of passive propagules in the simulated flow discrete phase modelling (DPM), a Lagrangian particle tracking technique, was used within Ansys Fluent (2021 R2).The discrete phase (here the propagule particle, parameters found in electronic supplementary material, table S6) is tracked as it interacts with the continuous fluid phase.Here, an uncoupled approach was employed, where the discrete phase does not affect the continuous phase.Steady-state particle tracking was implemented, where the solution of the continuous phase flow was used to calculate the trajectories of the particle.A specified length scale determines the distance a particle travels between trajectory updates, acting as a pseudo-time step.A surface particle injection was defined at the inlet with a flow rate of 1 × 10 −20 m −1 .To account for the turbulent dispersion of particles, the discrete random walk model was used with random eddy lifetime (electronic supplementary material, §2.3 and table S6).When a particle hits the floor, it will settle if the local shear stress is below a critical value (τ crit ), otherwise, it will be reflected.This condition was specified using a custom function (https:// github.com/kmdelahooke/DPM),where τ crit = 0.1 Pa, as has been found experimentally for weakly adherent larvae [63].If a particle settles, its position and velocity are recorded.This information is used to assess how the probability of passive pelagic propagule settlement changes along the x-axis in each simulation.DPM simulations were replicated for the range of the frond heights (1, 2, 4, 6, 8 and 10 cm) and velocities (1, 5, 10 and 15 cm s −1 ) used in the CFD simulations, as well as for a range of particle sizes (50 μm, 200 μm, 500 μm, 1 mm and 5 mm), encompassing the range of sizes of extant parenchymella (demosponge) and planula (cnidarian) larvae [64,65].

3.1.
Were conga lines formed by chance?

Comparison of analytical and numerical results
The relationship of the predicted number of three-point alignments given a CSR background distribution (E pois ) with the spacing (r) and angle (θ) between three-point alignments and the density of the point pattern (λ), as well as the relationship of the predicted number of three-point alignments given a TC background distribution (E TC ) with r, θ and TC parameters κ (density of the parent points), μ (mean points per cluster) and σ (standard deviation of distance of an offspring point to its parent), as determined by both analytical and numeric methods, was compared between these methods (electronic supplementary material, figures S4 and S5).The analytical and numerical data are consistent with one another, with slight discrepancies likely due to the stochastic nature of numerical methods and possibly unaccounted for edge effects [66,67], indicating the validity of analytical methods (electronic supplementary material, figures S4 and S5).

Conga lines-are they real and where are they found?
Of the mapped areas, four surfaces have significantly more three-point alignments than are predicted by random chance (H5, MUN, Brasier, CPC), while nine do not (table 1; figure 5; electronic supplementary material, figures S6 and S7; n.b. additional apparent conga lines were observed outside the mapped area on the LMP surface, electronic supplementary material, figure S8).A large number of observed conga lines are composed of effaced (i.e.lacking well-preserved morphological features necessary for taxonomic assignment) fronds on MUN, H5 and Brasier, which comprise a high proportion of overall fossils on each of these surfaces (electronic supplementary material, table S1).The best-preserved conga lines on MUN and H5 both comprise Charnia spp.(figure 6a,b,e), while on the CPC surface, conga lines are only expressed as linear arrangements of circular raised structures inferred as holdfast discs.One of these conga lines contains a faint filamentous structure running through four of these features (figure 6c).These physical observations and probabilistic approaches together indicate that the conga lines likely reflect reproductive processes and are unlikely to have arisen by chance alone, and that they appear to reflect original features of these fossilized bedding plane assemblages.

3.2.
Were conga lines formed by waterborne propagule settling?

Computational fluid dynamics
Contour plots of flow parameters show that velocity, pressure and wall shear stress drop behind the obstacle (see also electronic supplementary material, figure S9).Two small eddies form either side of the frond, as well as a larger one in the lee (figure 7a).These eddies can also be seen as negative stream-wise velocity just above the floor (electronic supplementary material, figure S10).Increasing the height of the obstacle/frond increases the size of the larger lee-side eddy, following a linear relationship (adjusted R 2 = 0.999, p = 0.0002, F 1,2 = 5130; electronic supplementary material, figure S11), thus larger fronds would lead to a larger area of propagule settlement behind them, with longer predicted conga lines.Increasing the velocity also increases the length of the eddy, but nonlinearly, again increasing the predicted length of conga lines.The length of the region of lowest velocity (which is inversely proportional to the probability of propagule settlement) is always much larger than the observed length of Ediacaran conga lines, across frond height and a range of plausible velocities (figure 7b).At finer mesh sizes, the length of this zone increases (electronic supplementary material, figure S12), therefore increasing support for this conclusion.When comparing the CFD outputs from different turbulence models, the k-ω SST model (used here) predicts a slightly longer eddy, and thus an even longer conga line behind the frond than the RSM model (electronic supplementary material, figures S13 and S14).Therefore, CFD analyses suggest that for a frond of height <10 cm, the length of the observed conga lines (<4 cm) is much smaller than would be predicted (>8 cm) across all examined flow velocities (figure 7b), suggesting that passive pelagic propagule settlement is unlikely to be their formative process.

Discrete phase modelling
Discrete phase modelling (DPM) was used to determine the frequency of particle settlement under different velocity conditions relative to the position of the frond.At low inlet velocities (<10 cm s −1 ), the shear stress is always less than the critical shear velocity (0.1 Pa), thus the frequency of larval   settlement is equally high both in front of and behind the frond (figure 7c).At high inlet velocities (>50 cm s −1 ), the shear stress is greater than the critical shear stress, except behind the frond.Consequently, the anticipated frequency of propagule settlement increases behind the frond (frequency of propagule settlement is four times higher ~0.5 m behind the frond compared with immediately before; see figure 7c).Under conditions typical of the inferred deep-marine setting of H5 [37,68], the flow velocities are considered to have been slow (i.e.<10 cm s −1 [62,69]) and hence shear stresses would have been too low to ever re-suspend or prevent the settlement of propagules.When DPM was repeated with different particle sizes, only the smallest particle size (500 μm diameter), showed preferential settling in the lee of the frond (electronic supplementary material, figure S15).Thus, the presence of an obstacle provides no critical benefit to propagule settlement.

Discussion
Our probabilistic models have shown that the conga lines observed on four late Ediacaran fossil-bearing surfaces in Newfoundland were not the result of chance alignments (table 1; figure 6).CFD and DPM outputs show the region of lowest velocity in the lee of the frond is not comparable to the length of the conga lines observed on H5, such that propagule settlement is not expected to increase in the lee of a frond under low-velocity environmental conditions (<10 cm s −1 ; figure 7c).However, it should be noted that the steady-state analyses utilized herein do not consider the impact of fine-scale, ephemeral turbulent structures on propagule settlement.The steady-state results are consistent with the observation that analogous modern features, such as sedimentary obstacle marks [62,70], and preferential larval settlement behind kelp or coral [71,72], are restricted to higher energy environments such as wave-dominated coastal regions and flood events, which provide a sufficient gradient in local shear stresses and velocities to cause preferential deposition behind an obstacle [63,70].The hypothesis that the observed conga lines are the result of differential passive pelagic propagule settlement around an obstacle can therefore be rejected.Consequently, conga lines are likely indicative of a dispersal-limited reproductive strategy not previously described in rangeomorphs: either through the interaction between philopatric propagules with a current, or a runner-like stolon.
Philopatric propagules are prevalent across many groups of sessile marine invertebrates, including poriferans and cnidarians [42,43,[73][74][75][76].In the presence of a current, it is feasible that there could be preferential lee-side recruitment of the propagule relative to its parent, thereby producing a linear feature.Indeed, many conga lines have an axis orthogonal to the felling direction of the fronds (figure 6), perhaps suggesting the influence of a contour current.However, physical evidence for rangeomorph production of buds or fragments is limited [27,30 compare 11], and the preservation of larvae has not been reported in Avalonian strata.
Alternatively, the conga lines could be the consequence of a runner-like stolon, such as those found in hydrozoans [77], anthozoans [78] and entoprocts [79], as well as in algae [80] and terrestrial plants [81].Stoloniferous reproduction has previously been described in Avalonian communities [11,27], but not in a closely spaced runner-like form.The patterning of conga lines is highly directional, in contrast to the isotropic, network-like pattern of stolons inferred for Fractofusus [11].Similarly, the filamentous connections previously observed between frondose taxa [27] lack the directionality and the close spacing of the conga lines described here.Such variability in stolon patterning mirrors the variability seen in extant stoloniferous metazoans, whose stoloniferous connections are often described on a spectrum between 'runner-like' (larger spacing between zooids and rare stolon branching) and 'sheet-like' (closely packed zooids and extensive stolon branching) morphs [82].Ontogeny can control this stolon patterning, e.g.Hydractinia colonies transition from more runner-like to more sheet-like as they develop [83].Runner-like forms are also typically associated with nutrient-poor, resource heterogeneous and disturbed modern environments [80,84,85].For example, the soft coral Efflatounaria shows rapid stolon extension, creating a more runner-like form, as a plastic response to environmental disturbance [84].Such stolon forms, with minimal branching and spaced-out zooids, enable efficient colonization of the bare substrate as well as escape from poor environmental conditions, a 'guerrilla' strategy that is akin to an optimal foraging strategy but for sessile organisms [86,87].By contrast, sheet-like forms, with extensive branching and close-spaced zooids, reflect a 'phalanx' strategy, allowing organisms to capitalize on amenable conditions in their local vicinity [86].
The close-spacing of individual fronds in the conga lines is at odds with the adaptive benefit of a 'guerrilla' strategy for runner-like stolon, however, this may suggest that they may reflect an alternative strategy, or indeed no strategy at all.Furthermore, a runner-like stolon is not preserved in any conga lines bearing fronds detected in this study, with the exception of the filamentous structure running through inferred holdfast structures on Clapham's Pigeon Cove surface (figure 6c).This lack of preservation is also the case for inferred Fractofusus stolons [11] and may reflect either a taphonomic absence or a transient existence of the original biological stolon structure.In contrast to most modern benthos, the Ediacaran seafloor was often covered by a microbial mat, preserved as disparate textures on bedding planes [88,89].A mat provides a possible explanation for a preservational bias: if the stolon was embedded in the mat (perhaps for stability) or overgrown by later mat growth, then it may lie beneath the plane of preservation [27].Alternatively, the stolon could be ephemeral and no longer be present at the time the surface was preserved.Although most extant metazoan stolons are retained throughout the lifespan of the colony and provide a means of nutrient transfer, some stolons may disappear after a daughter colony is established (e.g.[90,91]).The adaptive benefit of a runner-like stolon may differ from modern analogues, but the Ediacaran conga lines are consistent with a stoloniferous formative mechanism.
Conga lines represent an additional strategy of rangeomorph reproduction, augmenting previously described propagules [11,26,28], putative buds and fragments [29,30] and previously described forms of stolons [11,27].These different reproductive modes may be solely associated with specific taxa, or may instead reflect mixed reproductive strategies within individual taxa capable of employing multiple modes simultaneously [3], reproductive plasticity [1,2] or a combination of both options.On the basis of available evidence, conga lines do not appear to represent a reproductive mode that defines a single taxon.There is no clear signal regarding the constituent taxa of conga lines because a large proportion of specimens are effaced (electronic supplementary material, table S1), but both the MUN and H5 surfaces contain seemingly monospecific conga lines formed of Charnia spp.(figure 6a,b,e).However, not all Charnia or indeed effaced frond specimens on a given surface are associated with conga lines.Charnia has been shown to form small clusters on Bed B, Charnwood, UK [25], and some Charnia specimens are associated with long filamentous structures on bed LC6 in Newfoundland [27], demonstrating that some of the taxa forming conga lines were also capable of other reproductive strategies.Based on this evidence, we cannot yet determine whether conga lines reflect reproductive plasticity or are one of multiple simultaneous reproductive modes employed by individual organisms.
Mixed reproductive strategies, whether plastic or simultaneous, have an adaptive benefit in variable environments.Plastic strategies conditional on environmental cues produce, in theory, an optimal phenotype for the conditions [92], and are adaptive when the timescale of environmental variability is greater than the lifetime of the organism [93].In contrast, stochastically mediated plasticity or multiple synchronous reproductive modes can be viewed as bet-hedging strategies [92].Such strategies produce a suboptimal phenotype, with lower instantaneous fitness but greater average fitness, thereby creating an evolutionary stable strategy.Modelling of multiple reproductive strategies in a fluctuating environment has found that the coexistence of two or more reproductive modes is favoured in all situations except when the local extinction risk is at the extremes (i.e.within variable environments), and thus multiple strategies can be viewed as indicative of variable environments [80].Therefore, mixed rangeomorph reproductive strategies may reflect an adaption to environmental variability, which could be spatial or temporal.
Limited data exist regarding spatial environmental variability within late Ediacaran deep-marine depositional environments.The sedimentology of bedding planes in instances where a fossiliferous bedding plane outcrops multiple times over several kilometres is largely invariant [94].Small-scale spatial environmental heterogeneity may exist, such as where a microbial mat or a local nutrient source is present [88,89], though previous work has shown a limited impact of inferred small-scale heterogeneities on community dynamics [25,33].Temporal environmental variation may have resulted from episodic disturbances by the turbidites and ash falls that characterize the succession [9,95,96] or from varying nutrient levels, temperatures or redox conditions, for which there is currently little bed-scale information from Avalonian settings.Although the frequency of disturbance from ash or turbidite events is not known, modern analogues have a periodicity of 10-200 years [97,98], which are likely to typically exceed rangeomorph lifespans.As such, it is plausible that temporal environmental variability may be a driver of mixed reproductive strategies, in particular reproductive plasticity, in rangeomorph populations.
Using mixed reproductive strategies would have enabled rangeomorphs to act as reproductive generalists under multiple conditions of environmental variation, allowing them to colonize large areas of the seafloor and achieve a global palaeogeographic distribution [11,28,99], as well as withstand environmental perturbations [96].Such reproductive strategies perhaps contributed to their persistence in late Ediacaran ecosystems for tens of millions of years [9,100].The use of mixed reproductive strategies by rangeomorphs demonstrates the ancient origins of this strategy, which characterizes the success of many modern marine invertebrate groups [4] and may have played a role in facilitating the initial success and global expansion of early macroscopic metazoans.

Conclusion
We describe novel 'conga lines', alignments of closely spaced frondose macro-organisms, from four bedding plane surfaces in the late Ediacaran of Newfoundland.Using probabilistic modelling, we show that these features did not form by random chance.Computational fluid dynamics, in particular discrete phase modelling, demonstrate that they did not result from passive pelagic propagule settling.We thus conclude that the conga lines were most likely formed by closely spaced, runner-like stolons, though current-influenced philopatric propagules cannot be entirely dismissed.Conga lines provide a new line of evidence in support of mixed rangeomorph reproductive strategies, in particular dispersal limitation associated with clonal organisms.Such strategies reflect evolutionary trade-offs in the earliest animal communities, perhaps in response to environmental instability, and may have contributed to the persistence and dominance of rangeomorphs in certain late Ediacaran marine settings.Data accessibility.Data can be found at https://doi.org/10.6084/m9.figshare.25314391and https://doi.org/10.6084/m9.figshare.24270499.

Ethics.
Electronic supplementary material is available online [101].
Declaration of AI use.

Figure 3 .
Figure 3. 'Conga lines' of rangeomorph fronds on the H5 surface, Discovery Geopark, Newfoundland, Canada, with larger frond on the right.(a, b) Scale bars = 1 cm.(c) Schematic illustration of specimens in (b), where the blue frond is larger than the orange fronds.

Figure 4 .
Figure 4. Schematic examples of different types of point processes.(a) Complete spatial randomness (CSR), modelled by a homogenous Poisson process.(b) Aggregation, here modelled by a Thomas cluster process.(c) Segregation, here modelled by a Gibbs process.(d) Visual representation of the method used to calculate alignments of three points (p 1 , p 2 , p 3 ), where r is the maximum spacing between points, θ is the maximum angle between points and b is the bearing between those points, such that | b′ − b |< θ.

Figure 5 .
Figure 5. (a) Observed (actual) versus predicted (expected) number of conga lines on the H5 fossil surface (θ = 30).Input data have been retrodeformed.(b) Probability of the actual number of conga lines with fronds spaced at a given radius (r) on the H5 surface (retrodeformed), with a Thomas Cluster background distribution.(c) Probability distribution of the number of conga lines on the H5 surface for a range of radius spacings.Points result from Monte-Carlo simulations; lines follow a Poisson distribution.

Figure 6 .
Figure 6.'Conga lines' of aligned frondose organisms on late Ediacaran fossil-bearing surfaces from Newfoundland, Canada.(a, b) H5, Bonavista Peninsula.(c) Clapham's Pigeon Cove, MPER, with a faint linear structure seemingly present connecting the raised circular impressions.(d) Effaced specimens on the Brasier/BR5 surface, MPER.(e) MUN surface, Bonavista Peninsula.Photographs in (a), (b) and (e) were taken of replica casts housed in the Sedgwick Museum of Earth Sciences, University of Cambridge, UK.Scale bars = 1 cm.

8 Figure 7 .
Figure 7. Results of the CFD and DPM analyses.(a) Contour plot of velocity magnitude in the lee of an obstacle (white vertical line).Free stream velocity is 10 cm s -1 , the obstacle is at 1 m and is 5 cm × 0.5 cm in height and width respectively.(b) Length of region of lowest velocity for a given frond height and flow velocity (blue), as compared to the observed length of H5 conga lines (orange).(c) Frequency of particle settlement along the floor under different velocity regimes.Obstacle is at 1 m.The frequency of settlement only increases after the obstacle when velocity is greater than 0.5 m s -1 .
Fieldwork and moulding were carried in the Discovery UNESCO Global Geopark in 2018 and 2021, under permits granted by the Provincial Archaeology office of the Department of Tourism, Culture, Arts and Recreation, Government of Newfoundland and Labrador to E.G.M. (P18.01) and K.M.D (P21.10), and at Mistaken Point Ecological Reserve UNESCO World Heritage Site under permits issued to E.G.M. (2016-19; 2021-2024).

Table 1 .
Range of distances (r) for which significantly more or fewer three-point alignments (conga lines) are found than predicted.
We have not used AI-assisted technologies in creating this article.