Josep Mercadal1*†, Susan Tremblay2, Leonie Kraska3, Martin A. Hutten4 and Pau Formosa-Jordan15*
- 1Department of Plant Developmental Biology, Max Planck Institute for Plant Breeding Research, Cologne, Germany
- 2New York State Museum, Madison Avenue, Albany, NY, 12230, USA
- 3LadHyX, CNRS, École Polytechnique, Institut Polytechnique de Paris, Palaiseau Cedex 91128, France
- 4DOI National Park Service, Glacier Bay National Park & Preserve, Gustavus, AK, USA
- 5Cluster of Excellence on Plant Sciences (CEPLAS), Max Planck Institute for Plant Breeding Research, Cologne, Germany
- *Corresponding authors. Email: jmercadal{at}mpipz.mpg.de (J.M.), pformosa{at}mpipz.mpg.de (P.F.-J.)
bioRxiv preprint DOI: https://doi.org/10.64898/2025.12.27.696693
Posted: December 27, 2025, Version 1
Copyright: † Lead author
Abstract
The incompleteness of the fossil record and the absence of genetic information limit our ability to understand ancient developmental processes and reconstruct their evolution. Here, we use comparative and mathematical analyses to uncover an evolutionarily conserved patterning mechanism underlying liverwort oil body cell specification. We demonstrate that patterns of oil body cell homologs in fossils of the oldest known liverwort (Metzgeriothallus sharonae) exhibit the spatial organization characteristic of lateral inhibition. A mechanism with local activation and lateral inhibition not only reproduces the patterns observed in fossils, but also oil body cell patterns in the extant liverworts Treubia lacunosa, Apotreubia nana, and Marchantia polymorpha. Our results point to an ancestral patterning mechanism that has been modulated by lineage-specific selective pressures over evolutionary time.
Evolutionary developmental biology has historically relied on comparative anatomy, morphometrics, and systematics, with fossil evidence playing a crucial role in inferring developmental patterns and their evolutionary relationships1–3. While modern approaches to understanding the evolution of development have gravitated towards molecular techniques, these remain inadequate when confronted with ancient fossil data, where genetic material is rarely preserved. In particular, reconstructing dynamical aspects of development—such as morphogenesis and cell-specification mechanisms—remains challenging and often underexplored. Recent advances in imaging and mathematical modeling have revived comparative approaches4–8, making it possible to recognize morphological fingerprints of developmental mechanisms in living and extinct organisms9–13. In animals, these include the evolution of teeth development14–17, the fin-to-limb transition18–20, and patterning of epidermal structures21–23, among others. While less studied in this context, plants offer a unique opportunity to investigate ancestral developmental mechanisms, as preserved cell walls and fixed cell positions facilitate the recovery of developmentally informative anatomy24–28. Evidence includes the presence of ancient polar auxin transport from growth patterns in wood and vascular tissues29,30, or polyploidy and conserved patterning in fossil stomata31–33.
A striking example of well-preserved cellular patterning can be observed in the Middle Devonian (c. 388 Ma) fossil Metzgeriothallus sharonae34, the oldest known liverwort. A defining feature of these fossils is the presence of dark, scattered cells that closely resemble oil body (OB) cells of modern liverworts (Fig. 1A), suggesting evolutionary homology35,36 and implicating these structures in the chemical defense of liverworts against arthropod herbivory37. Dark cells are distributed more or less evenly throughout the thallus wings, often forming characteristic salt-and-pepper patterns (Fig. 1). This spatial arrangement is associated with the phenomenon of lateral inhibition, where adjacent cells compete to differentiate into distinct types by inhibiting the same fate as their neighbors, generating regularly spaced patterns38–40.

