Non-technical Summary
Coastal tropical wetlands are some of the earliest recorded ecosystems on land and likely played a crucial role in the diversity of behaviors and adaptations seen in modern terrestrial animals. In the late Carboniferous (ca. 323–299 Ma), the relatives of modern tetrapods (reptiles, mammals, and amphibians) rapidly diversified in terms of the kinds of ecological niches they occupied, with their remains being predominantly found in coastal wetland environments. What is not clear is whether tetrapod communities in these late Carboniferous wetlands approached levels of ecological diversity similar to those of modern coastal tropical wetlands. We investigate differences in ecological diversity between late Carboniferous and modern wetland communities by comparing their body-mass distributions. Body mass is related to a wide range of aspects of an animal’s niche, so patterns in the array of body masses from a community can be useful as a broad metric for estimating the diversity of ecological roles in that ecosystem. Body masses cannot be directly measured in fossil tetrapods, so we expanded upon a preexisting method that relates the dimensions of back vertebrae to body mass in living tetrapods to estimate fossil body masses. We estimated body masses from two late Carboniferous fossil tetrapod communities from coastal wetland deposits in Nova Scotia and compared the arrays of fossil body masses to those of modern coastal tropical wetland communities. We found that, after accounting for differences resulting from fossil preservation, most differences between fossil and modern wetland communities could be explained by the limited dietary diversity of Carboniferous tetrapods.
Introduction
Coastal tropical wetlands and estuaries are vitally important ecosystems and have existed for the majority of the Phanerozoic (Greb et al. Reference Greb, DiMichele, Gastaldo, Greb and DiMichele2006; Ramsar Convention Secretariat 2013). The exact inhabitants have changed dramatically over the course of the Phanerozoic (Greb et al. Reference Greb, DiMichele, Gastaldo, Greb and DiMichele2006), but exactly how different early terrestrial ecosystems were from the modern coastal tropics in terms of functional diversity remains an open question. Though modern terrestrial vertebrates are not as well studied in coastal tropical wetlands as they are in inland regions, many are adapted to tropical wetland environments and are crucial for their conservation (Luther and Greenberg Reference Luther and Greenberg2009). Through the lens of body-mass distributions, we aim to strengthen our understanding of tetrapod community structure in coastal tropical wetland ecosystems by examining similarities among some of the most ancient examples of coastal tropical wetland tetrapod communities and their modern counterparts.
Some of the earliest known, predominantly terrestrial, coastal tropical wetland ecosystems are found in the unique Pennsylvanian fossil forests of Joggins, Florence, and Point Aconi, Nova Scotia, Canada, aged in European chronostratigraphy to the Langsettian and Cantabrian Pennsylvanian substages, respectively (Waldron et al. Reference Waldron, Giles and Thomas2017). The Joggins and Point Aconi fossil localities record estuarine swamp environments with multi-tiered canopy structures (Zodrow Reference Zodrow2003; Calder et al. Reference Calder, Gibling, Scott, Davies, Hebert, Greb and DiMichele2006; Falcon-Lang Reference Falcon-Lang2006; Falcon-Lang et al. Reference Falcon-Lang, Benton, Braddy and Davies2006) oft compared to modern mangrove ecosystems (e.g., Raymond and Phillips Reference Raymond, Phillips and Teas1983; Greb et al. Reference Greb, DiMichele, Gastaldo, Greb and DiMichele2006; Falcon-Lang et al. Reference Falcon-Lang, DiMichele, Elrick and Nelson2009; Cleal et al. Reference Cleal, Uhl, Cascales-Miñana, Thomas, Bashforth, King and Zodrow2012). The coastal wetland forests of Pennsylvanian Nova Scotia hosted early representatives of modern amniote lineages, as well as numerous non-amniote tetrapods (Carroll Reference Carroll1964, Reference Carroll1966, Reference Carroll1967a; Reisz Reference Reisz1972; Grey and Finkel Reference Grey and Finkel2011; Pardo and Mann Reference Pardo and Mann2018; Maddin et al. Reference Maddin, Mann and Hebert2019; Mann et al. Reference Mann, Gee, Pardo, Marjanović, Adams, Calthorpe, Maddin and Anderson2020). The Late Pennsylvanian–Early Permian collapse and disappearance of coastal wetlands in Western Pangaea, most likely a result of widespread aridification and increased seasonal variability (Montañez et al. Reference Montañez, Tabor, Niemeier, DiMichele, Frank, Fielding, Isbell, Birgenheier and Rygel2007), is recognizable in the macrofloral record of Atlantic Canada (Whittingham et al. Reference Whittingham, Fraser and Maddin2026) and has been invoked as a primary driver of phylogenetic and functional diversification in amniotes and lissamphibians (Benton Reference Benton2010; Sahney et al. Reference Sahney, Benton and Falcon-Lang2010; Pardo et al. Reference Pardo, Small, Milner and Huttenlocker2019; So et al. Reference So, Kufner, Pardo, Edwards, Price, Bevitt and LeClair-Diaz2024). Environments like those of Joggins and Point Aconi therefore represent a major stepping stone in the evolution of terrestrial ecosystems and can provide a window into the origins of modern tropical wetland community structures.
To determine the degree of functional similarity between ancient and modern ecosystems, it is critical to use comparable metrics that are informative of organismal niches. Historically, examinations of functional diversity in Carboniferous tetrapods have focused on global-scale trends (e.g., Benton Reference Benton2010; Sahney et al. Reference Sahney, Benton and Falcon-Lang2010), or the characterization of specific functional traits in exceptionally preserved specimens (e.g., Reisz and Fröbisch, Reference Reisz and Fröbisch2014; Modesto et al. Reference Modesto, Scott, MacDougall, Sues, Evans and Reisz2015; Maddin et al. Reference Maddin, Mann and Hebert2019; Mann et al. Reference Mann, Dudgeon, Henrici, Berman and Pierce2021; Mann et al. Reference Mann, Xiong, Calthorpe, Sues and Maddin2026). However, fossil assemblages like those from the Carboniferous of Nova Scotia typically lack well-preserved specimens in abundances sufficient for the direct estimation of the community-wide diversity of traits like diet or locomotion. We therefore look to body-mass distributions as the most practical measurement for examining community structure, as body mass can be estimated using commonly preserved skeletal elements (Romer and Price Reference Romer and Price1940). In tetrapods, body mass covaries with many other functional traits, including locomotion, thermoregulation, life history, home range size, metabolic rate, and diet (Western Reference Western1979; Nagy et al. Reference Nagy, Girard and Brown1999; Haskell et al. Reference Haskell, Richie and Olff2002; Perry and Garland Reference Perry and Garland2002; Sibly and Brown Reference Sibly and Brown2007; McNab Reference McNab2008; Clarke and Pörtner Reference Clarke and Pörtner2010; Sandel Reference Sandel2013; Lovegrove and Mowoe Reference Lovegrove and Mowoe2013; Pineda-Munoz et al. Reference Pineda-Munoz, Evans and Alroy2016). This diversity of covarying traits makes body mass particularly useful in broad estimations and comparisons of functional diversity, community structure, and niche occupation (Grossnickle Reference Grossnickle2020; Slater Reference Slater2022). Community-scale body-mass distributions may therefore be impacted by a wide variety of factors relevant to community assembly and structure (Brown and Nicoletto Reference Brown and Nicoletto1991; Brown et al. Reference Brown, Marquet and Taper1993; Allen et al. Reference Allen, Garmestani, Havlicek, Marquet, Peterson, Restrepo, Stow and Weeks2006) and show patterns that reflect the diversity of other covariant traits. Differences in the body-mass distributions of modern and fossil communities are therefore likely to be reflective of overall differences in functional diversity.
While tetrapod body-mass distributions are often influenced by variations in climate and habitat (Hoekstra and Fagan Reference Hoekstra and Fagan1998; Silva et al. Reference Silva, Brimacombe and Downing2001), those relationships may not be especially consistent or widely applicable (Rodríguez Reference Rodríguez1999). Furthermore, there is an abundance of literature relating Pennsylvanian coal swamps to modern coastal environments (Falcon-Lang Reference Falcon-Lang2000; Falcon-Lang et al. Reference Falcon-Lang, Rygel, Calder and Gibling2004; Greb et al. Reference Greb, DiMichele, Gastaldo, Greb and DiMichele2006; Montañez Reference Montañez2016), allowing us to control for the climatic effects by comparing ancient communities to modern communities with similar abiotic conditions. There are also evident differences between fossil and modern coastal wetland communities, which are excluded from comparison here, such as the apparent absence of volant or fossorial tetrapods from the fossil assemblages and the absence of giant arthropods from modern communities (Greb et al. Reference Greb, DiMichele, Gastaldo, Greb and DiMichele2006). This work aims to assess the degree of differentiation among the shapes of body-mass distributions of Pennsylvanian and modern terrestrial tetrapod communities in coastal tropical wetlands and to determine whether any of taphonomy, dietary ecology, Bauplan, or phylogenetic affiliation can account for the differences observed.
Methods
Material for Body-Mass Estimation
To date, the only published attempt at estimating the body masses of Carboniferous tetrapods from Nova Scotia comes from Modesto et al. (Reference Modesto, Scott, MacDougall, Sues, Evans and Reisz2015: supplemental material), who inferred the body masses of four Carboniferous amniotes along with three unidentified fragmentary synapsids. The body-mass estimates of Modesto et al. (Reference Modesto, Scott, MacDougall, Sues, Evans and Reisz2015) were based on femoral and humeral measurements using equations derived from modern tetrapods by Campione and Evans (Reference Campione and Evans2012). While Campione and Evans’s (Reference Campione and Evans2012) methods are robust, the Carboniferous tetrapods from Nova Scotia rarely preserve sufficient limb elements. Therefore, to obtain a suite of body-mass estimates more thoroughly representative of the fossil assemblages at Joggins and Point Aconi, it is necessary to look to alternative proxies for body-mass estimation.
To acquire fossil body-mass distributions, we expanded upon Romer and Price’s (Reference Romer and Price1940) method of estimating body mass in fossil tetrapods from the centra of trunk vertebrae (defined as the vertebrae between the pectoral and pelvic girdles), or what Romer and Price referred to as dorsal vertebral centra. This correlation was established by estimating the mass of Dimetrodon limbatus, using the volume of a three-dimensional reconstructed model of a specimen of D. limbatus with a presumed specific gravity of 1.025, and then relating that mass estimate to the dimensions of the thoracic and lumbar vertebrae of the original fossil. While Romer and Price (Reference Romer and Price1940) never explicitly stated the equation used to relate their vertebral dimensions to body mass, it has since been deduced based on their published measurements of D. limbatus (N. Brocklehurst personal communication 2021) as:
where BM is the estimated body mass (in kilograms), and r is the mean of the radii of the trunk centra (in millimeters). Based on Romer and Price’s (Reference Romer and Price1940) reporting, it appears that they measured the diameters of the posterior ends of the trunk centra, halving that number to produce a radius value. Being based on a single modeled specimen of D. limbatus, it is unclear how broadly applicable equation (1) is. However, given that much of the stress of bearing and moving the trunk mass of terrestrial quadrupedal tetrapods is borne by the thoracic and lumbar parts of the axial skeleton (Rockwell et al. Reference Rockwell, Evans and Pheasant1938; Carrier Reference Carrier1987; Giffin Reference Giffin1990; Galbusera and Bassani Reference Galbusera and Bassani2019), a relationship between quadruped body mass and trunk vertebral dimensions should logically exist. Indeed, Ijima and Kubo (Reference Ijima and Kubo2019) were able to identify a scaling relationship between vertebral dimensions and body size in some modern crocodylians, though this was never translated into a method for estimating fossil body sizes. We therefore reexamined the relationship between the dimensions of trunk vertebrae and body mass using modern tetrapods to improve the robustness and broad applicability of this body-mass proxy.
To examine the relationship between the dimensions of trunk centra and body mass, the centra of thoracic and lumbar vertebrae were measured from an array of 81 modern tetrapod taxa (see Supplementary Table A1) from the comparative osteology collections of the Canadian Museum of Nature (CMN) and Royal Ontario Museum (ROM). Among the measured specimens were 32 mammals, 45 squamates, 2 crocodylians, and 2 lissamphibians. Measured specimens were all adults, primarily quadrupedal, non-volant, and either semiaquatic or fully terrestrial. For each specimen, wet body mass at time of collection was recorded. Radii of trunk centra were collected by halving the dorsoventral and lateral diameters of the anterior end of each centrum, measured using digital vernier calipers. In large centra or those with substantially protrusive lateral processes, the outer lateral extent of the centrum was considered to be the point of maximum relief with respect to the plane of the anterior end of the centrum (Fig. 1).
Diagram defining measurements taken for trunk vertebral centra. Shown is the anterior view of a thoracic vertebra of Iguana iguana, modified from Williston (Reference Williston1925). Both diameters were measured between the points of maximum relief with respect to the plane of the anterior end of the centrum.
From the above measurements, we calculated an ovate approximated surface area (OSA) for each centrum as the product of the lateral and dorsoventral radii. All trunk centra were measured for each measured specimen, and maximum, minimum, and mean values were determined for each array of OSA measurements taken from each specimen. As a result, for any given specimen, maximum OSA corresponds to the approximate surface area of the largest trunk centrum, minimum OSA corresponds to the approximate surface area of the smallest trunk centrum, and mean OSA is a mean of all approximated surface areas across all trunk centra. The full dataset is available in Supplementary Information A.
Twenty-one fossil tetrapod specimens were measured for the purpose of body-mass estimation, numbering 12 from the Langsettian Joggins locality and 9 from the Cantabrian Point Aconi locality, both in Nova Scotia. While the Cantabrian substage is now considered defunct by some authors (Nelson and Lucas Reference Nelson, Lucas, Lucas, Hunt and Lichtig2021; Nelson et al. Reference Nelson, Lucas, Elrick, Lucas, DiMichele, Opluštil and Wang2023), the Cantabrian section of the Sydney Mines Formation, which includes the Point Aconi and associated Florence localities, remains biostratigraphically distinct and has not yet been chronostratigraphically reassessed (Zodrow Reference Zodrow1989a,Reference Zodrowb; Cleal et al. Reference Cleal, Dimitrova and Zodrow2003; Gibling et al. Reference Gibling, Saunders, Tibert, White, Pashin and Gastaldo2004). For this work, fossils reported from the Cantabrian Florence locality of Cape Breton, Nova Scotia, are included in the Point Aconi assemblage, as the two fossil sites are less than 5 km apart and are from very similar horizons just above the Lloyd Cove coal seam (Reisz Reference Reisz1972; Maddin et al. Reference Maddin, Mann and Hebert2019). For each fossil species, the largest available representative of that species with a preserved axial skeleton was selected to minimize the potential accidental inclusion of juveniles. In some cases, vertebral elements were measured from specimens yet to be taxonomically identified. Where available, specimens were measured using digital vernier calipers from museum specimens provided by the Redpath Museum (RM), Canadian Museum of Nature (CMN), Nova Scotia Museum (NSM), and British Museum (BM). Specimens not accessible in-person were measured from figures provided in Carroll (Reference Carroll1964, Reference Carroll1966, Reference Carroll1967b, Reference Carroll1969), Reisz (Reference Reisz1972), and Klembara (Reference Klembara1985). Specimens measured from literature figures were measured in the image editing software ImageJ (Schneider et al. Reference Schneider, Rasband and Eliceiri2012). References for specimen images can be found in Supplementary Table A2, specimen images can be found in Supplementary Information H.
Surface areas of vertebral centra were approximated by calculating the product of the dorsoventral and lateral radii for each preserved thoracic and lumbar centrum in each fossil specimen. Maximum values of those surface area approximations are given in Supplementary Table A1.
Linear Modeling and Estimation of Fossil Body Masses
Following the methodology of Campione and Evans (Reference Campione and Evans2012), linear regression models of body masses and approximated surface areas of dorsal vertebral centra were fit using both ordinary least squares (OLS) and semi-major axes (SMA) in R v. 2024.04.2 (R Development Core Team 2024) in the R package smatr (Warton et al. Reference Warton, Duursma, Falster and Taskinen2012). OLS and SMA models were made between the base-10 logarithm of body masses (in grams) and the base-10 logarithm of mean, maximum, and minimum approximated centrum surface areas (in square millimeters) of the extant specimens (Supplementary Table A1). OLS regression was used to determine the most explanatory model among the vertebral dimension variables measured. In addition to comparing r 2 values among models, we also compared standard errors of the estimate (SEEs), Akaike information criterion corrected scores (AICc scores), Akaike weights, and percent prediction errors (PPEs), as the latter metrics are more robust than r 2 with respect to correlations with sample size and when comparing log-log relationships (Smith Reference Smith1984; Van Valkenburgh Reference Van Valkenburgh, Damuth and McFadden1990). Akaike weights were determined using the R package AICcmodavg (Mazerolle Reference Mazerolle2023). PPEs were calculated following Wu et al. (Reference Wu, Baraldo and Forlanut1995) by taking the average of percent errors produced by dividing absolute residuals by their corresponding measured log-transformed body-mass values. Because all linear regressions were performed on log-transformed explanatory and response variables, the OSA–body mass relationships correspond to power-law relationships, and not true linear relationships. Prediction ranges were calculated based upon log-transformed results and are therefore not equivalent to PPE values for untransformed linear models, as the error ranges for power-law relationships will broaden as OSA increases.
To test for phylogenetic autocorrelation in our linear models, we obtained a tip-dated molecular phylogeny of all 81 measured modern tetrapod species from TimeTree of Life 5 (Kumar et al. Reference Kumar, Suleski, Craig, Kasprowicz, Sanderford, Li, Stecher and Hedges2022). We assessed phylogenetic autocorrelation by calculating Pagel’s lambda scores with the R package phytools (Revell Reference Revell2024), using the phylogeny acquired from TimeTrees and the residuals of each OLS model.
Additionally, differences in the slopes and intercepts of each SMA model were compared among different subgroups (defined by clade, terrestriality, and body-mass category) of the measured dataset to determine differences in model fit among each of the tested subgroups. This analysis was done through the slopetest and elevtest functions in the smatr R package (Warton et al. Reference Warton, Duursma, Falster and Taskinen2012). For each of the SMA linear models, slope and intercept tests were made to compare between mammals and squamates, semiaquatic and terrestrial species, and specimens from five discrete body-mass categories. Body-mass categories were defined as <10 g, 10–100 g, 100 g to 1 kg, 1–10 kg, 10–100 kg, and >100 kg. Lissamphibian and crocodylian SMAs could not be compared as the two clades were represented by only three and two specimens in our dataset, respectively.
The best-fit and most generally applicable model was used to predict body masses of fossil tetrapods from Joggins and Point Aconi fossil tetrapods and to re-estimate the “pelycosaur” body masses reported by Romer and Price (Reference Romer and Price1940). While Romer and Price (Reference Romer and Price1940) only reported a single radius value for each centrum, the centra of their measured specimens were reported to be nearly circular in cross section, such that the square of the radius of the centrum should still be roughly equivalent to the OSA calculated herein. Detransformed body-mass predictions were corrected for smearing biases following Duan (Reference Duan1983) and Smith (Reference Smith1993). For the body-mass distributions of the Joggins and Point Aconi assemblages, summary statistics were collected in R using the packages moments (Komsta and Novometsky Reference Komsta and Novometsky2015) and multimode (Ameijeiras-Alonso et al. Reference Ameijeiras-Alonso, Crujeiras and Rodríguez-Casal2021).
Compilation of Modern Body Masses
For comparison with the estimated body-mass distributions of the Joggins and Point Aconi fossil assemblages, occurrences of terrestrial and semiaquatic reptiles and mammals were compiled from 29 tropical estuarine ecosystems using the Global Biodiversity Information Facility (GBIF). Selected modern locations were all between 15°S and 15°N, were centered on a coastal estuary or river delta, and had minimal impact from human development as determined based on Venter et al.’s (Reference Venter, Sanderson, Magrach, Allan, Beher, Jones and Possingham2016) Human Footprint map. Each of the 29 modern sites was defined by a 2500 km2 square centered around an estuary or delta and the freshwater body that feeds it. GBIF occurrence data of mammal and reptile species were then compiled within each 2500 km2 square, for which occurrence download records, including site coordinates, can be found in Supplementary Information B. Fully aquatic, volant, domestic, and extinct species were removed from each modern assemblage.
For each species in the modern occurrence dataset, species-wide maximum wild body mass was recorded to most closely approximate the specimen selection methods used for the fossil tetrapods. Reptile body masses were taken from the reptile body-mass database compiled by Slavenko et al. (Reference Slavenko, Tallowin, Itescu, Raia and Meiri2016: supporting information 2). Most mammal body masses were taken from the PanTHERIA database (Jones et al. Reference Jones, Bielby, Cardillo, Fritz, O’Dell, Orme and Safi2009). Mammal body masses were acquired from other sources in the cases of Neofelis diardi (Hearn et al. Reference Hearn, Ross, Pamin, Bernard, Hunter and MacDonald2013), Hylobates albibarbis (Cheyne Reference Cheyne, Gursky and Supriatna2010), Crocidura macarthuri (Mittermeier and Wilson Reference Mittermeier and Wilson2018), Pelomys campanae (Freudenthal and Martín-Suárez Reference Freudenthal and Martín-Suárez2013), Sylvisorex moro (Mittermeier and Wilson Reference Mittermeier and Wilson2018), Oecomys auyantepui (Wilson et al. Reference Wilson, Mittermeier and Lacher2017), Hyladelphys kalinowskii (Voss and Jansa Reference Voss and Jansa2009), Neacomys paracou (Voss et al. Reference Voss, Lunde and Simmons2001), Hylomyscus walterverheyeni (Nicolas et al. Reference Nicolas, Wendelen, Barriere, Dudu and Colyn2008), Leopardus tigrinus (MacDonald et al. Reference MacDonald, Loveridge, Nowell, Macdonald and Loveridge2010), and Galictis vittata (Yensen and Tarifa Reference Yensen and Tarifa2003), all species absent from the PanTHERIA database. The reptile species absent from the Slavenko et al. (Reference Slavenko, Tallowin, Itescu, Raia and Meiri2016) database also lacked published maximum body masses and so had to be excluded from analyses altogether. Broad dietary categories were assigned to each mammal and reptile species, generalized to herbivory, omnivory, and carnivory. For reptiles, dietary categories were pulled from the ReptTraits database (Oskyrko et al. Reference Oskyrko, Mi, Meiri and Du2024), and from the PanTHERIA database (Jones et al. Reference Jones, Bielby, Cardillo, Fritz, O’Dell, Orme and Safi2009) for mammals. Full data tables of body masses and relevant traits for the modern datasets are available in Supplementary Information C.
Taphonomic Modeling
Our current understanding of taphonomic trends suggests that small-bodied organisms are preferentially lost from the fossil record (Dodson Reference Dodson1971; Kidwell and Aigner Reference Kidwell, Aigner, Bayer and Seilacher1985; Behrensmeyer Reference Behrensmeyer1988, Reference Behrensmeyer, Allison and Briggs1991; Behrensmeyer et al. Reference Behrensmeyer, Kidwell and Gastaldo2000), and that expectation is a point of primary consideration when examining fossil body-mass distributions (e.g., Alroy Reference Alroy2003; Eberth et al. Reference Eberth, Xing and Clark2010; Brown et al. Reference Brown, Evans, Campione, O’Brien and Eberth2013). Additionally, many of the tetrapod fossils from Joggins and all tetrapod fossils from Point Aconi were found preserved in the hollowed-out stumps of arborescent lycopsids (Dawson Reference Dawson1894; Carroll Reference Carroll1966; Reisz Reference Reisz1972; Maddin et al. Reference Maddin, Mann and Hebert2019). These unique depositional environments might have incurred additional bias against the preservation of articulated tetrapods too large to fit within the stumps, which typically measured between 10 and 50 cm in diameter (Falcon-Lang et al. Reference Falcon-Lang, Benton, Braddy and Davies2006). To reasonably compare Joggins and Point Aconi to modern sites, it is then necessary to simulate multiple modes of taphonomic size bias in modern assemblages and determine which modes (if any) best improve the comparability of modern and fossil body-mass distributions.
To account for differences in sample size between the fossil and modern datasets, each site in the full modern dataset was randomly subsampled to diversities equivalent to those of the fossil assemblages using the sample function in base R without replacement. Modern sites were also subsampled down to sample size thresholds of 60, 40, and 20 species to account for type 2 error risks for Kolmogorov-Smirnov (KS) tests based on small samples. Each modern locality was randomly subsampled 100 times, with each replicate being compared with the fossil assemblage. Three different variations of this subsampling approach were constructed to test for the effects of taphonomic size bias on the body-mass distributions of the fossil assemblages. To simulate bias against smaller-bodied species, modern sites were subsampled such that each species with a body mass less than the estimated lower bound of the smallest animal from either fossil assemblage (Ricnodon sp.) had only a 5% probability of being sampled, with each remaining species having a 95% probability of being sampled. The same bias was applied in reverse to simulate bias against the preservation of larger-bodied species, with the 5% probability of being sampled being applied to each species with a body-mass greater than the estimated upper bound of the largest animal from either fossil assemblage (Carbonoherpeton carrolli). Finally, bias against both large and small extremes of body size was simulated by combining both sampling probability restrictions. Sampling probability thresholds were based on the maximum and minimum estimated fossil body masses to best approximate bias against the preservation of body masses beyond the range of those recorded from either fossil assemblage.
In addition, alternate versions of the modern occurrence data were also randomly subsampled (with and without accounting for taphonomic size bias) and tested for differences in similarity to the full modern assemblages. To imitate the effect of allochthony in the studied fossil assemblages, additional modern datasets with wider geographic ranges were accessed from GBIF. Both fossil tetrapod assemblages studied herein have been considered to be autochthonous or para-autochthonous (Falcon-Lang et al. Reference Falcon-Lang, Benton, Braddy and Davies2006; Maddin et al. Reference Maddin, Mann and Hebert2019), but fluvially fed estuarine fossil assemblages can be highly allochthonous (Behrensmeyer and Hook Reference Behrensmeyer, Hook, Behrensmeyer, Damuth, DiMichele, Potts, Sues and Wing1992). Increasing the spatial scale of sampling can have variable effects on body-mass distributions (Brown and Nicoletto Reference Brown and Nicoletto1991; Arita and Figueroa Reference Arita and Figueroa1999; Bakker and Kelt Reference Bakker and Kelt2000; Coetzee et al. Reference Coetzee, le Roux and Chown2012). We therefore investigated whether altering the geographic scale of sampling for modern sites increases comparability with the fossil assemblages. These larger sites were selected based on the same criteria as the twenty-nine 2500 km2 squares, but were defined by twelve 10,000 km2 and twelve 40,000 km2 squares. Each of the twelve 10,000 km2 squares encompasses a separate preexisting 2500 km2 square, with each 40,000 km2 square likewise encompassing one of the 10,000 km2 squares. Geographic ranges were designed to include as much of the upstream watershed as possible, primarily extending inland relative to their smaller counterparts. Full data tables of body masses and relevant traits for the broadened datasets are available in Supplementary Information D. At each level of geographic scale, all 12 modern sites were iteratively randomly subsampled and collectively tested for relative degree of fit with the fossil body-mass distributions in the same way as outlined for the tests of fit with taphonomic size-bias simulations.
Tests of Fit between Modern and Fossil Body-Mass Distributions
Similarity between the body-mass distributions of fossil and modern tetrapod assemblages was primarily determined via KS test (Chakravarti et al. Reference Chakravarti, Laha and Roy1967). The KS test produces a D-statistic, representing the degree of dissimilarity between the two tested distributions, and a p-value, representing the probability that the two datasets are sampled from the same distribution, effectively providing a goodness-of-fit test. KS tests were performed using the R package Matching (Sekhon Reference Sekhon2011), each test being bootstrapped 1000 times. Both of the two estimated fossil body-mass distributions were compared via KS test to 100 random subsamples of the full range of body masses from each modern site, and then against various taphonomic subsampling simulations (as described earlier). For both fossil assemblages, the arrays of D-statistics from the KS tests against each subsampling method were compared amongst each other via multiple Mann-Whitney U-test in PAST 4.13 (Hammer et al. Reference Hammer, Harper and Ryan2001) to determine whether there were significant differences in the results of each set of tests. For all comparisons between multiple arrays of D-statistics, significance of Mann-Whitney U-tests was subjected to sequential Bonferroni correction to avoid multiple-testing biases (Holland and Copenhaver Reference Holland and Copenhaver1987). Along with KS tests, medians were collected from the subsampled modern sites for comparison with the fossil distributions, as well as test statistics from Anscombe-Glynn tests of kurtosis (Anscombe and Glynn Reference Anscombe and Glynn1983) and D’Agostino-Pearson tests of skewness (D’Agostino and Pearson Reference D’Agostino and Pearson1973). Test statistics were acquired using the moments (Komsta and Novometsky Reference Komsta and Novometsky2015) R package. The distributions of the summary statistics obtained from the fossil distributions were then compared with the arrays of summary statistics obtained from the variously subsampled modern distributions by single-case test in PAST 4.13. These methods allowed for the determination of the subsampling approach most appropriate for accounting for potential taphonomic size bias in each fossil assemblage.
KS tests were also conducted on various subsets of the primary full modern dataset to test for the effects of cladistic, physiological, and dietary diversity on the goodness of fit between fossil and modern body-mass distributions. Dietary effects on body-mass distributions can often be complex. Modern mammals exhibit some partitioning of body-mass space along dietary lines, with small (<1 kg) body sizes commonly being occupied by omnivores, herbivores, and insectivorous carnivores (Pineda-Munoz et al. Reference Pineda-Munoz, Evans and Alroy2016), while concentration of large body sizes may reflect a higher abundance of non-insectivorous carnivores (Hoekstra and Fagan Reference Hoekstra and Fagan1998). Similarly, herbivory in squamates has been associated with increased body size (Pough Reference Pough1973; Cooper and Vitt Reference Cooper and Vitt2002), though this association is not consistent across taxonomy and climate (Van Damme Reference Van Damme1999; Espinoza et al. Reference Espinoza, Wiens and Tracy2004; Vitt Reference Vitt2004). In Paleozoic tetrapods, the evolution of herbivory is hypothesized to have driven increased tetrapod body masses (Reisz and Fröbisch Reference Reisz and Fröbisch2014; Brocklehurst Reference Brocklehurst2016). If the Joggins and Point Aconi fossil assemblages lack herbivores, or are generally less dietarily diverse, they may better fit body-mass distributions of modern assemblages restricted to exclude herbivores and omnivores.
When considering taxonomic affinity, various mammal clades (Brown et al. Reference Brown, Marquet and Taper1993; Kelt Reference Kelt1998; Bakker and Kelt Reference Bakker and Kelt2000), as well as snakes (Boback and Guyer Reference Boback and Guyer2003), have been proposed to show unique peaks in their body-mass distributions, though the validity and ubiquity of clade-specific body-mass optima have been repeatedly challenged (Roy et al. Reference Roy, Jablonski and Martien2000; Meiri Reference Meiri2008; Raia et al. Reference Raia, Carotenuto and Meiri2010; Cox et al. Reference Cox, Boback and Guyer2011). Such patterns imply that phylogenetic history may constrain the evolution of body masses, producing shared optima among closely related clades that share similar shapes and life histories (Allen et al. Reference Allen, Garmestani, Havlicek, Marquet, Peterson, Restrepo, Stow and Weeks2006). There may be further constraint toward smaller body mass for abducted-limb quadrupeds like the amniotes and microsaurs found at Joggins and Point Aconi, as narrower, more vertical-oriented gaits have been proposed to have evolved to support larger amniote body sizes, based on evidence from trackways (Buchwitz et al. Reference Buchwitz, Jansen, Renaudie, Marchetti and Voigt2021). If the fossil assemblages are, in general, more closely comparable to mammalian or reptilian amniotes, they may better fit body-mass distributions of modern assemblages consisting only of mammals or only of reptiles. Further constraining modern reptile assemblages to just those with similar Bauplans to the fossil tetrapods may yet reveal whether there might be general, constrained optima in the body masses of squat, quadrupedal vertebrates.
To test for differences associated with dietary diversity, summary statistics and KS test results were compared among the complete set of modern body-mass distributions and two subsets of those distributions, one which included only species recorded as omnivorous or carnivorous, and one which included only species recorded as carnivorous. Differences associated with cladistic diversity were tested for by comparing summary statistics and KS test results among those produced by the complete modern body-mass distributions and subsets with either exclusively contained mammals or exclusively contained reptiles. Finally, the results of summary statistics and KS tests using the complete, mammalian, and squamate modern body-mass distributions were compared with those from a highly restrictive subset of only quadrupedal squamates and crocodylians with abducted-limb Bauplans. All subsets were subjected to random iterative subsampling simulating the best-fitting model of taphonomic size bias. In some cases, subsetting reduced some modern sites to fewer occurrences than those recorded from the fossil sites, preventing subsampling. In those cases, the affected sites were removed from that comparative subset to ensure consistency in sampling method. This issue also prevented simulations of sites subsampled to mimic multiple concurrent restrictions (e.g., sites subsampled to have only carnivores with abducted-limb Bauplans), as including multiple rounds of subsampling resulted in there being too few modern sites with sufficient diversity for comparison. As a result, combinatorial effects of diet and cladistic affinity could not be tested for. Test results were collected for all subsets at multiple modern subsample size thresholds as had been done for the taphonomic subsampling. Some subsets lacked sites with sufficient species diversity for all levels of subsampling to be applied. For these subsets, test results from those higher sample size thresholds could not be obtained. Specifically, a subsample size threshold of 60 species could not be met for the subsets consisting of only carnivores, only mammals, or only tetrapods with abducted-limb Bauplans. The subset of tetrapods with only abducted-limb Bauplans also lacked sites with high enough diversity to be subsampled at a sample size of 40 species.
Beyond cladistic affiliation, diet, and Bauplan, there are multiple additional factors that may be responsible for the differentiation we observe between the body-mass distributions of the Point Aconi and Joggins communities and those of the modern. Functional characteristics like locomotor mode, home range, and thermoregulation, among others, can be variably collinear with body mass (e.g., Reiss Reference Reiss1988; Kelt and Van Vuren Reference Kelt and Van Vuren1999; Lovegrove and Mowoe Reference Lovegrove and Mowoe2013; Sandel Reference Sandel2013; Kohli and Rowe Reference Kohli and Rowe2019), are generally expected to have been more constrained among late Carboniferous tetrapod species (Dunne et al. Reference Dunne, Close, Button, Brocklehurst, Cashmore, Lloyd and Butler2018; Buchwitz et al. Reference Buchwitz, Jansen, Renaudie, Marchetti and Voigt2021; Lucas Reference Lucas2022; Faure-Brac et al. Reference Faure-Brac, Woodward, Aubier and Cubo2024, but see Sumida and Modesto Reference Sumida and Modesto2001; Grigg et al. Reference Grigg, Nowack, Bicudo, Bal, Woodward and Seymour2022), and therefore may contribute to differences in body-mass distributions. Consequently, species interactions within Carboniferous tetrapod communities are likewise expected to differ substantially from the those of the modern (Olson Reference Olson1966; Otoo Reference Otoo2025), and may also have resulted in body-mass distributions that lack modern analogues. The impact of these additional factors on our results could not be assessed, as comprehensive data on additional functional traits and species interactions were not available for the sampled modern communities beyond diet and body mass.
Furthermore, while the depositional environments of the Joggins and Point Aconi localities evidence abiotic conditions mostly analogous to modern mangrove swamps or other low-latitude estuarine environments (Greb et al. Reference Greb, DiMichele, Gastaldo, Greb and DiMichele2006; Falcon-Lang et al. Reference Falcon-Lang, DiMichele, Elrick and Nelson2009; Cleal et al. Reference Cleal, Uhl, Cascales-Miñana, Thomas, Bashforth, King and Zodrow2012), the floral communities of those habitats were dominated by plants with life histories and physiologies that strongly differ from those of any modern forests (Wilson et al. Reference Wilson, White, Montañez, DiMichele, McElwwain, Poulsen and Hren2020; Whittingham et al. Reference Whittingham, Fraser and Maddin2026). Catastrophic changes in plant community structure in the late Carboniferous and early Permian have been interpreted as being instrumental in shaping the functional diversity of tetrapods (Benton Reference Benton2010; Sahney et al. Reference Sahney, Benton and Falcon-Lang2010; Brocklehurst and Benton Reference Brocklehurst and Benson2021), in turn impacting tetrapod body-mass distributions, though these effects cannot be readily isolated from concurrent changes in global climate and the shrinking or loss of wetland ecosystems (Dunne et al. Reference Dunne, Thompson, Butler, Rosindell and Close2023; Whittingham Reference Whittingham2025). While we did implement controls on abiotic conditions through restricting geographic criteria for selecting modern sites, there is no way to restrict modern communities in a way that controls for the differences between modern and late Carboniferous plant community structure, rendering hypotheses regarding the influence of floral communities on the body-mass distributions of Carboniferous tetrapod communities untestable within the methodological framework of this study.
R code for all analyses is available through the GitHub link provided in Supplementary Information I.
Results
Body Mass-Vertebral Centra Relationships
Linear regression of log10-transformed body mass against log10-transformed approximations of trunk centrum surface areas revealed strong positive power-law relationships between all centrum measurements and body mass (Table 1), with all r 2 values exceeding 0.93 and no mean PPEs greater than 12% of log10-transformed body mass. OLS results revealed that mean OSA has the strongest correlation with body mass, with an r 2 of 0.971, AICc of −22.114, and PPE of 9.1 ± 2.6% log10 grams. Maximum OSA also showed strong correlations with body mass with an r 2 value of 0.968, AICc of −14.157, and PPE of 9.3 ± 2.6% log10 grams. For untransformed body masses, PPEs of log-transformed data translate to median prediction errors of 30.2% for mean OSA, 32.23% for maximum OSA, and 34.3% for minimum OSA. Akaike weights heavily favored the mean OSA–body mass relationship (Table 1).
Results of ordinary least squares (OLS) regressions comparing log10-transformed body mass (in grams) with log10-transformed ovate approximations of trunk centrum surface area (OSA; in mm2). For each regression, slopes, intercepts, and r-squared values are given, along with standard errors of the estimate (SEE), Akaike information criterion corrected scores (AICc), Akaike weights percent prediction errors (PPE), and Pagel’s lambda values. The 95% confidence intervals for PPEs and p-values for Pagel’s lambda tests of phylogenetic autocorrelation are also given
Table 1. Long description
The table contains 11 columns and 4 rows. The columns are: Test variable, Intercept, Slope, r-squared, S E E, A I C c, Akaike weight, P P E transformed, P P E untransformed, Pagel’s lambda, and p Pagel’s lambda.
Row 1: Mean O S A. Intercept 1.896968666, Slope 1.388818494, r-squared 0.9707, S E E 0.20560, A I C c minus 22.114, Akaike weight 0.9816, P P E transformed 0.09125, P P E untransformed 0.30537, Pagel’s lambda 0.517600, p 0.016158.
Row 2: Minimum O S A. Intercept 2.057136694, Slope 1.395533174, r-squared 0.9539, S E E 0.25847, A I C c 15.098, Akaike weight 8.15 times 10 super minus 9, P P E transformed 0.10624, P P E untransformed 0.34331, Pagel’s lambda 0.606163, p 0.000008.
Row 3: Maximum O S A. Intercept 1.750874791, Slope 1.378311132, r-squared 0.9677, S E E 0.21624, A I C c minus 14.157, Akaike weight 0.0184, P P E transformed 0.09315, P P E untransformed 0.32232, Pagel’s lambda 0.000073, p 1.000000.
While OLS results showed the strongest correlations between mean OSA and body mass, SMA regression showed significant (p-value = 0.004) variation in the intercept of the mean OSA–body mass relationship between mammals and reptiles (Fig. 2A–C). This differentiation between mammals and reptiles is even stronger (p-value < 0.0001) in the intercept of the minimum OSA–body mass relationship (Table 2). Of the three OSA metrics, only the maximum OSA showed no significant deviations in its relationship with body mass among the tested subgroups (Table 2) and showed the least difference in residuals among the tested subgroups (Fig. 3). The linear regression equation for estimated body mass and maximum OSA is as follows:
where OSAmax is the maximum of all products of each trunk centrum’s lateral and dorsoventral radii in square millimeters. Equation (2) represents the strongest tested relationship between the dimensions of trunk vertebrae and body mass to show no significant deviations with respect to clade, terrestriality, or body-mass category (Fig. 2D–F). This relationship between raw body mass and maximum OSA without log transformation can also be transformed into a power law, which appears as follows after applying a smearing correction of 1.119561 (sensu Smith, Reference Smith1993):
Semi-major axes (SMA) relating log10-transformed body mass (in log10 grams) to mean (A–C) and maximum (D–F) log10-transformed ovate approximated surface area (OSA: in log10 mm2) of trunk centra in 82 modern tetrapods. A and D display OSA–body mass relationships grouped by clade, with mammals in blue and reptiles in yellow. As only three amphibians were included in the analysis, they cannot be included in SMA comparison between clades, and their data points are omitted from graphs A and D. B and E display OSA–body mass relationships grouped by terrestriality, with semiaquatic animals in blue and fully terrestrial animals in yellow. C and F display body-mass relationships grouped by body-mass category, represented in dark blue (<10 g), light blue (10–100 g), red (100 g to 1 kg), green (1–10 kg), gray (10–100 kg), and yellow (>100 kg). In all cases, the best-fit line for the full dataset is displayed in black.
Figure 2. Long description
A six-panel grid of scatter plots labeled A through F. All plots share a common Y-axis of Log B M (body mass in grams) ranging from 0 to 5 and an X-axis of Log surface area (O S A in millimeters squared) ranging from negative 1.0 to 2.5. Every plot shows a tight, positive linear relationship with a black best-fit line running diagonally from the bottom-left to the top-right.
* Top Row (Mean O S A):
- Panel A: Data points grouped by clade. Blue circles (mammals) and yellow circles (reptiles) are interspersed along the regression line.
- Panel B: Data points grouped by terrestriality. Blue circles (semiaquatic) and yellow circles (terrestrial) show overlapping distributions.
- Panel C: Data points grouped by mass categories. Colors transition from dark blue at the bottom-left (less than 10 grams) through light blue, red, green, and gray, to yellow at the top-right (greater than 100 kilograms).
* Bottom Row (Maximum O S A):
- Panel D: Similar to A, showing Max O S A grouped by clade with blue and yellow points.
- Panel E: Similar to B, showing Max O S A grouped by terrestriality with blue and yellow points.
- Panel F: Similar to C, showing Max O S A grouped by mass categories with the same color-coded gradient from dark blue to yellow.
In all panels, the individual regression lines for the subgroups closely track the primary black best-fit line, indicating a consistent scaling relationship across all biological categories.
Results of semi-major axes (SMA) regressions comparing log10-transformed body mass (BM; in grams) with log10-transformed ovate approximations of trunk centrum surface area (OSA; in mm2). Values depicted are p-values produce by testing for differences in slope and intercept between multiple subgroups for each regression. Subgroups are divided by clade (mammals vs. reptiles), terrestriality (terrestrial vs. semiaquatic), and body mass category (<10 g, 10–100 g, 100 g to 1 kg, 1–10 kg, 10–100 kg, >100 kg)
Table 2. Long description
The table consists of seven columns and four rows.
Column headers from left to right are: Test variable, Clade slopes p, Terrestriality slopes p, B M category slopes p, Clade intercepts p, Terrestriality intercepts p, and B M category intercepts p.
Row 1: Mean O S A. Values are 0.122467004, 0.278401793, 0.982835507, 0.004127, 0.726750922, and 0.736150343.
Row 2: Minimum O S A. Values are 0.699044137, 0.651696477, 0.904979766, 2.554 times 10 super minus 15, 0.332100471, and 0.124487296.
Row 3: Maximum O S A. Values are 0.495930689, 0.084153529, 0.9705266, 0.3923, 0.948495566, and 0.23198928.
Box-and-whisker plots of residuals organized by clade for the regressions of log10-transformed ordinary least square (OLS) of body mass against mean (A), maximum (B), and minimum (C) ovate approximated surface area (OSA).
Figure 3. Long description
A three-panel multi-plot showing residuals on the y-axis and clades on the x-axis. Each panel includes a horizontal zero line and three box plots for Lissamphibia, Mammalia, and Reptilia. Each box plot features a central median line, a red dot representing the mean, and whiskers extending to outliers.
Panel A, Mean O S A. Lissamphibia shows a very tight box entirely below the zero line at approximately negative 0.3. Mammalia is centered above the zero line with a mean near 0.1. Reptilia is centered on the zero line with a wider distribution.
Panel B, Maximum O S A. Lissamphibia is below the zero line with a mean near negative 0.2. Mammalia is centered just above the zero line. Reptilia has the highest mean and largest vertical spread, extending from negative 0.6 to 0.4.
Panel C, Minimum O S A. Lissamphibia is below the zero line at approximately negative 0.4. Mammalia is entirely above the zero line with a mean near 0.2. Reptilia is centered slightly below the zero line with a mean near negative 0.1.
Prediction errors for body-mass estimates for equation (3) can be produced as follows:
where BME is the prediction error to be added and subtracted from the body-mass estimate to produce upper and lower bounds for the estimate. Errors on untransformed body masses are therefore logarithmically scaled, meaning that PPEs for log-transformed body masses are preferable for use in calculating error ranges.
Pagel’s lambda tests only failed to reveal evidence of phylogenetic autocorrelation in the relationship between maximum OSA and body mass (Table 1), which showed a Pagel’s lambda of 0.0000733 (p-value = 1.0). The residuals of the relationships for both mean and minimum OSA were shown to be highly phylogenetically autocorrelated, with Pagel’s lambda values of 0.5176 (p-value = 0.016) and 0.6062 (p-value < 0.001), respectively.
Fossil Body-Mass Estimates
Equation (3) was used to estimate body masses of the 12 fossil specimens from Joggins and 9 fossil specimens from Point Aconi, each with an error range of ±9.3% log10 grams, based on the mean prediction error value obtained for the regression of maximum OSA and body mass (Table 3). The same equation and error range were also used to produce new body-mass estimates for the species measured by Romer and Price (Reference Romer and Price1940) (Table 4). A test of modality (sensu Ameijeiras-Alonso et al. Reference Ameijeiras-Alonso, Crujeiras and Rodríguez-Casal2019) revealed that the log10 distribution of body masses estimated from the Joggins tetrapods is two-peaked (p-value [greater than 1 peak] = 0.008, p-value [greater than 2 peaks] = 0.522), while a D’Agostino-Pearson test (D’Agostino and Pearson Reference D’Agostino and Pearson1973) showed no skewness (skewness = 0.0658, p-value = 0.902), and an Anscombe-Glynn test (Anscombe and Glynn Reference Anscombe and Glynn1983) showed low kurtosis (k = 1.70). Median body mass at the Joggins locality is 0.377 kg. The same tests showed that the log10 distribution of body masses estimated from the Point Aconi tetrapods is single peaked (p-value [greater than 1 peak] = 0.952), is not significantly skewed (skewness = −0.578, p-value = 0.319), has near-normal kurtosis (k = 3.04), and a median body mass of 1.164 kg. A bootstrapped KS test revealed no significant differences (D-statistic = 0.417, p-value = 0.251) in the body-mass distributions of the two fossil assemblages (Fig. 4).
Body masses of fossil tetrapods from the Pennsylvanian of Nova Scotia estimated using maximum ovate approximated surface area (OSA). Estimated masses are given along with error ranges derived from the equation (2) percent prediction error of 9.4% applied to the log10-transformed body mass (BM) estimate. Specimens measured from the literature are marked with an asterisk (*). For Hylonomus lyelli, Protoclepsydrops haplous, Echinerpeton intermedium, and Archaeothyris florensis, BM estimates made by Modesto et al. (Reference Modesto, Scott, MacDougall, Sues, Evans and Reisz2015) are also given
Table 3. Long description
The table consists of six columns: Specimen number, Locality, Name, B M estimate in kilograms, B M estimate range in kilograms, and Modesto et al. 2015 estimate in kilograms.
Data for Joggins locality:
- C M N 10041, Steenerpeton silvae: 0.475 kg (range 0.46 to 0.491).
- C M N 10048, Hylonomus lyelli: 0.015 kg (range 0.015 to 0.016), Modesto estimate 0.1.
- N S M 987 G F 99.1, Dendrerpeton acadianum: 2.575 kg (range 2.47 to 2.684).
- B M R.5775, Trachystegos megalodon: 1.165 kg (range 1.122 to 1.209).
- N S M 994 G F 1.1, Calligenethlon watsoni: 8.810 kg (range 8.396 to 9.246).
- R M 12056, Archerpeton anthracos: 0.039 kg (range 0.038 to 0.04).
- R M 12202, Protoclepsydrops haplous: 0.789 kg (range 0.761 to 0.817), Modesto estimate 10 to 100.
- R M 12091, Asaphestera intermedia: 0.370 kg (range 0.359 to 0.382).
- B M R.477, Ricnodon species: 0.013 kg (range 0.013 to 0.013).
- R M 21128, Novascoticus multidens: 0.021 kg (range 0.021 to 0.021).
- R M 21129, Indet.: 0.627 kg (range 0.606 to 0.648).
- R M 21191, Leiocephalikon problematicum: 0.034 kg (range 0.033 to 0.035).
Data for Point Aconi locality:
- M C Z 4082, Archaeothyris florensis: 1.304 kg (range 1.255 to 1.354), Modesto estimate 2.
- M C Z 4090, Echinerpeton intermedium: 0.590 kg (range 0.57 to 0.61), Modesto estimate 0.25.
- M C Z 4088, Indet.: 0.789 kg (range 0.761 to 0.817).
- M C Z 3481, 3482, Paleothyris acadiana: 0.392 kg (range 0.38 to 0.405).
- M C Z 3034, Limnostygis relictus: 2.881 kg (range 2.762 to 3.005).
- M C Z 2773, Carbonoherpeton carrolli: 13.475 kg (range 12.812 to 14.173).
- N S M 017 G F 020.001, Dendromaia unimakiensis: 5.330 kg (range 5.093 to 5.578).
- N S M 017 G F 020.022, Indet. amniote: 0.040 kg (range 0.04 to 0.041).
- N S M 017 G F 020.004, Indet. tetrapod: 1.994 kg (range 1.915 to 2.076).
Mean body masses of tetrapod species measured by Romer and Price (Reference Romer and Price1940) estimated using maximum ovate approximated surface area (OSA). Estimated masses are given along with error ranges derived from the equation (2) percent prediction error of 9.3% applied to the log10-transformed body mass (BM) estimate. Romer and Price’s mean BM estimates are also provided for each species
Table 4. Long description
The table consists of four columns and 28 rows of species data. The columns are labeled: Species, Max. O S A B M estimate (kg), Max. O S A B M estimate range (kg), and Romer B M estimate.
Key data points include:
* Aerosaurus greenleeorum: 8.81 kg (range 8.4 to 9.25), Romer estimate 17.
* Baldwinonus trux: 33.56 kg (range 31.75 to 35.47), Romer estimate 52.
* Casea broilii: 12.29 kg (range 11.69 to 12.92), Romer estimate 21.
* Clepsydrops vinslovii: 137.29 kg (range 128.91 to 146.22), Romer estimate 14.
* Dimetrodon grandis: 339.36 kg (range 317.05 to 363.23), Romer estimate 254.
* Dimetrodon natalis: 17.02 kg (range 16.16 to 17.92), Romer estimate 28.
* Edaphosaurus pogonias: 248.83 kg (range 232.88 to 265.88), Romer estimate 186.
* Mycterosaurus longiceps: 1.46 kg (range 1.4 to 1.51), Romer estimate 3.
* Ophiacodon major: 268.94 kg (range 251.62 to 287.47), Romer estimate 230.
* Sphenacodon ferocior: 95.94 kg (range 90.26 to 101.99), Romer estimate 83.
The O S A estimates generally show lower values for smaller species and higher values for larger species compared to the Romer estimates, though significant discrepancies exist, such as in Clepsydrops vinslovii where the O S A estimate is nearly ten times higher than the Romer estimate.
Frequency distributions of log10-transformed body masses (in log10 kg) for the Joggins (black) and Point Aconi (red) fossil assemblages. Log10-transformed body-mass estimates for individual specimens are shown for specimens from Joggins (bottom, black) and Point Aconi (top, red). The p-value displayed represents the result of a bootstrapped Kolmogorov-Smirnov test comparing the two distributions.
Figure 4. Long description
The graph features a title at the top: Estimated log B M s from Joggins vs Point Aconi (p-value: 0.251). The x-axis is labeled Log B M and ranges from negative 4 to 4. The y-axis is labeled Density and ranges from 0.0 to 1.0.
Two density curves are shown:
- A black curve representing Joggins shows a bimodal distribution with a primary peak near 0 and a secondary, lower peak near negative 1.5.
- A red curve representing Point Aconi shows a taller, narrower primary peak centered just below 0 and a much smaller secondary peak near negative 1.5.
Rug plots at the top and bottom indicate individual specimen data points:
- Red vertical ticks for Point Aconi are clustered at the top of the plot between negative 1.5 and 1.0.
- Black vertical ticks for Joggins are clustered at the bottom of the plot between negative 2.0 and 1.0.
The curves overlap significantly between negative 2 and 2 on the x-axis.
Comparison of Fossil and Modern Body-Mass Distributions
Comparison of D-values resulting from KS tests between the fossil assemblages and variously subsampled full modern assemblages showed significantly different results between some subsampling modes (Supplementary Table F1). Subsampling the full modern sites at the same sample sizes as the fossil sites, without simulating size bias, produced body-mass distributions that significantly differed from the Joggins assemblage in 26.9% of repetitions (mean D-statistic = 0.504) and significantly differed from Point Aconi in 7.0% of repetitions (mean D-statistic = 0.377). Raw modern body-mass distributions differed from Joggins in 36.1% (mean D-statistic = 0.477), 48.8% (mean D-statistic = 0.474), and 57.1% of repetitions (mean D-statistic = 0.462), and differed from Point Aconi in 8.2% (mean D-statistic = 0.331), 2.1% (mean D-statistic = 0.292), and 1.2% of repetitions (mean D-statistic = 0.284), when subsampled at sample sizes of 20, 40, and 60, respectively. Direct site-by-site graphical comparison of fossil body-mass distributions with unaltered modern body-mass distributions is available in Supplementary Information E. Accounting for taphonomic size bias significantly increased goodness of fit between the Joggins assemblages and the modern (Fig. 5; Supplementary Table G). Simulating taphonomic size bias in the modern datasets revealed that the Joggins assemblage was best fit by simulations of size bias against small (<12.6 g) vertebrates, which significantly differed from the fossil distributions in 18.9% of repetitions (mean D-statistic = 0.473). Taphonomic simulations did not improve the fit of the Point Aconi assemblage to modern distributions of any sample size. KS-test results are mirrored by comparison of median body masses between taphonomically biased modern distributions and Point Aconi, where only the unaltered modern distributions show median body masses consistent with that of Point Aconi at all sample sizes (Supplementary Information F). Based on these best-fit models of size bias, taphonomic bias against small animals was simulated for all comparisons with Joggins, while no taphonomic bias was simulated for comparisons with Point Aconi.
Box-and-whisker plots of D-statistics produced by Kolmogorov-Smirnov tests of goodness of fit between fossil body-mass distributions (Joggins in A, C, E, G; Point Aconi in B, D, F, H) and modern body-mass distributions altered to simulate taphonomic size bias. Modern distributions tested were subsampled at sample sizes equivalent to fossil distributions (A, B) and at increased sample sizes of 20 (C, D), 40 (E, F), and 60 species (G, H). Lower D-statistic values represent greater similarity between the simulated and fossil datasets. Red dots represent mean values for each array of D-statistics.
Figure 5. Long description
The figure consists of eight panels arranged in a two-column by four-row grid. The left column (A, C, E, G) is titled Modern Simulations versus Joggins, and the right column (B, D, F, H) is titled Modern Simulations versus Point Aconi.
Each panel shares a common X-axis labeled Taphonomic Simulation with four categories: Unbiased, Bias Against Small Size, Bias Against Large Size, and Bias Against Extremes. The Y-axis represents the D-Value, with scales varying by row. Each plot contains four vertical box-and-whisker elements. A red dot within each box indicates the mean value, while a horizontal black line indicates the median. Outliers are represented by small black dots above and below the whiskers.
* Row 1 (A and B): Sample sizes equivalent to fossil distributions. In A, Bias Against Large Size shows the highest D-values (mean near 0.6), while Bias Against Small Size is lowest (mean near 0.4). In B, values are more compressed between 0.3 and 0.5.
* Row 2 (C and D): Increased sample size of 20 species. Trends remain similar to Row 1 but with slightly tighter distributions.
* Row 3 (E and F): Increased sample size of 40 species. In E, the mean for Bias Against Small Size drops below 0.4. In F, Bias Against Small Size shows a higher mean than the other three categories.
* Row 4 (G and H): Increased sample size of 60 species. In G, Bias Against Large Size remains the highest outlier with a mean near 0.55. In H, all categories show means clustered between 0.28 and 0.33.
Across all Joggins panels, Bias Against Large Size consistently produces the highest D-values, indicating the least similarity to the fossil data. In Point Aconi panels, the differences between simulation types are less pronounced than in the Joggins panels.
After accounting for taphonomic size bias, comparison of the Joggins and modern assemblages of increasing geographic scale failed to reveal any significant improvements in goodness of fit with scale compared to the 2500 km2-scale assemblages (Fig. 6). In comparison with Joggins, the 12 modern sites sampled from 2500 km2, 10,000 km2, and 40,000 km2 squares significantly differed in 4.8% (mean D-statistic = 0.402), 7.4% (mean D-statistic = 0.401), and 4.9% (mean D-statistic = 0.387) of repetitions, respectively, patterns that are reflected at each sample size. By contrast, in comparison with Point Aconi, modern sites sampled from 2500 km2, 10,000 km2, and 40,000 km2 squares significantly differed in 6.0% (mean D-statistic = 0.366), 3.2% (mean D-statistic = 0.357), and 2.3% (mean D-statistic = 0.346) of repetitions, respectively, and increases in geographic scale resulted in statistically significant increases in fit to the Point Aconi distribution (Supplementary Table G). KS test patterns are consistent across sample sizes for the Point Aconi site, but skewness and median body mass significantly differ from modern sites of increased scale at sample sizes of 40 and 60 species, while skewness also differs from 2500 km2 modern sites at the same sample sizes (Supplementary Table F2).
Box-and-whisker plots of D-statistics produced by Kolmogorov-Smirnov tests of goodness of fit between fossil body-mass distributions (Joggins in A, C, E, G; Point Aconi in B, D, F, H) and modern body-mass distributions produced by changing geographic scale of sampling. Modern distributions tested were subsampled at sample sizes equivalent to fossil distributions (A, B) and at increased sample sizes of 20 (C, D), 40 (E, F), and 60 species (G, H). Lower D-statistic values represent greater similarity between the simulated and fossil datasets. Red dots represent mean values for each array of D-statistics.
Figure 6. Long description
A multi-panel attribute space diagram consisting of eight box-and-whisker plots arranged in two columns and four rows. The left column (panels A, C, E, G) compares modern simulations to Joggins, and the right column (panels B, D, F, H) compares them to Point Aconi.
* Axes: The Y-axis represents the D-Value (D-statistic), where lower values indicate higher similarity. The X-axis represents the Scale of Sampling in kilometers-squared, with three categories: 2500, 10000, and 40000.
* Data Structure: Each plot contains three box-and-whisker elements. A red dot indicates the mean value, a horizontal black line indicates the median, and black dots represent outliers.
* Row 1 (A and B): Subsampled at fossil sample sizes. D-values are highest here, ranging mostly between 0.3 and 0.5. Means remain relatively stable across geographic scales.
* Row 2 (C and D): Increased sample size of 20 species. D-values show a slight downward trend as geographic scale increases, particularly for Point Aconi.
* Row 3 (E and F): Increased sample size of 40 species. D-values are lower than previous rows, with Point Aconi means dropping below 0.3 at larger scales.
* Row 4 (G and H): Increased sample size of 60 species. This row shows the lowest overall D-values. In panel H (Point Aconi), the mean D-value decreases significantly from the 2500 scale to the 10000 and 40000 scales, indicating that larger modern sampling areas more closely resemble the fossil distribution at higher species counts.
Controlling for diet improved goodness of fit between modern body-mass distributions and both sites (Fig. 7). Removing herbivores from the modern occurrences caused them to significantly differ from the Joggins body-mass distribution in 12.8% of repetitions (mean D-statistic = 0.411), while further restricting the modern occurrences to only include carnivores caused them to significantly differ from the Joggins body-mass distribution in 11.8% of repetitions (mean D-statistic = 0.427). KS tests consequently revealed that modern sites simulated to lack herbivores most significantly improved in goodness of fit to the Joggins distribution at three of the four tested sample sizes (Bonferroni-corrected p-values < 0.001). Removing herbivores from the modern dataset did not appear to improve fit between the modern body-mass distributions and those from Point Aconi at a sample size matching that of the fossil site, causing them to significantly differ in 11.2% of repetitions (mean D-statistic = 0.390). However, modern sites without herbivores show significant improvements of fit at higher sample sizes (Bonferroni-corrected p-values < 0.001), resulting in significant differences in 10.1% (mean D-statistic = 0.304) and 1.5% of repetitions (mean D-statistic = 0.298) at sample sizes of 40 and 60, respectively. Furthermore, restriction of modern occurrences to only include carnivores significantly increased fit with the Point Aconi assemblage compared with both the herbivore-removed and full modern datasets when compared at higher sample sizes (Bonferroni-corrected p-values < 0.001). After restricting modern occurrences to only carnivores, they only significantly differed from the Point Aconi assemblage in 8.3% (mean D-statistic = 0.316) and 0% of repetitions (mean D-statistic = 0.239) at sample sizes of 20 and 40, respectively.
Box-and-whisker plots of D-statistics produced by Kolmogorov-Smirnov tests of goodness of fit between fossil body-mass distributions (Joggins in A, C, E, G; Point Aconi in B, D, F, H) and simulated modern body-mass distributions with increasingly restricted dietary diversity. Modern distributions tested were subsampled at sample sizes equivalent to fossil distributions (A, B) and at increased sample sizes of 20 (C, D), 40 (E, F), and 60 species (G, H). Lower D-statistic values represent greater similarity between the simulated and fossil datasets. Red dots represent mean values for each array of D-statistics.
Figure 7. Long description
A multi-panel figure containing eight box-and-whisker plots arranged in two columns and four rows. The left column (A, C, E, G) compares modern simulations to Joggins, and the right column (B, D, F, H) compares them to Point Aconi.
All panels use a Y-axis labeled D-Value and an X-axis labeled Dietary Subset with categories: All Diets, No Herbivores, and Only Carnivores. Red dots indicate the mean D-Value for each group.
* Row 1 (A and B): Sample size equivalent to fossil distributions. In A, D-Values are relatively stable across subsets around 0.4. In B, D-Values are lower, centered near 0.35.
* Row 2 (C and D): Sample size of 20 species. In C, the Only Carnivores group shows a higher median and wider range than other subsets. In D, values remain consistent across subsets.
* Row 3 (E and F): Sample size of 40 species. In E, the Only Carnivores group shows a slight decrease in D-Value. In F, there is a sharp decline in D-Value from All Diets (near 0.65) to No Herbivores and Only Carnivores (below 0.3).
* Row 4 (G and H): Sample size of 60 species. These panels only show All Diets and No Herbivores. In G, No Herbivores shows a significantly lower D-Value (near 0.25) compared to All Diets (near 0.38). In H, a similar trend occurs with No Herbivores dropping to approximately 0.3 from a high of 0.6 in All Diets.
Controlling for clade and physiology in modern assemblages yielded variable similarities with fossil body-mass assemblages (Fig. 8). When compared with Joggins, modern mammalian body-mass distributions significantly differed in 17.8% of repetitions (mean D-statistic = 0.453), a significant reduction in goodness of fit with respect to the full modern distributions (Bonferroni-corrected p-value < 0.001). Modern reptilian assemblages also show significantly less similarity with the Joggins assemblage than the full assemblages (Bonferroni-corrected p-value <0.001), with body-mass distributions significantly differing from Joggins in 22.1% of repetitions (mean D-statistic = 0.478). Modern assemblages consisting only of sprawling, abducted-gait squamates and crocodylians most strongly differed from the Joggins body-mass distribution, with 54.1% of KS tests producing significant results (mean D-statistic = 0.588). These patterns are consistent at higher sample sizes. Conversely, the Point Aconi assemblage shows no significant differences in similarity between the full modern distributions and the modern reptile distributions (Bonferroni-corrected p-value = 1), with the latter significantly differing from Point Aconi in 8.3% of repetitions (mean D-statistic = 0.379). Compared with the full modern assemblages, both of the mammalian and abducted-gait modern assemblages are significantly less similar to the Point Aconi body-mass distribution (p-values < 0.001), with the mammals significantly differing from Point Aconi in 44.9% of repetitions (mean D-statistic = 0.499) and the abducted-gait squamates and crocodylians significantly differing from Point Aconi in 34.2% of repetitions (mean D-statistic = 0.534). As with the Joggins site, patterns for the Point Aconi site are also consistent at higher sample sizes, with the exception of the modern sites consisting of only reptiles, which show a significant improvement in fit from the full distributions at a sample size of 60 (p-value <0.001), but at no other sample size.
Box-and-whisker plots of D-statistics produced by Kolmogorov-Smirnov tests of goodness of fit between fossil body-mass distributions (Joggins in A, C, E, G; Point Aconi in B, D, F, H) and modern body-mass distributions simulated to include only mammals, only reptiles, or only tetrapods with abducted-limb body postures. Modern distributions tested were subsampled at sample sizes equivalent to fossil distributions (A, B) and at increased sample sizes of 20 (C, D), 40 (E, F), and 60 species (G, H). Lower D-statistic values represent greater similarity between the simulated and fossil datasets. Red dots represent mean values for each array of D-statistics.
Figure 8. Long description
A multi-panel grid of eight box-and-whisker plots labeled A through H. The y-axis for all panels is D-Value, representing the Kolmogorov-Smirnov test statistic. The x-axis is Physiological Subset. Panels in the left column (A, C, E, G) compare modern simulations to Joggins fossils. Panels in the right column (B, D, F, H) compare them to Point Aconi fossils.
* Row 1 (A and B): Sample sizes equivalent to fossil distributions. Panel A shows four categories: All Groups, Reptiles, Mammals, and Abducted-limb Tetrapods. D-values range from approximately 0.2 to 0.9, with Abducted-limb Tetrapods showing the highest mean (red dot) near 0.6. Panel B shows similar categories with lower overall D-values for All Groups and Reptiles.
* Row 2 (C and D): Sample size increased to 20. The distribution of D-values tightens. In both panels, Abducted-limb Tetrapods maintain the highest D-values, while All Groups and Reptiles show the lowest mean D-values, indicating higher similarity to the fossils.
* Row 3 (E and F): Sample size increased to 40. These panels exclude the Abducted-limb Tetrapods category. In panel E, All Groups has the lowest mean D-value (approx 0.38). In panel F, All Groups and Reptiles have similar low means near 0.3.
* Row 4 (G and H): Sample size increased to 60. These panels only show All Groups and Reptiles. In panel G, All Groups has a lower mean (approx 0.38) than Reptiles (approx 0.46). In panel H, both groups show very similar mean D-values near 0.28.
Across all panels, red dots indicate the mean D-statistic, horizontal black lines indicate the median, and black dots above the whiskers represent outliers.
Discussion
Body-Mass Estimation
The existence of a strong positive relationship between maximum OSA and body mass enabled complex comparisons of fossil assemblages to modern communities. While not as highly correlated as mean OSA with body mass, the absence of significant phylogenetic autocorrelation as well as absence of significant slope and intercept differences among subgroups (e.g., clade, diet) makes maximum OSA more applicable to estimation of the body masses of Pennsylvanian tetrapods. The lack of phylogenetic autocorrelation and differences among groups is curious, given that mean OSA, which incorporates the largest trunk vertebral centra, shows the opposite pattern. This difference implies that the smaller trunk vertebrae drive the variability in correlation between mean OSA and body mass among clades. In fact, minimum OSA shows the most variable correlation with body mass among clades. While there is no documented mechanism to explain the greater variability of the mean and minimum OSA–body mass correlations among groups, it may be indicative of variation among mammals and reptiles in the weight-bearing functions of the smallest and largest trunk vertebrae. Mammals and reptiles show significant differences in the primary orientation of flexion in vertebral joints (Jones et al. Reference Jones, Dickson, Angielczyk and Pierce2021). In many of the measured specimens, the largest trunk vertebrae were situated at the extreme ends of the trunk column, near the pelvic or pectoral girdles (see Supplementary Table A3). Joints in the pelvic region of the spine of both mammals and reptiles are often among the least flexuous joints in the trunk column, and the same is true for joints in the pectoral regions of some reptiles, (Molnar et al. Reference Molnar, Pierce and Hutchinson2014; Jones et al. Reference Jones, Dickson, Angielczyk and Pierce2021, Reference Jones, Angielczyk and Pierce2024). Thus, locomotor-related differences in functioning of vertebrae between mammals and reptiles should be weakest for the largest vertebrae, potentially explaining the absence of phylogenetic autocorrelation for maximum OSA. However, not all of measured specimens showed the greatest OSA values from the extreme ends of the trunk column, and the largest trunk centra of some of the measured mammals were found at the anterior end of the trunk column (Supplementary Table A3), which is expected to experience high flexion in mammals (Jones et al. Reference Jones, Angielczyk and Pierce2024). The potential effect of limb posture and spinal flexion on the mean OSA–body mass relationship cannot, therefore, be confidently extricated from any other physiological differences between mammals and reptiles.
Furthermore, the phylogenetic affinities of many Pennsylvanian tetrapods are either unresolved or the topic of ongoing debate (e.g., Modesto et al. Reference Modesto, Scott, MacDougall, Sues, Evans and Reisz2015; Ford and Benson Reference Ford and Benson2019; Adams et al. Reference Adams, Mann and Maddin2019; Mann and Paterson Reference Mann and Paterson2020; Mann et al. Reference Mann, Gee, Pardo, Marjanović, Adams, Calthorpe, Maddin and Anderson2020), and the trunk centra of both modern reptiles and modern mammals frequently show highly diverse, derived morphologies with irregular shapes (Rockwell et al. Reference Rockwell, Evans and Pheasant1938; Boszczyk et al. Reference Boszczyk, Boszczyk and Putz2001; Houssaye et al. Reference Houssaye, Mazurier, Herrel, Volpato, Tafforeau, Boistel and De Buffrénil2010; Čerňanský et al. Reference Čerňanský, Yaryhin, Ciceková, Werneburg, Hain and Klembara2018), especially in contrast to the generally circular centra of Pennsylvanian tetrapods (e.g., Romer and Price Reference Romer and Price1940; Carroll Reference Carroll1966; Reisz Reference Reisz1972; Klembara Reference Klembara1985). The clade-dependent differentiation in the body mass–mean OSA relationship means that the most appropriate use of the mean OSA estimation model would involve using separate models for reptiles and mammals. Given the uncertainty of phylogenetic affinity of the fossil taxa in this study, especially in the case of specimens with indeterminate taxonomy, we find that we lack sufficient evidence to definitively assign our fossil taxa to either of a mammalian- or reptilian-specific mean OSA–body mass estimation model, making the mean OSA model inappropriate for use in this particular case. It is for this reason that we chose to use maximum OSA for our body-mass estimates. Additionally, the reported equation should not be applied to the estimation of body mass in tetrapods suspected to be predominantly bipedal, which demonstrate substantial differences from quadrupeds in terms of axial skeletal morphology and the loading of weight upon the axial skeleton (Jiang et al. Reference Jiang, Moreau, Raso, Russell and Bagnall1995; Boszczyk et al. Reference Boszczyk, Boszczyk and Putz2001; Smit et al. Reference Smit2002; Balasubramanian et al. Reference Balasubramanian, Peters, Robinson, Singh and Kent2016; Russo et al. Reference Russo, Marsh and Foster2020).
Based on equation (3), the body-mass estimates we produced for Archaeothyris florensis and Hylonomus lyelli are lower than those of Modesto et al. (Reference Modesto, Scott, MacDougall, Sues, Evans and Reisz2015), while our body-mass estimate for Echinerpeton intermedium is higher. Due to the absence of direct comparison between body-mass estimations based on maximum OSA and Campione and Evans’s (Reference Campione and Evans2012) limb measurements, it is not clear which of the body-mass estimates for these taxa may best approximate biological reality. To ensure consistency of measurement through both fossil body-mass distributions, we opt to use the estimate provided by equation (3). Romer and Price’s (Reference Romer and Price1940) “pelycosaur” body-mass estimates are also generally higher than those produced by equation (3). Given that Dimetrodon limbatus was the model specimen for Romer and Price’s (Reference Romer and Price1940) estimation method, it is perhaps not surprising that their estimated body masses for large-bodied species of Dimetrodon and Edaphosaurus are the ones that best align with estimates from equation (3). In addition to our reappraisal, the body mass of Edaphosaurus boanerges was also estimated by Hurlburt (Reference Hurlburt1999) on the basis of graphic double integration and combined total length and tail girth. Hurlburt’s (Reference Hurlburt1999) best estimates of body mass from each method produced body masses of 92.79 kg and 75.05 kg, respectively, compared with the mean estimate of 79.10 kg produced by equation (3) and Romer and Price’s (Reference Romer and Price1940) mean estimate of 83 kg. Hurlburt’s estimate from total length falls within the equation (3) prediction error range for E. boanerges (86.52–97.75 kg) indicating that maximum OSA can provide a proxy for body mass that is repeatable with respect to graphic double integration methods based on total length. While Romer and Price’s (Reference Romer and Price1940) estimate of 83 kg falls outside of the predicted error range for E. boanerges, it is still reasonably close to the equation (3) estimate of 91.97 kg, indicating that the original estimates or Romer and Price are comparable to those of equation (3) in the specific case of large-bodied pelycosaurs.
Somewhat surprisingly, OSA–body mass relationships for fully terrestrial and semiaquatic animals were not significantly differentiated in our SMA results. Buoyant support of the trunk greatly reduces load stress on the tetrapod axial skeleton, resulting in substantial differences in vertebral morphology between fully aquatic and fully terrestrial tetrapods, particularly in terms of bone density (Buchholtz and Schur Reference Buchholtz and Schur2004; Long and Gordon Reference Long and Gordon2004; Galbusera and Bassani Reference Galbusera and Bassani2019). Intuitively, this should imply some differentiation of OSA–body mass relationships between fully terrestrial and semiaquatic tetrapods. However, it is possible that the requirement of supporting an organism’s full body weight while on land prevents vertebral characteristics of semiaquatic tetrapods from closely resembling those of fully aquatic tetrapods. The applicability of OSA–body mass relationships to both fully terrestrial and semiaquatic tetrapods allows us to include the embolomeres Calligenethlon watsoni and Carbonoherpeton carrolli, both possibly semiaquatic in life habit (Steen Reference Steen1934; Klembara Reference Klembara1985), in our estimates of fossil body-mass distributions.
Comparing Fossil and Modern Body-Mass Distributions
Comparison of body-mass distributions from fossil assemblages and modern communities yielded substantial apparent similarity, though the degree of similarity varies widely across modern sites (see Supplementary Information E). Similarity is somewhat enhanced when taking potential taphonomic size biases into account. For Joggins, modeling size bias produced results in line with expected taphonomic biases against small-bodied individuals (Behrensmeyer Reference Behrensmeyer1988, Reference Behrensmeyer, Allison and Briggs1991; Behrensmeyer et al. Reference Behrensmeyer, Kidwell and Gastaldo2000). However, no taphonomic bias significantly improved the fit of modern sites to the Point Aconi locality, which is odd given the expected action of stump preservation as a contributor to size bias. While the Joggins site includes some specimens found outside hollowed lycopsids, all of the Point Aconi fossils used herein were found within tree stumps (Carroll Reference Carroll1967b; Reisz Reference Reisz1972; Klembara Reference Klembara1985; Maddin et al. Reference Maddin, Mann and Hebert2019), which does not appear to have had a substantial impact on body-mass distributions. Absence of bias against large-bodied species at Joggins is somewhat supported by the ichnofossil record, where the largest gleno-acetabular distance (a close approximation of trunk length) measured is of a Hylopus hardingi trackway estimated at 129.8 mm (Marchetti et al. Reference Marchetti, Voigt, Buchwitz, MacDougall, Lucas, Fillmore, Stimson, King, Calder and Fröbisch2021), though ichnological evidence does exist of a potentially large-bodied trackmaker for which body lengths cannot be approximated (Prescott et al. Reference Prescott, Stimson, Dafoe, Gibling, MacRae, Calder and Hebert2014). In the body fossil record, skeletal fragments from Joggins putatively identified as cf. Baphetes, an aquatic or semiaquatic tetrapod with possible tetrapod or temnospondyl affinities (Milner et al. Reference Milner, Milner and Walsh2009; Sookias et al. Reference Sookias, Böhmer and Clack2014), may represent individuals larger than those with measurable vertebral elements, but definitive identification of large baphetids in the Joggins assemblage so far remains elusive (Hebert and Calder Reference Hebert and Calder2004; Grey and Finkel Reference Grey and Finkel2011; Sookias et al. Reference Sookias, Böhmer and Clack2014). Size filtering was shown to be the predominant source of taphonomic bias of the two sources of taphonomic bias that we accounted for when comparing to the Joggins site, as no increases of geographic scale precipitated consistent, significantly positive changes to goodness of fit between Joggins and the modern body-mass distributions. By contrast, the Point Aconi site did see significant improvement in K-S test results when compared with modern wetlands with broadened sampling areas. This is a surprising result, given assertions that the Joggins and Point Aconi tetrapod assemblages were primarily autochthonous or para-autochthonous in deposition (e.g., Falcon-Lang et al. Reference Falcon-Lang, Benton, Braddy and Davies2006). This may lend credence to Dawson’s (Reference Dawson1894) idea that some lycopsid stump preservation may have been a product of animals being washed into the hollowed stumps postmortem, instead of using the stumps as shelter (sensu Maddin et al. Reference Maddin, Mann and Hebert2019). However, multiple summary statistics showed significant deviation between the shapes of the Point Aconi distributions and those of all three geographic scales (see Supplementary Information F). Furthermore, less than 5% of modern sites at any geographic scale showed significant deviation from the Point Aconi distribution, making it seem unlikely that any effect of allochthony on body-mass distributions for the site is of substantial consequence. While the omission of some taxa as a result of taphonomic processes can be accounted for by taphonomic simulations, sample sizes are still likely to have a substantial impact on the overall shape of fossil body-mass distributions. The addition of future specimens may have significant effects on the shapes of the body-mass distributions, and so the body mass distributions produced for Point Aconi and Joggins should not be considered to be definitive representations of total community body-mass distribution.
Given that the taphonomic models specifically select against preservation of species that lie far outside of the body-mass ranges observed in the fossil sites, it should be expected that some aspects of modern body-mass distributions have been made artificially similar to the fossil distributions. Of the statistics covered herein, this would particularly affect median body masses, and comparison of median body-mass values between fossil and taphonomically controlled modern sites should therefore be considered highly likely to be subject to overfitting. However, we do not expect that the overall fit of fossil to taphonomically controlled modern body-mass distributions is a product of overfitting. If this were the case, both the Joggins and Point Aconi assemblages should clearly best fit the modern distributions, which were simulated to be taphonomically biased against both very large and very small species. This is not the case for either assemblage, as Joggins is best fit to taphonomic bias against only very small species, and Point Aconi fits as well to taphonomically uncontrolled sites as it does to taphonomic bias against extremely large species. Neither site showed strong fit with modern sites subsampled to simulate bias against both very large and very small species. These results therefore indicate that these KS-test results are not a result of overfitting in the taphonomic models.
Taking into account the results of the taphonomic models, we could then identify if any further ecological factors substantially influenced differences between the fossil and modern distributions. The major significant improvement in fit between modern sites and either fossil assemblage arose from modeling body-mass distributions of modern communities with increasingly restricted dietary diversity. The significant improvement of fit implies that some degree of difference in body masses might be attributable to low dietary diversity in both the Joggins and Point Aconi assemblages. Of the tetrapods known from Point Aconi and Florence, two (A. florensis and E. intermedium) have been directly interpreted as insectivorous or generally carnivorous (Reisz Reference Reisz1972), while the most evidence of early herbivory in tetrapods comes from younger fossils than those found at Point Aconi (Sues and Reisz Reference Sues and Reisz1998; Brocklehurst et al. Reference Brocklehurst, Kammerer and Benson2020). However, the recent description of the pantylid Tyrannoroter heberti by Mann et al. (Reference Mann, Xiong, Calthorpe, Sues and Maddin2026) suggests that there may have been herbivorous tetrapods in the Point Aconi fossil community, though this remains uncertain (Whittingham Reference Whittingham2025). Tyrannoroter heberti is known only from skull material, and so is not included in our body-mass distribution for Point Aconi. Had the nine measured specimens from Point Aconi all been carnivores, it might have impacted the resulting body-mass distributions (Price and Hopkins Reference Price and Hopkins2015; Pineda-Munoz et al. Reference Pineda-Munoz, Evans and Alroy2016). The high functional diversity of modern omnivore species (Reuter et al. Reference Reuter, Hopkins and Price2023), in particular, may explain the increase in goodness of fit in the exclusively carnivorous subset, as absence of omnivorous species may have substantially influenced fossil body-mass distributions. It is important to note that removing herbivores and omnivores from the modern tetrapod body-mass distributions is not the same as modeling environments with only carnivores, but rather modeling environments where other dietary niches are filled, but not represented in body-mass distributions. Therefore, the improved fit of these subsets to the Point Aconi distribution does not imply that the Point Aconi assemblage lacked herbivores or omnivores, but rather that those roles may not have been taken up by the nine measured species.
Carboniferous coal swamp ecosystems possessed arthropod fauna with a substantial diversity of diets and body masses (Scott and Taylor Reference Scott and Taylor1983; DiMichele and Hook Reference DiMichele, Hook, Behrensmeyer, Damuth, DiMichele, Potts, Sues and Wing1992; Greb et al. Reference Greb, DiMichele, Gastaldo, Greb and DiMichele2006); the Sydney Mines Formation, which hosts the Point Aconi assemblage, boasts evidence of such an arthropod fauna (e.g., Dunlop Reference Dunlop2018), including trackways attributable to the large millipede-like Arthropleura (Falcon-Lang Reference Falcon-Lang2006). Even more so than Point Aconi, the Joggins assemblage hosts diverse arthropod trace and body fossils (Copeland Reference Copeland1957; Calder Reference Calder, Blundell and Scott1998; Falcon-Lang et al. Reference Falcon-Lang, Rygel, Calder and Gibling2004, Reference Falcon-Lang2006; Grey and Finkel Reference Grey and Finkel2011; Buhler and Grey Reference Buhler and Grey2017) and similarly shows improved fit to modern distributions with restricted dietary diversity, though the Joggins site is best fit to modern sites with herbivores removed, instead of those restricted to only carnivorous tetrapods. Therefore, it is plausible that herbivorous and omnivorous niches now occupied by tetrapods in modern assemblages were instead mostly occupied by arthropods in these communities. It is common for coastal wetland assemblages to have their food webs incorporate majority arthropod or aquatic primary consumers (Polis and Hurd Reference Polis and Hurd1996; Douglas et al. Reference Douglas, Bunn and Davies2005; Abrantes and Sheaves Reference Abrantes and Sheaves2008; Abrantes et al. Reference Abrantes, Johnston, Connolly and Sheaves2015), but the results of our modeling appear to indicate that tetrapods may still have played less diverse of a dietary role than expected in Pennsylvanian communities. These results imply that the tetrapod communities of Joggins and Point Aconi were in line with terrestrial communities expected to be typical of the Pennsylvanian, with terrestrial herbivory restricted to invertebrates and a combination of marine and terrestrial nutrients being filtered up to tetrapod faunivores through piscivory and insectivory (Olson Reference Olson1966; Otoo Reference Otoo2025).
The phylogenetic placement of many of the Joggins tetrapods on the amniote stem (e.g., Mann et al. Reference Mann, Gee, Pardo, Marjanović, Adams, Calthorpe, Maddin and Anderson2020; Arbez et al. Reference Arbez, Atkins and Maddin2022) may explain why we did not find improved goodness of fit among body-mass distributions using various taxonomic subsets of the modern data. Of the species included in the body-mass distribution from Joggins, only H. lyelli is confidently identified as a crown amniote (Carroll Reference Carroll1964), while Protoclepsydrops haplous and Asaphestera intermedia have both been dubiously associated with Pelycosauria (Reisz Reference Reisz and Westphal1986; Mann et al. Reference Mann, Gee, Pardo, Marjanović, Adams, Calthorpe, Maddin and Anderson2020). Body-mass distributions are not likely to be influenced by the squat, abducted-limb Bauplans of the studied fossil tetrapods, given that restriction of modern body-mass distributions to reptiles of similar shapes yielded the least similarity to either fossil assemblage of all biologically defined modern subsets. This last result also challenges the idea that such Bauplans are a substantial source of body-size constraint (sensu Buchwitz et al. Reference Buchwitz, Jansen, Renaudie, Marchetti and Voigt2021). While taxonomic affinity and Bauplan did not improve goodness of fit between the fossil and modern distributions, it is possible that applying both dietary and taxonomic or body-shape restrictions to the modern datasets may do so. Unfortunately, stacking multiple subsampling methods for the modern sites results in there being too few sites with sufficient diversity for adequate comparison. Given the significant negative impacts of taxonomic and Bauplan restrictions on goodness of fit between modern and fossil sites, it is unlikely that combining those restrictions with dietary restrictions would have improved fit, though this expectation cannot be quantitatively verified with the data available.
Conclusions
The results presented herein demonstrate similarity between the body-mass distributions of modern and Pennsylvanian tetrapods in tropical deltaic and estuarine settings. While size-bias simulations also showed that common taphonomic biases against small organisms likely impact the array of species preserved at Joggins, taphonomy alone appears unable to explain the differentiation in body-mass distributions. While few modern and fossil body-mass distributions could be significantly distinguished between ecologically unaltered distributions, emulating restrictions in dietary diversity significantly improved the fit of modern and fossil distributions. Our results therefore reflect differences in tetrapod functional diversity between modern and Pennsylvanian wetland communities that likely stem from the limited dietary diversity of Pennsylvanian tetrapods. Finally, the predictive utility of the maximum OSA–body mass relationship in extant quadrupedal tetrapods arms paleontologists with an improvement on a preexisting estimator of body mass based on commonly preserved skeletal elements, making it easier for future studies to help improve our understanding of the evolution of terrestrial ecosystems through the lens of body mass.
Acknowledgments
The authors would like to thank J. Atkins, H. Larsson, and T. Fedak of the Joggins Fossil Institute, Redpath Museum, and Nova Scotia Museum, respectively, for allowing access to and loans of fossil materials for this study. We would also like to thank D. Evans and S. Sugimoto of the Royal Ontario Museum for providing access to modern skeletal materials and our reviewers for their constructive criticism. Special thanks to N. Brocklehurst, who provided his reverse-engineered estimation of Romer and Price’s (Reference Romer and Price1940) original body mass–DVC relationship. The Natural Sciences and Engineering Research Council of Canada funds D.F. (NSERC RGPIN-2018-05305) and H.C.M. (NSERC RGPIN-2020-06442). We also acknowledge that many of the fossils used in this study were collected on the traditional territory of the Mi’kmaq and Maliseet First Nations without their consent or knowledge as a product of centuries of colonial repression and exclusion.
Competing Interests
The authors declare no competing interests.
Data Availability Statement
Data available from the Dryad and Zenodo Digital Repositories: https://doi.org/10.5061/dryad.tdz08kqg0 and https://doi.org/10.5281/zenodo.22262672.
Open access