Figure 1:Light micrographs of the Middle Devonian fossil liverwort Metzgeriothalus sharonae.
A. Thallus comprised of midrib and unistratose wings with scattered dark cells. Scale bar = 1 mm. B–E. Representative sections of M. sharonae showing sections of different thalli, where scattered dark cells can be readily discerned. These dark cells are regularly distributed in space, mostly isolated (idioblasts) or forming small clusters of two to three cells (B, C). Occasionally, larger, mostly linear clusters of dark cells can also be observed (D). In a few sections, dark cells appear only occasionally and in idioblastic form (E). Scale bars = 0.1 mm.
Here, we show that dark cell patterns in M. sharonae exhibit the characteristic hallmarks of lateral inhibition, including an overrepresentation of isolated dark cells (idioblasts) and a distinctive spatial regularity. We demonstrate that the spatial statistics of dark cells are consistent with cell specification mediated by local activation and lateral inhibition, accounting for the full diversity of patterns observed in fossil samples. Finally, we show that the same cell specification mechanism can reproduce the spatial organization of OB cells in the extant liverworts Treubia lacunosa, Apotreubia nana, and Marchantia polymorpha, supporting the evolutionary conservation of an ancestral patterning process. Our framework explains lineage-specific phenotypic transitions—such as the complete loss of OBs or the shift from ubiquitous to idioblastic OB distribution—as outcomes of a conserved patterning module being repurposed over evolutionary time.
Dark cell pattern diversity and spatial statistics
We collected M. sharonae material from the Cairo quarry in the Catskill Delta deposit of eastern New York state, well-known for its diverse macroflora, and the only known liverwort-dominated community from the Paleozoic34,36,37. From a total of 120 collected M. sharonae samples, we selected 58 sections based on the size (> 1 mm2), good preservation of cellular structures, and clearly defined dark cells. The total number of cells in each section ranged from ≈ 400 to 900. Due to the variability in frequency and regularity of dark cell patterns, we classified sections into four categories: I) regular patterns with high dark cell frequency (ρI ≥ 0.19), where dark cells form salt-and-pepper patterns with a separation of 1 to 2 epidermal cells (Fig. 1B); II) regular patterns with low dark cell frequency (ρII < 0.19), where dark cells are separated by 2 to 5 epidermal cells (Fig. 1C); III) mixed patterns, where the regularity of dark cell patterns is not obvious due to the presence of gaps and large clusters (Fig. 1D); and IV) patterns with very few scattered dark cells (Fig. 1E). Most sections belong to category II (nII = 43), whereas only a few belong to categories I (nI = 7), III (nIII = 5), and IV (nIV = 3). Because dark cells are rare in type IV sections, we excluded them from the statistical analyses, leaving a total of n = 55 fossil sections.
We used the following metrics to characterize the spatial distribution of dark cells (Supplementary Materials; Figs. S1 and S2): dark cell frequency ρ, i.e., the proportion of dark cells relative to the total number of cells in the tissue; the distribution of sizes of dark cell clusters fk, where k represents the cluster size; the distribution of cell distances between pairs of neighboring dark cells Dj, where j stands for the number of separating epidermal cells; and the spatial autocorrelation via Moran’s I statistic. We compared these metrics with those obtained with a random model of cell specification, in which each cell becomes dark with probability p or epidermal with probability 1 − p (Fig. 2; Fig. S3). We performed simulations of the random model on irregular cell lattices with cell numbers similar to those in fossil samples (Fig. 2A). In fossils, dark cells appear with a frequency ρ = 0.16 ± 0.2 (Fig. 2A,B). The distribution of cluster sizes shows a strong excess of idioblasts (clusters of size k = 1), significantly more than what would be expected by random specification (Fig. 2D). In contrast, larger clusters occur at considerably lower proportions than expected by chance (Fig. 2D), and typically take the form of linear arrays, which are distinctive hallmarks of lateral inhibition41,42. Most dark cells are separated by 1 to 3 cells, with an average of ≈ 2.7 cells (Fig. 2E). This regularity contrasts with the random model, where dark cell separation is significantly more variable (Fig. 2E; Fig. S2).

Figure 2:Non-randomness of dark cell patterns in M. sharonae fossils.
A. Representative section of M. sharonae showing a dark cell pattern (top) and a simulated random pattern with similar dark cell frequency (bottom). B. Boxplots showing the frequency of dark cells for fossil sections (brown, n = 55) and for random patterns (gray, n = 55) with dark cell probability p = 0.16 × (1 + 0.7𝒰 [−0.5, 0.5]), where 𝒰 [−0.5, 0.5] denotes a uniform distribution of random numbers between −0.5 and 0.5. With this choice, random patterns exhibit a frequency distribution similar to that of fossil samples. C. Metrics used to determine the spatial statistics of dark cell patterns, including the cluster sizes in different colors and the shortest paths between neighboring dark cells in red (see Fig. S1 and Supplementary Materials for a complete description of these metrics). D. Distribution of cluster sizes in fossil patterns (brown) compared to random patterns (gray). The inset shows the same data on a log-log scale. Nc stands for the total number of clusters analyzed. E. Distribution of distances (shortest paths) between dark cells in fossils (brown) and random patterns (gray). F. Morphospace showing the proportion of idioblastic dark cells versus Moran’s I in fossil (brown) and random (gray) patterns. G. Morphospace showing the average distances between dark cells versus their corresponding coefficient of variation of fossil (brown) and random (gray) patterns. In panels B, F, and G, box widths represent interquartile ranges (IQR = Q3-Q1), horizontal lines represent medians, and whiskers extend from Q1-1.5IQR to Q3+1.5IQR. Error bars and shaded areas in D and E denote ± standard deviations. The shaded regions in panels F and G represent the convex hulls of each dataset. Statistical significance was assessed using the Mann-Whitney U-test; *** p < 0.01; ns: non-significant.
To better understand how dark cell patterns compare to a random model of cell specification, we located each fossil sample in the morphospaces of Moran’s I statistic versus idioblast proportion, and average dark cell distance ⟨λ⟩ versus CVλ —the coefficient of variation of these distances (see Supplementary Materials for how these metrics are defined). Our results show that in these spaces, dark cell patterns cluster without overlapping with random patterns (Fig. 2F,G). Dark cell patterns have a high proportion of idioblasts and negative values of Moran’s I, revealing a spatial organization that the random model cannot reproduce (Fig. 2F). While the average distance between dark cells is similar between the fossils and the random model, the CVλ is significantly lower in the fossils (Fig. 2G), reflecting lower variability in these distances and a characteristic pattern regularity. Altogether, these results show that the statistical properties of dark cell patterns in fossils are incompatible with random cell specification. This suggests that the mechanism of dark cell specification in M. sharonae was not random—as would be expected by a cell-autonomous process—but likely involved cell-to-cell interactions via lateral inhibition, preventing adjacent dark cells from adopting the same fate and generating the observed spatial regularity.
A model with local activation and lateral inhibition reproduces the diversity of dark cell patterns
Developmental mechanisms generating cell-type-specific patterning typically rely on cell-to-cell communication through variants of lateral inhibition43–45. While the genetic and molecular factors involved in patterning might differ significantly across biological systems, the dynamics and spatial statistics are often telltale signs of the underlying developmental processes and offer crucial information for distinguishing between mechanisms46,47. Due to the lack of genetic and dynamical information, we compared the spatial statistics of dark cell patterns in M. sharonae with those obtained with a minimal cellular model of pattern formation, where a local activator u induces the activity of a mobile inhibitor v42,48 (Fig. 3A; Supplementary Materials). At each time point, the state of each cell is described by the values of u(t) and v(t). We assume the activator is constantly produced, self-activates, and is degraded linearly and by the inhibitor (Fig. 3B). The inhibitor is also constantly produced, is induced by the activator, is linearly degraded, and can diffuse to adjacent cells (Fig. 3B). We explored the steady-state patterning properties of the model in an irregular cellular tissue by varying the parameters γ—the maximum rate at which u activates v—and α0—the basal production rate of u (Fig. 3C; Fig. S4). We use the stationary values of u as a proxy for cell fate, with high u indicating dark cells and low u indicating epidermal cells.

Figure 3:A model with local activation and lateral inhibition reproduces the diversity of dark cell patterns in M. sharonae fossils.
A. Sketch of the interactions in a cellular model of local activation and lateral inhibition. The state of two adjacent undifferentiated cells can be destabilized due to diffusive lateral inhibition, priming them to adopt opposite fates. B. Equations of the mathematical model corresponding to the interactions depicted in A. The key regulations are highlighted in blue. C. Stability diagram showing different patterning regimes obtained with the mathematical model as a function of the parameters γ and α0 (Supplementary Materials; Fig. S4). The red square denotes a cusp point where two saddle-node bifurcations collide and annihilate; the green square denotes the collision between a subcritical Turing bifurcation and two saddle-node bifurcations. Snapshots of representative simulated patterns are displayed across the stability diagram. D. Comparison of the spatial statistics of high-density fossil samples (brown; n = 7) with the mechanistic (pink; n = 55) and random (gray; n = 100) models (see Fig. S5). E. Comparison of low-density fossil samples (brown; n = 43) with the mechanistic (pink; n = 826) and random (gray; n = 100) models (see Fig. S5). F. In a few sections, dark cells are inhomogeneously distributed (left) or appear only occasionally (right). Statistics were not performed due to the small number of sections of these types. Nevertheless, the mathematical model can qualitatively reproduce the corresponding patterns (Fig. S5). G. Several patterns predicted by the model are not observed in the fossils. Cell colors represent the values of the variable u in each cell, ranging from u ∈ [0.2, 1.2] (see Fig. S4). The values of α0 and γ used for panels D–G are specified in Fig. S5. The rest of the parameters are detailed in Table S1.
Analyzing the model through bifurcation theory and linear stability analysis38,49 reveals different patterning regimes (Fig. 3C; Fig. S4). If α0 is small, the activator is not strong enough to induce cell differentiation, and no dark cells emerge (Fig. 3C; region 1). Conversely, for large α0, all cells become dark (region 2). Patterns with different dark cell frequencies and regularities can emerge for intermediate values of (γ, α0). The system has a bistable region where cellular fates are strongly influenced by initial conditions (region 3). Here, lateral inhibition does not play a significant role because all possible cell combinations are stable (Fig. S4). The model can also generate localized patterns with different densities (regions 4 and 5) that are triggered only by specific initial conditions (Fig. S4). Finally, for a wide range of (γ, α0), the system can generate self-organized patterns (region 6), whereby ordered spatial patterns emerge spontaneously due to the amplification of inhomogeneous perturbations. Patterns across this region differ in their spatial statistics, including dark cell frequency, idioblast proportion, and average dark cell distance (Fig. S4). Mixed patterns—where periodic domains coexist with more random-looking patterns—can be obtained when the values of (γ, α0) are close to the bistable region (regions 4,5,6 close to region 3; Fig. S4).
We performed a parameter sweep over the pairs (γ, α0) to test whether the model could reproduce the statistical metrics of fossil patterns (dark cell frequency, idioblast proportion, Moran’s I, average dark cell distance ⟨λ⟩, and CVλ) across the four categories of our classification (Fig. 3D,F; Fig. S5). For each category, we identified the pairs (γ, α0) that generated patterns in which all metrics simultaneously fell within their observed experimental ranges (Supplementary Materials; Fig. S5). We found that the model could reproduce all the observed phenotypes, generating patterns in morphospace regions that virtually overlap or contain the experimental results (Fig. 3D,E). Notably, the model yields multiple parameter pairs that satisfy all these conditions (Fig. S5), underscoring the mechanisms’ ability to generate similar phenotypes through different routes. Uncommon patterns, such as those with no well-defined regularity or very few localized dark cells, could also be reproduced by the model (Fig. 3F; Fig. S5). Additionally, the model can generate phenotypes that are not observed in fossil samples. These include tissues completely covered in dark cells, tissues with mostly dark cells and few localized light cells, labyrinthine patterns, and tissues without dark cells (Fig. 3G; Fig. S4,S5).
Taken together, these results demonstrate that a cell specification mechanism through local activation and lateral inhibition can explain the full variation of dark cell patterns in M. sharonae. This flexibility in generating diverse patterning phenotypes highlights the evolutionary potential of the mechanism, enabling adaptive responses to varying selective pressures.
An evolutionarily conserved patterning mechanism for oil body cell specification
Our results suggest a cell patterning mechanism underlying the specification of dark cells in M. sharonae, the oldest known fossil liverwort. Dark cells in M. sharonae and other Paleozoic liverworts are considered homologous to the OB cells of modern liverworts36,37. Hence, to investigate whether this mechanism has been conserved across the liverwort lineage, we analyzed the spatial distribution of OB cells in thalli of the extant liverworts Treubia lacunosa, Apotreubia nana, and Marchantia polymorpha (Fig. 4). In T. lacunosa and A. nana, cells harboring OBs are on average larger than non-OB cells, and most appear as idioblasts (Fig. 4A,C). Notably, they occur at a similar frequency as in M. sharonae (ρ ≈ 0.15; Fig. 4B,D). Most of the OB cell patterns in T. lacunosa and A. nana are inconsistent with random cell specification, despite some overlap in the regions occupied by the morphospaces (Fig. 4B,D). Instead, both the idioblast proportion, Moran’s I, average OB cell distance ⟨λ⟩, and its coefficient of variation CVλ are consistent with patterns shaped by lateral inhibition. To test whether our model could reproduce the observed statistics, we looked for pairs (γ, α0) whose statistical metrics reproduced OB cell patterns in T. lacunosa and A. nana. We found several pairs that fulfilled all the conditions (Fig. 4B,D, pink; Fig. S6), spanning a morphospace that virtually overlapped or contained the experimental results. Therefore, OB cell patterns in T. lacunosa and A. nana are consistent with cell specification governed by local activation and lateral inhibition.

Figure 4:Lateral inhibition patterning of oil body cells is conserved across liverwort taxa.
A. Distribution of OB cells in T. lacunosa. B. Spatial statistics comparing experimental OB cell patterns in T. lacunosa (brown; n = 50) with patterns obtained with the mechanistic (pink; n = 777) and random (gray; n = 55) models. C. Distribution of OB cells in A. nana. D. Spatial statistics comparing OB cell patterns in A. nana (green; n = 16) with patterns obtained with the mechanistic (pink; n = 305) and random (gray; n = 55) models. E. M. polymorpha thallus stained with BODIPY, exposing oil bodies and cell boundaries ( Fig. S7). F. Spatial statistics comparing OB cell patterns in M. polymorpha (gold; n = 21) with patterns obtained with the mechanistic (pink; n = 84) and random (gray; n = 55) models. Statistical differences in panels B, D, and F were tested using ANOVA followed by Tukey HSD (p < 0.01); letters indicate statistically significant groups. Scale bars in A, C, and E: 0.1 mm.
Finally, we tested whether our model could reproduce OB cell patterns in M. polymorpha. In M. polymorpha thalli, OB cell frequency is lower than in fossils, T. lacunosa and A. nana, and almost exclusively in idioblastic form (Fig. 4E,F; Fig. S8). Due to this lower frequency, the average distance between OB cells is higher than in the other species. By contrast, the coefficient of variation spans the same scale, suggesting that despite differences in frequency and OB cell distance, the pattern variability is conserved. Comparing these distributions with what a random model with the same OB cell frequency reveals that patterns in M. polymorpha are also shaped by lateral inhibition. Indeed, our model reproduced all the spatial statistics associated with these patterns (Fig. 4F; Fig. S6).
Taken together, our results show that a mechanism with local activation and lateral inhibition can reproduce the spatial patterns of OB cells not only in M. sharonae fossils, but also in the extant liverworts T. lacunosa, A. nana, and M. polymorpha. This consistency across distantly related liverwort taxa points to a deeply conserved developmental mechanism underlying liverwort OB cell specification.
Discussion
Oil bodies play important roles in lipid and secondary metabolite storage, thereby contributing significantly to plant defense mechanisms50–54. The regular and idioblastic nature of OB cells might reflect a strategy to minimize unnecessary overlap, potentially optimizing their protective or signaling functions. Despite their putative defensive role, the frequency and spatial distribution of OB vary among liverwort taxa. In Haplomitriopsida, taxa in the genus Haplomitrium have several OBs in all cells. In contrast, those in the Treubiales (like T. lacunosa and A. nana) have a single OB in idioblastic cells52,55 (Fig. 5A). In Marchantiopsida, most taxa have a single OB in idioblastic cells, while in Jungermanniopsida (where M. sharonae is currently assigned34,36), most taxa have OBs in all cells52,55 (Fig. 5A). The three liverwort classes are estimated to have diverged from a common ancestor no later than the Late Silurian, and last shared a common ancestor with mosses in the Late Ordovician56,57. These divergence times highlight the ancient existence of OBs as a crucial adaptation to life on land.

Figure 5:Oil body cell evolution through the lens of a flexible patterning module.
A. Oil body cell distribution across the liverwort phylogeny. Branch colors indicate the general spatial arrangement of OB cells in the corresponding lineages: Haplomitriopsida (green), Marchantiopsida (magenta), and Jungermanniopsida (blue), together with a description of how OB cells are distributed. MRCA: most recent common ancestor. B. Patterning regimes (Fig. 3C) showing where OBs appear in all cells (light blue), in idioblastic cells (magenta), or are absent (yellow). Phenotypic transitions during evolution can be interpreted as excursions in this space. For instance, the transition from having OBs in all cells to their disappearance (as seen in Nowellia or Anthelia) can be interpreted either as a smooth change (e.g, through slowly changing environments) or as a sharp transition (e.g., through strong mutations). Similarly, the transition from OBs in idioblasts to OB loss (as occurs in many genera of Marchantiopsida) can be interpreted as the crossing of a bifurcation boundary.
The variation in type and spatial organization of OBs suggests a developmental flexibility that can be co-opted to produce more specialized adaptations suited to changing ecological contexts. Notably, most of the phenotypic variation in OB cell patterns was already present in M. sharonae, suggesting that developmental plasticity could have supported timely responses to the changing environments and emerging biotic stresses of Devonian landscapes. In our modeling framework, phenotypic variation reflects the mechanisms’ ability to explore different regions of the patterning space (Figs. 3C and 5B). Modulating the relative effect of the activator and inhibitor constitutes a route for evolutionary change by shifting the system between different stable states4,7. These adjustments can lead to bifurcations in the developmental dynamics, where small variations in the parameters result in large phenotypic changes. In this regard, phenotypic transitions such as the loss of OBs or the shift from ubiquitous to idioblastic OB distribution reflect incursions into qualitatively distinct patterning domains (Fig. 5B), which can occur as a result of, e.g., genetic mutations, biotic stresses, or environmental perturbations. We thus hypothesize that the patterning mechanism underlying OB cell specification arose early in liverwort evolution and has persisted despite extensive diversification, being repurposed over time due to lineage-specific selective pressures.
Despite being one of the defining synapomorphies of the group, OBs have rarely been used to explore evolutionary relationships within the liverwort clade, except for a few studies35,58. This is striking considering that the internal and spatial variability of OB cells is taxonomically informative and could offer insights into the adaptive strategies of early land plants. Idioblastic dark cells, very similar in structure and distribution to those of M. sharonae, have also been described in fossils from the early liverworts Pallaviciniites devonicus36,59 and Treubiites kidstonii60,61, and in Paleozoic mosses62–64, suggesting OBs may have originated much earlier than currently assumed and could even predate the divergence of bryophyte lineages.
Recent studies in M. polymorpha are shedding light on the molecular and developmental mechanisms underlying OB cell specification. The transcription factors MpC1HDZ, MpERF13, and MYB02 act as activators of OB cell specification, as loss-of-function mutants lack OB cells altogether53,65–67. Notably, MpC1HDZ and MpERF13 activate the expression of MYB0267. The transcription factor MpTGA represses the expression of both MpC1HDZ and MpERF13, and loss-of-function mutants of MpTGA show a dramatic increase in OB cells. Our model can qualitatively reproduce all these phenotypes (Fig. S9). In addition to modulating the spatial distribution of OB cells, several of these genes are necessary for the synthesis of OB-specific compounds. TGA represses MpSYP12B, an OB-resident SNARE protein that localizes in the OB membrane53,68, and upregulates terpenoid biosynthesis enzymes encoded by TPS and MTPSL68. The transporter MpABCG1 was recently shown to mediate the OB accumulation of sesquiterpenes, compounds that contribute to herbivore resistance69. Finally, the immune regulator MpNPR has been shown to control the expression of MpERF13, MpSYP12B, MpMYB02, and MpABCG1, thus linking the defensive function of OBs with the developmental processes underlying OB cell distribution. Despite the lack of any feedback mechanism that could function as an activator-inhibitor system capable of cell-type patterning, these studies are beginning to link the components of what appears to be a complex gene regulatory network (Fig. S10), providing a foundation for future reconstructions of the developmental mechanisms underlying OB formation and evolution.
Supporting information
Supplementary Materials (including Materials and Methods, Supplementary Figures and Supplementary Tables)[supplements/696693_file02.pdf]
Funding
J.M. acknowledges financial support from the Alexander von Humboldt-Stiftung (ESP 1236193 HFST-P) and the European Union’s Horizon 2023 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 101153033. S.T. acknowledges financial support from the University of California, Berkeley’s Department of Integrative Biology, the Erwin Resetko Family scholarship fund, and the Evolving Earth Foundation. L.K. is supported by the Agence Nationale de la Recherche grant ANR-23-CE13-0035-01. J.M. and P.F.-J. also received funding from the DFG through the Cluster of Excellence CEPLAS (EXC 2048/1 Project ID: 390686111) and from a Core Grant from the Max Planck Society.
Author contributions
J.M. designed the research, performed the image analysis, and analyzed the mathematical models. S.T. collected the samples of M. sharonae and T. lacunosa, and generated the images. L.K. produced the confocal images of M. polymorpha. M.H. collected the samples of A. nana and generated the images. J.M. wrote the manuscript with input from the other authors.
Competing interests
The authors declare no competing interests.
Data and materials availability
All data reported in this paper will be shared by the lead contact upon request. All original code for data analysis and simulations is available in https://github.com/josepmercadal/fossils.
ACKNOWLEDGEMENTS
J.M. thanks the members of the Formosa-Jordan group for valuable discussions. S.T. thanks Taehee Han and Victoria Chu for providing laboratory assistance for the selection and preparation of fossils, David Glenny and Bill Malcolm for field support in New Zealand, and Linda Hernick for guidance and advice during field work in New York. L.K. thanks Arezki Boudaoud for hosting and support, and Pierre Mahou for providing access to the microscope facilities and technical support. M. H. thanks Edward Kulack for help in preparing LIDAR data and James Walton for help with fieldwork.
Footnotes
This pre-print is available under a Creative Commons License (Attribution 4.0 International), CC BY 4.0, as described at http://creativecommons.org/licenses/by/4.0/
References
- [1]. B. K. Hall, “Palaeontology and evolutionary developmental biology: a science of the nineteenth and twenty–first centuries,” Palaeontology, vol. 45, no. 4, pp. 647–669, 2002.
- [2]. S. F. Gilbert, “The morphogenesis of evolutionary developmental biology,” International Journal of Developmental Biology, vol. 47, no. 7-8, p. 467, 2003.
- [3]. R. Amundson, The changing role of the embryo in evolutionary thought: roots of evo-devo. Cambridge University Press, 2005.
- [4]. P. Alberch, “From genes to phenotype: dynamical systems and evolvability.,” Genetica, vol. 84, no. 1, pp. 5–11, 1991.
- [5]. D. Sepkoski, “Stephen Jay Gould, Jack Sepkoski, and the ‘Quantitative Revolution’ in American Paleobiology,” Journal of the History of Biology, vol. 38, pp. 209–237, 2005.
- [6]. J. A. Cunningham, I. A. Rahman, S. Lautenschlager, E. J. Rayfield, and P. C. Donoghue, “A Virtual World of Paleontology,” Trends in Ecology & Evolution, vol. 29, no. 6, pp. 347–357, 2014.
- [7]. J. Jaeger and N. Monk, “Bioattractors: Dynamical Systems Theory and the Evolution of Regulatory Processes,” The Journal of Physiology, vol. 592, no. 11, pp. 2267–2281, 2014.
- [8]. J. DiFrisco and J. Jaeger, “Homology of Process: Developmental Dynamics in Comparative Biology,” Interface Focus, vol. 11, no. 3, p. 20210007, 2021.
- [9]. H. Sanders, G. W. Rothwell, and S. Wyatt, “Paleontological context for the developmental mechanisms of evolution,” International Journal of Plant Sciences, vol. 168, no. 5, pp. 719–728, 2007.
- [10]. A. D. Chipman and G. D. Edgecombe, “Developing an integrated understanding of the evolution of arthropod segmentation using fossils and evo-devo,” Proceedings of the Royal Society B, vol. 286, no. 1912, p. 20191881, 2019.
- [11]. G. W. Rothwell and A. M. Tomescu, “Structural Fingerprints of Development at the Intersection of Evo-Devo and the Fossil Record,” in Evolutionary Developmental Biology: a Reference Guide, pp. 573–602, Springer, 2021.
- [12]. C. Haug and J. T. Haug, “Methods and practices in paleo-evo-devo,” Evolutionary Developmental Biology: A Reference Guide, pp. 1151–1164, 2021.
- [13]. A. M. Tomescu and G. W. Rothwell, “Fossils and Plant Evolution: Structural Fingerprints and Modularity in the Evo-Devo Paradigm,” EvoDevo, vol. 13, no. 1, p. 8, 2022.
- [14]. I. Salazar-Ciudad and J. Jernvall, “A computational model of teeth and the developmental origins of morphological variation,” Nature, vol. 464, no. 7288, pp. 583–586, 2010.
- [15]. E. Harjunmaa, K. Seidel, T. Häkkinen, E. Renvoisé, I. J. Corfe, A. Kallonen, Z.-Q. Zhang, A. R. Evans, M. L. Mikkola, Salazar-Ciudad, et al., “Replaying evolutionary transitions from the dental fossil record,” Nature, vol. 512, no. 7512, pp. 44–48, 2014.
- [16]. A. Ortiz, S. E. Bailey, G. T. Schwartz, J.-J. Hublin, and M. M. Skinner, “Evo-devo models of tooth development and the origin of hominoid molar diversity,” Science advances, vol. 4, no. 4, p. eaar2334, 2018.
- [17]. R. Zimm, F. Berio, M. Debiais-Thibaud, and N. Goudemand, “A shark-inspired general model of tooth morphogenesis unveils developmental asymmetries in phenotype transitions,” Proceedings of the National Academy of Sciences, vol. 120, no. 15, p. e2216959120, 2023.
- [18]. J. Raspopovic, L. Marcon, L. Russo, and J. Sharpe, “Digit patterning is controlled by a bmp-sox9-wnt turing network modulated by morphogen gradients,” Science, vol. 345, no. 6196, pp. 566–570, 2014.
- [19]. K. Onimaru, L. Marcon, M. Musy, M. Tanaka, and J. Sharpe, “The fin-to-limb transition as the re-organization of a turing pattern,” Nature communications, vol. 7, no. 1, p. 11582, 2016.
- [20]. T. A. Stewart, J. B. Lemberg, N. K. Taft, I. Yoo, E. B. Daeschler, and N. H. Shubin, “Fin ray patterns at the fin-to-limb transition,” Proceedings of the National Academy of Sciences, vol. 117, no. 3, pp. 1612–1620, 2020.
- [21]. R. L. Cooper, A. P. Thiery, A. G. Fletcher, D. J. Delbarre, L. J. Rasch, and G. J. Fraser, “An ancient turing-like patterning mechanism regulates skin denticle development in sharks,” Science advances, vol. 4, no. 11, p. eaau5484, 2018.
- [22]. M. Staps, P. W. Miller, C. E. Tarnita, and R. Mallarino, “Development shapes the evolutionary diversification of rodent stripe patterns,” Proceedings of the National Academy of Sciences, vol. 120, no. 45, p. e2312077120, 2023.
- [23]. C.-C. Tseng, T. E. Woolley, T.-X. Jiang, P. Wu, P. K. Maini, R. B. Widelitz, and C.-M. Chuong, “Gap junctions in turing-type periodic feather pattern formation,” PLoS Biology, vol. 22, no. 5, p. e3002636, 2024.
- [24]. C. K. Boyce, “The Fossil Record of Plant Physiology and Development—What Leaves Can Tell Us,” The Paleontological Society Papers, vol. 14, pp. 133–146, 2008.
- [25]. P.-M. Delaux, A. J. Hetherington, Y. Coudert, C. Delwiche, C. Dunand, S. Gould, P. Kenrick, F.-W. Li, H. Philippe, S. A. Rensing, et al., “Reconstructing Trait Evolution in Plant Evo–Devo Studies,” Current Biology, vol. 29, no. 21, pp. R1110–R1118, 2019.
- [26]. E. Petrone-Mendoza, F. Vergara-Silva, and M. E. Olson, “Plant morpho evo-devo,” Trends in Plant Science, 2023.
- [27]. A. J. Hetherington, J. G. Dubrovsky, and L. Dolan, “Unique Cellular Organization in the Oldest Root Meristem,” Current Biology, vol. 26, no. 12, pp. 1629–1633, 2016.
- [28]. A. J. Hetherington, “The Role of Fossils for Reconstructing the Evolution of Plant Development,” Development, vol. 151, no. 20, 2024.
- [29]. G. Rothwell and S. Lev-Yadun, “Evidence of polar auxin flow in 375 million-year-old fossil wood,” American Journal of Botany, vol. 92, no. 6, pp. 903–906, 2005.
- [30]. G. W. Rothwell, H. Sanders, S. E. Wyatt, and S. Lev-Yadun, “A fossil record for growth regulation: The role of auxin in wood evolution1,” Annals of the Missouri Botanical Garden, vol. 95, no. 1, pp. 121–134, 2008.
- [31]. J. Masterson, “Stomatal Size in Fossil Plants: Evidence for Polyploidy in Majority of Angiosperms,” Science, vol. 264, no. 5157, pp. 421–424, 1994.
- [32]. J. C. McElwain and M. Steinthorsdottir, “Paleoecology, Ploidy, Paleoatmospheric Composition, and Developmental Biology: a Review of the Multiple Uses of Fossil Stomata,” Plant Physiology, vol. 174, no. 2, pp. 650–664, 2017.
- [33]. P. J. Rudall and R. M. Bateman, “Leaf Surface Development and the Plant Fossil Record: Stomatal Patterning in Bennettitales,” Biological Reviews, vol. 94, no. 3, pp. 1179–1194, 2019.
- [34]. L. V. Hernick, E. Landing, and K. E. Bartowski, “Earth’s oldest liverworts—metzgeriothallus sharonae sp. nov. from the middle devonian (givetian) of eastern new york, usa,” Review of Palaeobotany and Palynology, vol. 148, no. 2-4, pp. 154–162, 2008.
- [35]. R. M. Schuster, “The hepaticae and anthocerotae of north america,” Vols 1–6: New York, NY, USA and London, UK: Columbia University Press; Vols 5 & 6: Chicago, IL, USA: Field Museum of Natural History., 1966.
- [36]. S. Tremblay and J. Mercadal, “Homology of the dark cells of Paleozoic liverworts with the specialized oil body cells of modern liverworts (Marchantiophyta),” bioRxiv, no. doi:10.64898/2025.12.24.696392, 2025.
- [37]. C. C. Labandeira, S. L. Tremblay, K. E. Bartowski, and L. VanAller Hernick, “Middle devonian liverwort herbivory and antiherbivore defence,” New Phytologist, vol. 202, no. 1, pp. 247–258, 2014.
- [38]. J. R. Collier, N. A. Monk, P. K. Maini, and J. H. Lewis, “Pattern formation by lateral inhibition with feedback: a mathematical model of delta-notch intercellular signalling,” Journal of theoretical Biology, vol. 183, no. 4, pp. 429–446, 1996.
- [39]. M. Cohen, M. Georgiou, N. L. Stevenson, M. Miodownik, and B. Baum, “Dynamic filopodia transmit intermittent delta-notch signaling to drive pattern refinement during lateral inhibition,” Developmental cell, vol. 19, no. 1, pp. 78–89, 2010.
- [40]. N. Guisoni, R. Martinez-Corral, J. Garcia-Ojalvo, and J. de Navascués, “Diversity of fate outcomes in cell pairs under lateral inhibition,” Development, vol. 144, no. 7, pp. 1177–1186, 2017.
- [41]. A. Thamm, T. E. Saunders, and L. Dolan, “MpFEW RHIZOIDS1 miRNA-Mediated Lateral Inhibition Controls Rhizoid Cell Patterning in Marchantia polymorpha,” Current Biology, vol. 30, pp. 1905–1915, 2020.
- [42]. J. Mercadal, M. Ferreira-Guerra, A. I. Caño-Delgado, and M. Ibañes, “Cell-Type Specification at Criticality Underlies Rhizoid Patterning in the Liverwort Marchantia polymorpha,” Cell Reports, vol. 45, no. 1, p. 116728, 2026.
- [43]. I. Salazar-Ciudad, J. Jernvall, and S. A. Newman, “Mechanisms of Pattern Formation in Development and Evolution,” Development, vol. 130, no. 10, pp. 2027–37, 2003.
- [44]. K. U. Torii, “Two-dimensional Spatial Patterning in Developmental Systems,” Trends in Cell Biology, vol. 22, no. 8, pp. 438–446, 2012.
- [45]. F. Schweisguth and F. Corson, “Self-Organization in Pattern Formation,” Developmental Cell, vol. 49, no. 5, pp. 659–677, 2019.
- [46]. T. W. Hiscock and S. G. Megason, “Mathematically guided approaches to distinguish models of periodic patterning,” Development, vol. 142, no. 3, pp. 409–419, 2015.
- [47]. H. D. Summers, J. W. Wills, and P. Rees, “Spatial Statistics is a Comprehensive Tool for Quantifying Cell Neighbor Relationships and Biological Processes via Tissue Image Analysis,” Cell Reports Methods, vol. 2, no. 11, 2022.
- [48]. H. Meinhardt and A. Gierer, “Pattern formation by local self-activation and lateral inhibition,” BioEssays, vol. 22, pp. 753–760, 2000.
- [49]. S. H. Strogatz, Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering. CRC press, 2018.
- [50]. X. He, Y. Sun, and R.-L. Zhu, “The Oil Bodies of Liverworts: Unique and Important Organelles in Land Plants,” Critical Reviews in Plant Sciences, vol. 32, no. 5, pp. 293–302, 2013.
- [51]. M. Tanaka, T. Esaki, H. Kenmoku, T. Koeduka, Y. Kiyoyama, T. Masujima, Y. Asakawa, and K. Matsui, “Direct evidence of specific localization of sesquiterpenes and marchantin a in oil body cells of marchantia polymorpha l.,” Phytochemistry, vol. 130, pp. 77–84, 2016.
- [52]. F. Romani, J. R. Flores, J. I. Tolopka, G. Suárez, X. He, and J. E. Moreno, “Liverwort oil bodies: diversity, biochemistry, and molecular cell biology of the earliest secretory structure of land plants,” Journal of Experimental Botany, vol. 73, no. 13, pp. 4427–4439, 2022.
- [53]. T. Kanazawa, H. Morinaka, K. Ebine, T. L. Shimada, S. Ishida, N. Minamino, K. Yamaguchi, S. Shigenobu, T. Kohchi, A. Nakano, and T. Ueda, “The liverwort oil body is formed by redirection of the secretory pathway,” Nature Communications, vol. 11, no. 1, p. 6152, 2020.
- [54]. F. Schweizer, I. Monte, R. Solano, and P. Reymond, “Marchantia polymorpha defense against snail herbivory,” Plant-Environment Interactions, vol. 6, no. 2, p. e70052, 2025.
- [55]. J. R. Flores, S. A. Catalano, J. Muñoz, and G. M. Suárez, “Combined phylogenetic analysis of the subclass marchantiidae (marchantiophyta): towards a robustly diagnosed classification,” Cladistics, vol. 34, no. 5, pp. 517–541, 2018.
- [56]. J. Heinrichs, J. Hentschel, R. Wilson, K. Feldberg, and H. Schneider, “Evolution of leafy liverworts (jungermanniidae, marchantiophyta): estimating divergence times from chloroplast dna sequences using penalized likelihood with integrated fossil evidence,” Taxon, vol. 56, no. 1, pp. 31–44, 2007.
- [57]. J. Bechteler, G. Peñaloza-Bojacá, D. Bell, J. Gordon Burleigh, S. F. McDaniel, E. Christine Davis, E. B. Sessa, A. Bippus, D. Christine Cargill, S. Chantanoarrapint, et al., “Comprehensive phylogenomic time tree of bryophytes reveals deep relationships and uncovers gene incongruences in the last 500 million years of diversification,” American Journal of Botany, vol. 110, no. 11, p. e16249, 2023.
- [58]. B. Crandall-Stotler, R. Stotler, and D. Long, “Phylogeny and classification of the marchantiophyta,” Edinburgh journal of botany, vol. 66, no. 1, pp. 155–198, 2009.
- [59]. F. M. Hueber, “Hepaticites devonicus: A New Fossil Liverwort from the Devonian of New York,” Annals of the Missouri Botanical Garden, pp. 125–131, 1961.
- [60]. J. Walton, “Carboniferous bryophyta: I. hepaticae,” Annals of Botany, vol. 39, no. 155, pp. 563–572, 1925.
- [61]. V. Krassilov and R. Schuster, “Paleozoic and mesozoic fossils,” New manual of bryology, vol. 2, pp. 1172–1193, 1984.
- [62]. E. V. M. Maslova, Y. V. Mosseichik, I. A. Ignatiev, O. V. Ivanov, and M. S. Ignatov, “On the leaf development in palaeozoic mosses of the order protosphagnales,” Arctoa, vol. 21, pp. 241–264, 2012.
- [63]. O. V. Ivanov, E. V. Maslova, and M. S. Ignatov, “Development of the sphagnoid areolation pattern in leaves of palaeozoic protosphagnalean mosses,” Annals of Botany, vol. 122, no. 5, pp. 915–925, 2018.
- [64]. M. S. Ignatov and E. V. Maslova, “Fossil mosses: what do they tell us about moss evolution?,” Bryophyte Diversity and Evolution, vol. 43, no. 1, pp. 72–97, 2021.
- [65]. F. Romani, E. Banić, S. N. Florent, T. Kanazawa, J. Q. Goodger, R. A. Mentink, T. Dierschke, S. Zachgo, T. Ueda, J. L. Bowman, M. Tsiantis, and J. E. Moreno, “Oil body formation in marchantia polymorpha is controlled by mpc1hdz and serves as a defense against arthropod herbivores,” Current Biology, vol. 30, no. 14, pp. 2815–2828, 2020.
- [66]. H. Kubo, S. Nozawa, T. Hiwatashi, Y. Kondou, R. Nakabayashi, T. Mori, K. Saito, K. Takanashi, T. Kohchi, and K. Ishizaki, “Biosynthesis of riccionidins and marchantins is regulated by r2r3-myb transcription factors in marchantia polymorpha,” Journal of plant research, vol. 131, pp. 849–864, 2018.
- [67]. H. Kubo, K. Sunhwa, H. Teramori, and K. Takanashi, “Mp r2r3-myb2 is a key regulator of oil body formation in marchantia polymorpha,” Planta, vol. 260, no. 3, p. 68, 2024.
- [68]. N. Gutsche, J. Koczula, M. Trupp, M. Holtmannspötter, M. Appelfeller, O. Rupp, A. Busch, and S. Zachgo, “Mp tga, together with mp npr, regulates sexual reproduction and independently affects oil body formation in marchantia polymorpha,” New Phytologist, vol. 241, no. 4, pp. 1559–1573, 2024.
- [69]. E. C. Forestier, P. Asprilla, F. Romani, I. Bonter, E. Frangedakis, and J. Haseloff, “The abcg1 transporter facilitates sesquiterpene accumulation in marchantia polymorpha oil bodies,” bioRxiv, pp. 2025–02, 2025.
