Fecundity in fossil Bryozoa: accounting for colony fragmentation and the spatial division of reproductive labor

Abstract

Our ability to measure evolution by natural selection in the fossil record is limited by the near impossibility of estimating the fecundity and thus relative fitness of most fossil organisms. Neocheilostome bryozoans are an important exception, because they have calcified larval brood chambers known as ovicells that provide an approximate estimate of the colony’s sexual fecundity. This clade has a rich fossil record dating back ~100 million years, providing potential opportunities to observe changes in relative fitness and natural selection through many past intervals of environmental change. However, neocheilostome fossil specimens are often highly fragmented, and fragments are not necessarily randomized subsets of a colony. To make use of the majority of the neocheilostome fossil record, we need to test the effect colony organization has on our methods of inferring colony fecundity from fragmented specimens.

In this study, we measure colony fecundity in a population of Recent neocheilostome bryozoan specimens of the species Parasmittina eccentrica Winston & Jackson, 2021 and quantify the nonrandom spatial arrangement of ovicells due to colony organization. We then simulate fragmenting these specimens and test the statistical robustness of standard methods one might use to reconstruct fecundity from fossil specimens. We find that ovicells are clustered and concentrated at mid-distances from the ancestrula (the oldest part of the colony). As a result, estimates of a colony’s fecundity from a single fragment have higher variance than would be expected if ovicells were randomly distributed. When estimating average population fecundity, observed variance among fossil fragments is a better estimator of sample variance than methods that assume spatial independence (such as a binomial distribution), especially for fragment sizes of 8 mm or less. While there is much to be learned about neocheilostome ovicell arrangement across taxa and environments, we can robustly estimate fecundity from small fossil fragments even in extinct neocheilostome species.

Non-technical Summary

Evolution by natural selection occurs when individuals in a population have differential reproductive success due to heritable traits. Reproductive success, or fecundity, is near impossible to estimate in most fossil organisms, so paleobiologists’ ability to conclusively demonstrate natural selection in fossil time series is limited. This prevents us from understanding the role natural selection plays in processes only observable over long timescales, such as major climate changes and mass extinctions. There is at least one exceptional study system, however; the neocheilostome bryozoans are a monophyletic clade of colonial, marine invertebrate suspension feeders. They brood their larvae in calcified (skeletonized) chambers known as ovicells, which provide an approximate estimate of sexual fecundity even in fossil specimens. The neocheilostome fossil record goes back ~100 million years, and they are abundant, globally distributed, speciose, and extant. These many advantages make them a rising model system for measuring fecundity, relative fitness, and natural selection in the fossil record. However, most fossil specimens are highly fragmented, creating the need for careful methods of inferring colony fecundity from fossil fragments.

In this study, we take a population of Recent neocheilostome bryozoan specimens of the species Parasmittina eccentrica and map the spatial arrangement of ovicells across these colonies. We quantify the nonrandom arrangement of ovicells, which is due to the integration of individuals in a colonial organism. We then simulate fragmenting these specimens, and we test the efficacy of a variety of statistical methods one might use to estimate the original colony fecundity from fragmented fossil specimens. We find that in the studied population, ovicells are clustered and concentrated at mid-distances from the ancestrula (the oldest part of the colony). As a result, estimates of a colony’s fecundity from a single fragment have higher error than expected. Furthermore, how we estimate average population fecundity gives variable results, especially for small fragment sizes.

Introduction

Fecundity is fundamental to measuring relative fitness and thus evolution by natural selection, and it is almost impossible to estimate in the fossil record. The result is a disconnect, wherein neontologists measure fitness but can only extrapolate the potential macroevolutionary consequences, while paleontologists observe macroevolutionary patterns but can only interpolate the possible microevolutionary mechanisms (Reznick and Ricklefs 2009; Kearney et al. 2024). Rare opportunities to measure fecundity in the fossil record allow us to connect relative fitness, the strength of natural selection, the rate of trait evolution, and the rates of speciation and extinction, testing the mechanistic assumptions made by both fields (O’Dea and Jackson 2009; Liow et al. 2024). Furthermore, reproductive investment is itself an evolving trait, and the ability to measure fecundity in the fossil record can help us understand how the diversity of reproductive strategies in the modern biosphere arose. The fossil record of neocheilostome bryozoans offers one such opportunity.

The Order Cheilostomatida within Phylum Bryozoa is a monophyletic clade of colonial marine invertebrate suspension feeders (Waeschenbach et al. 2012). The cheilostome fossil record begins in the Late Jurassic, and cheilostomes became the most abundant and speciose bryozoan order in the Late Cretaceous (Jablonski et al. 1997; Lidgard et al. 2012, 2021; Sepkoski et al. 2000). With more than 600 extinct and extant genera (Lidgard et al. 2021), it is an ecologically broad clade with diverse functional traits that are preserved in their rich fossil record, making them a rising model system for studies bridging micro- and macroevolution (Di Martino and Liow 2021; Liow et al. 2024). Cheilostome bryozoan colonies are composed of individuals known as zooids. Zooids can reproduce sexually as simultaneous hermaphrodites, sequential hermaphrodites, or gonochronistically, or they can reproduce by asexual budding (Ostrovsky 2013). Colonies can thus reproduce asexually by accidental fragmentation or self-induced fragmentation (O’Dea et al. 2008). The neocheilostomes, the dominant monophyletic clade in the order, have synapomorphic larval brood chambers, which are usually a calcified, external structure developed by the zooid distal to the zooid that is producing the zygote; these are ovicells sensu stricto (Grant et al. 2023). Reproduction in neocheilostomes thus requires colony integration, so that fecundity and reproductive strategy are colony-level as well as zooid-level traits. Ovicell density is therefore theoretically a proxy for colony investment in sexual reproduction (Di Martino and Liow 2021; Liow et al. 2024), and ovicell arrangement on a colony is a potential proxy for the strategic timing of sexual reproduction. Because ovicells are calcified, a record of sexual reproduction can be preserved in neocheilostome specimens, whether they be living, deceased, or fossil.

A small number of studies have begun the use of ovicell density on fossil specimens as a proxy for colony sexual fecundity. In cheilostome bryozoans from the Cretaceous/Paleogene boundary of Denmark, Håkansson and Thomsen found a trade-off between a colony’s investment in sexual and asexual reproduction; colonies with more fragile (arborescent) growth forms had lower ovicell densities than more robust colonies, a pattern consistent across all species studied over a 3 million year time series (Thomsen and Håkansson 1995; Håkansson and Thomsen 2001). Di Martino and Liow (2021) and Liow et al. (2024) advanced this work by using ovicell density to estimate trait-fitness associations in Antarctothoa and Microporella from the Pleistocene–Recent fossil record of the Wanganui Basin, New Zealand. However, our ability to further expand the use of this proxy is hindered by the fact that many neocheilostomes appear to have an organized arrangement to their ovicells, and fossil specimens are often highly fragmented. In species where ovicell locations are not random, the smaller the fragment, the greater the error in our estimate of fecundity, and the less able we are to discern trends in the data. We can leverage living neocheilostomes to better estimate fecundity from fragmented fossil specimens with unknown degrees of nonrandom ovicell arrangement.

Observations and experiments on living neocheilostomes have found ovicell density and ovicell arrangement can vary throughout colony development and are influenced by environment. Isolated case studies have found a colony’s reproductive investment and timing can depend on feeding currents (Chimonides and Cook 1981; Hunter and Hughes 1993; Shunatova and Ostrovsky 2002), zooid generation (Harvell 1984), substrate (Yagunova and Ostrovsky 2010), season (Lombardi et al. 2006), symbionts (Mathew et al. 2016), or stress such as ocean acidification (Swezey et al. 2017). Nonrandom distribution of ovicells in colonies is also sometimes qualitatively noted (O’Dea et al. 2010). Arborescent colony growth forms in particular are often constrained as to which zooids can develop ovicells due to the three-dimensional colony structure (Fairall 2004). Thus, we know anecdotally that ovicell arrangement is nonrandom, and a fragment of a colony is likely not a representative subsample of the entire colony. Yet to our knowledge no quantitative models of ovicell arrangement exist.

Here we examine a population in the neocheilostome genus Parasmittina from a previously published study (Jackson and Cheetham 1990), and we map the presence/absence of ovicells across 53,166 zooids belonging to 27 colonies. With an image dataset that uniquely combines high magnification and coverage of whole colonies, we contribute the most data-intensive quantitative analysis of ovicell density and arrangement in a bryozoan population to date. We then evaluate the statistical robustness of different methods of estimating fecundity from simulated colony fragments and suggest a most robust protocol.

Methods

Study System

The collection of specimens used in this study are from the experiments of Jackson and Cheetham (1990). At the time of their study, these specimens were identified as Parasmittina areolata Canu & Bassler, 1927. The species taxonomy within this genus was revised by Winston and Jackson (2021) and then revised again by Farias et al. (2024). We reassign the specimens of this collection to Parasmittina eccentrica Winston & Jackson, 2021, based on their Caribbean origin, 3 avicularia types, 20+ ovicell pores, 2–3 oral spines, and size measurements (see Winston and Jackson 2021; Farias et al. 2024; Supplementary Fig. 1). The original study was conducted at the Smithsonian San Blas Field Station (Jackson and Cheetham 1990). The parent generation was collected on coral fragments on the Caribbean coast of Panama. Non-bryozoan fouling organisms were removed from the coral fragments, and then colonies with visible embryos were each isolated in tanks with running seawater that had been filtered for larvae. Bare coral fragments were left with each parent colony for 5–10 days. Then coral fragments bearing the F₁ generation were set back in the sea attached to cement blocks and cleaned of non-bryozoan fouling organisms each month. Thus, the parent generation may have experienced some competition before collection, but competition was eliminated for the F₁ generation. The F₁ colonies have a known maternal colony but an unknown paternal colony. While the exact age of each colony is not known, the F₁ generation were all allowed to reach sexual maturity. After the original study, specimens were desiccated and deposited in the collections of the Smithsonian National Museum of Natural History (Collection USNM PAL 787), where they remain available. The colony IDs used in this study are consistent with the Jackson and Cheetham (1990) and USNM specimen labels; “Pa##” indicates the parent colony, and “M###” indicates the offspring number, so that Pa55 is the parent of Pa55M11.

Data Collection

Scanning electron microscopy (SEM) has been the traditional method for bryozoan taxonomy, but the specimen size limitations of most SEM setups have prohibited high-magnification images of entire large bryozoan colonies. The newest SEM technology is capable of large-area image stitching, but these facilities are still prohibitively rare and expensive to access. This is likely in part why no studies to our knowledge have quantified the arrangement of polymorphic zooids across entire colonies. To address this gap, we image 27 P. eccentrica colonies with a Keyence VHX-7000 digital light microscope. The Keyence takes both focal stack and XY-stitched images, allowing us to take composite images of entire colonies (Fig. 1A).

Figure 1. Illustration of how colonies were mapped. A, Unlabeled image of colony Pa55 (green dye added by Jackson and Cheetham [1990]), B, Red arrows indicate zooids with ovicells, and blue arrows indicate zooids without ovicells. Ovicells are visible as a raised crescent of pores distal to the orifice. C,D, Examples of colony regions containing zooids with and without ovicells. E, Fully labeled colony.

We used the free software 3D Slicer v. 5.8.0 for image annotation (Fedorov et al. 2012; Kikinis et al. 2014). For each colony image, we placed binary point labels on zooids with ovicells or without ovicells (Fig. 2). We place a third point class on the ancestrula of each colony (the oldest part of the colony, where the larvae metamorphosed into the three original zooids) in colonies where the ancestrula is visible. In some of the parent generation, the ancestrula is not contained in the portion of the colony that was collected from the wild. In other cases, the ancestrula is covered by secondary frontal budding. Zooids that could not be classified, whether due to limitations of our imaging methods or colony abnormalities, were left unlabeled. In total, we classified 53,166 zooids from 5 parent colonies and 22 offspring colonies.

Figure 2. Ovicell maps of 27 Parasmittina eccentrica colonies. Blue points indicate zooids without ovicells. Red points indicate zooids with ovicells. Blank areas within a colony represent unclassified zooids. The axes are centered on the approximate location of the ancestrula. In some colonies, the direction of the ancestrula lies outside the specimen, so the location of the ancestrula was extrapolated.

In the following analyses, we refer to the ovicell density of the entire observed colonies as the “true” fecundity (“true” in quotation marks, because there are some omitted zooids and some specimens with broken margins). Even for the full colonies, ovicell density is still an estimate of the true value. Compared with the fossil fragments we simulate, however, these near-entire colonies are a good approximation of the true ovicell density.

All analyses are conducted with R v. 4.4.2 in RStudio v. 2024.12.0+467. Slicer data are imported into R with the Morpho v. 2.12 and geomorph v. 4.0.9 libraries (Schlager 2017; Baken et al. 2021; Adams et al. 2024). Simulations are parallelized with R libraries doParallel v. 1.0.17 and foreach v. 1.5.2 (Microsoft Corporation and Weston 2022a,b). Figures are made with R libraries tidyverse v. 2.0.0, ggplot2 v. 3.5.1, ggtext v. 0.1.2, ggforce v. 0.4.2, ggpubr v. 0.6.0, and scales v. 1.3.0 (Wickham 2016; Wickham et al. 2019, 2023; Wilke and Wiernik 2022; Kassambara 2023; Pedersen 2024).

Analyses: Quantifying Ovicell Arrangement

We measure fecundity in terms of ovicell density, the proportion of labeled zooids with ovicells in a certain area. If ovicell arrangement is nonrandom, we expect ovicell density within a given area to have higher variance than if ovicell arrangement were random. Methods exist to test whether a spatial point pattern is random, uniform, clustered, or dispersed (Perry et al. 2006). However, these methods are applicable to populations of organisms that can occur at any distance from one another. In the case of a bryozoan colony, zooids are of a limited size range, so ovicells can only occur fixed distances apart. Furthermore, each observed colony has a different shape and size and different arrangements of unlabeled zooids. Therefore, rather than applying an existing method to test whether the distribution of ovicells among zooids in a colony is nonrandom, we simulate a null distribution based on the real colony geometries.

Our null hypothesis is that ovicells are randomly arranged among zooids in a colony. To generate this null distribution, we take each colony’s number of observed ovicells and randomly assign those ovicells among the colony’s labeled zooids. Unclassified zooids are still omitted. We repeat this simulation many times for each colony (the exact number of iterations depends on the analysis, described below) to generate the null distribution for each analysis. We then compare the observed data with each null distribution. Where the observed data fall outside the 95% quantiles of the null distribution, we conclude ovicell arrangement is nonrandom. We apply this approach to the following two analyses.

First, we test the hypothesis that ovicells are clustered; if so, as distance from an ovicell increases, the density of other ovicells should decrease. To quantify this, we measure the density of other ovicells as the distance from each ovicell increases in 2 mm increments. We repeat the process for our randomized null dataset to generate a null distribution. We then compare the observed correlation between ovicell density and distance with the null distribution to test whether ovicells are significantly clustered.

Second, we test the hypothesis that ovicell density varies systematically with distance from the ancestrula. In the subset of colonies in the parental and F₁ generations with an identifiable ancestrula (16 colonies), we measure ovicell density as the distance from the ancestrula increases in 2 mm increments. We repeat this process for our randomized null dataset to generate a null distribution. We then compare the observed pattern with the null distribution to test whether ovicell density varies significantly with distance from the ancestrula in each colony. We also fit a generalized additive model (GAM) to estimate how much ovicell density changes with each 2 mm increment from the ancestrula. A GAM is appropriate, because we do not expect the relationship to necessarily be linear or monotonic, and we want to model the relationship without hypothesizing any specific polynomial relationship (Fieberg 2022; R library gam v. 1.22-5 [Hastie 2024]).

Analyses: Simulating Sampling the Fossil Record

Our ultimate goal is to understand how nonrandom ovicell arrangement affects our estimate of colony fecundity (i.e., ovicell density) in fossil neocheilostome specimens of variable size. To simulate specimen fragmentation, we take the (X,Y) range of each colony (with units in mm) and randomly sample an X coordinate and a Y coordinate from a uniform distribution 500% wider along each axis than the colony itself. Sampling from a uniform distribution wider than the colony prevents the center of the colony from being oversampled (an edge effect). We then sample the zooids that fall within a radius r of the (X,Y) point, where r is drawn from a uniform distribution between 2 and 10 mm. This generates one simulated fossil fragment. We continue to resample the colony without replacement until the entire colony is sampled. Each zooid can only be sampled once. If the next random circle overlaps with previously sampled parts of the colony, the resulting fragment will be noncircular in shape. Therefore, the random fragments have a variety of shapes. This process is one simulation, which we repeat 500 times.

We then apply standard statistical approaches that a paleobiologist might use to (1) estimate colony fecundity from a fragment and (2) estimate the difference in average fecundity between samples of fragments. For each of these tests, we compare the mean and standard error of the estimates to the “true” value calculated from the complete colonies. If the estimate and its 95% confidence interval do not contain the “true” value 95% of the time, then the method is either biased or misestimating error. We also want to know whether the failure of a method is due to ovicell clustering, so we need a null distribution. We follow the same procedure to simulate fragments of colonies with random ovicell locations. We repeat each test for these randomized fragments to generate a null distribution for each analysis.

Analyses: Estimating Colony Fecundity from a Fragment

If ovicell arrangement is random, then each zooid has an equal and independent probability of bearing an ovicell. The number of ovicells on a fragment will therefore be a binomial random variable. The ovicell density of the fragment, $ \hat{p} $ , is our estimate of the ovicell density of the entire colony. The normal approximation of binomial standard error is $ \sqrt{\left(\hat{p}\left(1-\hat{p}\right)/n\right)} $ , where n is the number of zooids in the fragment. The 95% confidence interval is $ \hat{p} $ $ \pm $ 1.96 * $ \sqrt{\left(\hat{p}\left(1-\hat{p}\right)/n\right)} $ . By definition, the 95% confidence interval will capture the “true” value in 95% of iterations. If ovicell arrangement is nonrandom, then we expect the “true” value will fall outside the 95% confidence interval for more than 5% of iterations. That is, the observed error will exceed the standard error of a binomial random variable.

Analyses: Estimating the Average Fecundity of a Population of Fragments

Although the previous analysis is useful to understand how nonrandom ovicell arrangement may inflate error, in a realistic study, a paleobiologist would sample a population of fragments and estimate an average fecundity. While in practice one would not know which or how many colonies a sample of fragments come from, for simplicity, we treat the fragments of each colony as a population. For each simulated iteration of colony fragmentation, we subsample 30 fragments and estimate the mean ovicell density and its standard error. We repeat this simulation for fragments of 2–4 mm, 2–6 mm, 2–8 mm, and 2–10 mm, because fragment size will affect the total fraction of each colony represented in 30 fragments. We refer to these as “size classes.” The size classes are overlapping, with only the maximum fragment size differing between classes. We repeat this for our randomized colonies to produce null distributions. A 2-mm-radius fragment is the smallest colony fragment with a reasonable sample size of zooids with which to measure ovicell density. A 10-mm-radius fragment is almost the entire colony. Therefore, we are representing the range of preservation quality one might find in the fossil record.

There are multiple methods one might use to estimate the mean and standard error of each population, so we repeat this analysis for three different methods and compare the statistical robustness of each method. Method 1 takes the mean ovicell density of sampled fragments and estimates standard error from a binomial distribution. This is the same procedure as using an unweighted binomial generalized linear model (GLM). Method 2 takes the proportion of zooids with ovicells pooled across all fragments (a weighted mean) and estimates standard error from a binomial distribution, akin to a binomial GLM weighted by the number of zooids in each fragment (Di Martino and Liow 2021). Method 3 takes the proportion of zooids with ovicells pooled across all fragments (a weighted mean) and estimates standard error from the observed error among fragments. Exact equations for each method are shown in Table 1. We compare the efficacy of the three methods for the four fragment size classes. A robust method will have a 95% confidence interval that captures the “true” value in 95% of iterations.

Table 1. The three methods used to estimate mean population ovicell density and its standard error (SE)

Results

Ovicell Arrangement Is Nonrandom

Across all colonies, mean ovicell density is higher near other ovicells than expected at random (Fig. 3). In other words, ovicells are clustered. Clustering is greatest within 2 mm of an ovicell and decreases with distance. Zooids in this population are 400–800 microns in length from orifice to distal orifice and 300–500 microns at their widest, so a 2 mm cluster is approximately 4–6 zooids in diameter.

Figure 3. A, There is a negative correlation between ovicell density and distance from another ovicell in the observed colonies (pink) but not the null distribution (gray). Therefore, in every colony, the density of ovicells is higher around another ovicell than expected at random. Clustering is most significant within 2 mm of an ovicell. B, To visualize the strength of this correlation averaged across colonies, ovicell density is detrended by subtracting the mean of the null distribution for each colony.

Across all colonies, ovicell density varies with distance from the ancestrula more than expected at random (Fig. 4). Ovicells are less dense within 8 mm of the ancestrula, and more dense within 10–18 mm of the ancestrula. For the subset of colonies that reached greater than 20 mm, ovicell density also decreases in the 24–26 mm interval.

Figure 4. A, Ovicell density in the observed colonies (pink) varies with distance to the ancestrula more than expected at random (gray). B, To visualize the average change in ovicell density with distance from the ancestrula across all colonies, ovicell density is detrended by subtracting the mean of the null distribution for each colony. Colored lines indicate where the individual colonies contribute to each box plot. C, Generalized additive model (GAM) fit to detrended data shows ovicell density is significantly lower within 8 mm of the ancestrula and significantly higher 10–18 mm from the ancestrula. The gray ribbon is the 95% confidence interval.

Estimated Colony Fecundity from Fossil Fragments Has High Error

Because ovicell arrangement is nonrandom, error in estimating colony ovicell density from a single fragment exceeds acceptable limits (Fig. 5). We define acceptable error as within $ \pm $ 1.96 standard errors of the “true” value. Even for fragments containing at least 100 zooids, a seemingly reasonable sample size, 71% produce unacceptable error. The median error is $ \pm $ 3.6 standard errors, and the 95^th percentile is $ \pm $ 15.5 standard errors from the “true” value. Real fossil fragments (and SEM images) are often less than 100 zooids. When ovicell arrangement is randomized, error does fall within expected limits (yellow band in Fig. 5). The question remains whether, with a sufficient sample size, fragments can still be used to estimate an unbiased mean.

Figure 5. A, Ovicell density of each fragment (points) compared with the “true” ovicell density of the entire colony (black line). B, The error in estimating colony ovicell density from a single fragment versus the number of zooids in the fragment. Error is measured in terms of the number of binomial standard errors the estimate is from the “true” value. The acceptable error limit is $ \pm $ 1.96 standard errors, a 95% confidence interval, highlighted in red. The null distribution, when ovicell arrangement is randomized, is highlighted in yellow. The black points indicate fragments of the observed colonies.

Method 3 Is the Most Robust Estimate of Average Population Fecundity across Fragment Sizes

We use each of the methods in Table 1 to estimate the mean and standard error of ovicell density for a random subsample of colony fragments of four nested size classes, 2–4 mm, 2–6 mm, 2–8 mm, and 2–10 mm. We then see in what proportion of iterations the 95% confidence interval captures the “true” value; if less than 0.95, the method is underestimating sample variance and has too high a Type I error rate, and if greater than 0.95, the method is overestimating sample variance and has too high a Type II error rate. For all three methods and all four size classes, we find the 95% confidence intervals capture the “true” value 100% of the time for the randomized colonies comprising the null distributions.

For the observed colonies, methods 1 and 3 are robust; the 95% confidence interval of the estimate captures the “true” value at least 95% of iterations for all four size classes. Method 2, however, is robust for the 2–10 mm size class, but has a success rate below 95% for fragment sizes of 2–4 mm, 2–6 mm, and 2–8 mm (Fig. 6). Method 1, though robust, performs worse when fragment sizes are more variable (2–10 mm), because the unweighted average ovicell density has higher variance than estimated for a binomial random variable. Method 2 performs worse when fragment sizes are smaller and less of the colony is represented in the sample, because sample variance is again underestimated by a binomial distribution. Method 3 is the most consistent and robust method, because the observed variance among fragments is the best estimate of the real sample variance. In other words, the observed variance among fragments is a better estimator of sample variance than a binomial distribution, because the assumptions of a binomial random variable are not met. The estimates and their 95% confidence intervals for all three methods and 500 iterations per colony are illustrated for the 2–4 mm size class in Supplementary Figure 2.

Figure 6. Robustness of three standard statistical methods for estimating average population fecundity from a random sample of colony fragments. The x-axis is the maximum fragment size (in mm), representing size classes 2–4 mm, 2–6 mm, 2–8 mm, and 2–10 mm. The y-axis is the proportion of accurate estimates, meaning the proportion of the 500 iterations per colony (500 × 27 colonies = 13,500 iterations per point) in which the 95% confidence interval captures the “true” value. A robust method will have a proportion >0.95, and this cutoff is indicated with a bolded line. Methods 1, 2, and 3 correspond to the methods described in Table 1. “Observed” refers to the observed colonies with their real ovicell locations, and “Null” refers to the simulated colonies with randomized ovicell locations.

Discussion

Neocheilostome Bryozoa are a promising model system due to their rich fossil record and calcified larval brood chambers, which provide the opportunity to estimate fecundity from fossil specimens (Di Martino and Liow 2021; Liow et al. 2024). The opportunity to measure fecundity in the fossil record raises the potential to measure natural selection during past episodes of rapid environmental change and to observe the evolution of reproductive strategy through past radiations and extinctions. However, in many neocheilostome species, ovicells are visibly clustered or otherwise nonrandomly arranged, and most neocheilostome fossil specimens are only partially preserved. Furthermore, new AI methods exist for the auto-recognition of ovicells from SEM images (Di Martino et al. 2023), but older SEM images achieved high magnification by capturing only a subset of the colony. To expand the utility of the neocheilostome study system to more taxa and worse fossil preservation, we must understand the spatial distribution of reproductive labor in a colony and its effect on statistical inferences.

We find that ovicells in the neocheilostome bryozoan Parasmittina eccentrica are clustered and most dense at mid-distances from the ancestrula. To our knowledge, ovicell arrangement has not previously been quantified, so we cannot infer how representative these results may be of other neocheilostomes. From our personal observations, however, we suspect P. eccentrica has a middling degree of nonrandom ovicell arrangement; we have observed that, for example, Dendrobeania lichenoides has sharply defined bands of ovicells in the seasonal North Pacific. Several developmental and environmental factors likely determine ovicell arrangement. Ovicell development is induced via hormone release, and hormones diffuse through adjacent zooids (Ostrovsky 2013). Hormones are most concentrated in zooids near the site of release, so ovicell clustering is expected even in the absence of adaptation. However, adaptive responses to environmental cues such as seasonality (Lombardi et al. 2006) or stress (Swezey et al. 2017) would only heighten a pattern of spatial clustering. Similarly, the relationship between ovicell density and distance to the ancestrula is likely due to a combination of developmental limitations and adaptation. The first several generations of zooids around the ancestrula in neocheilostomes are typically smaller and cannot sexually reproduce (Ostrovsky 2013). However, the P. eccentrica colonies used in this study exhibit secondary frontal budding, sometimes around the ancestrula, and still the pattern of low ovicell density within 8 mm of the ancestrula persists. These specimens are also from the Caribbean coast of Panama, which has very little seasonal variation in temperature or nutrient concentrations (Jackson and O’Dea 2023), so seasonality is unlikely to be a strong control on their reproductive timing. Instead, the concentration of ovicells at 10–18 mm from the ancestrula could result from the optimization of resource sharing via diffusion among zooids. Future work could model the optimal ovicell arrangement for energy flow in the absence of external environmental selective pressures. We demonstrate that nonrandom ovicell arrangement is quantifiable and statistically significant, but we do not expect the pattern of ovicell arrangement found in this population of P. eccentrica to be generalizable across taxa or habitats.

Nonrandom ovicell arrangement means that the presence or absence of an ovicell on each zooid is non-independent. Standard methods of estimating sampling error of a proportion (e.g., a binomial 95% confidence interval) assume independence. While statistical methods for dealing with non-independence do exist (Zimmerman et al. 1993), ovicell arrangement and therefore the nature of the non-independence vary too much among species and environments to be predictable. Further experimental, observational, and simulation work on the reproductive physiology of extant neocheilostomes under diverse biotic and abiotic conditions would help constrain the amount of ovicell non-independence we should expect in fossils.

For this population of P. eccentrica, we find estimates of colony fecundity from a single fragmented specimen have higher variance (and therefore error) than estimated by methods that assume independence. As a result, method 3, which takes the weighted average ovicell density among fragments and uses the observed variance among these fragments to estimate sample variance, is the best estimator of population average fecundity among the three methods we evaluated (Table 1; Fig. 6). These results can provide guidance on whether robust inference is possible for a given research question. For example, if one wanted to compare the fecundity of two or more fragmented (or only partially imaged) colonies engaged in competition for space, the estimate of the fecundity of each individual colony will have erroneously high error. The degree of error will depend on the strength of ovicell clustering and the fragment size. However, if one wants to compare the average fecundity of two or more groups (the groups could be defined by taxon, functional morphology, paleoenvironment, overgrowth status, time bin, etc.), then robust estimates are attainable using method 3. Most research questions can be framed as a comparison of trait means between populations.

Limitations and Future Directions

We have dealt with what we see as the most pressing obstacle to utilizing this model system, but there are further methodological puzzles. In this study, we always sampled an equal number of fragments from each colony. When sampling the fossil record, we rarely know whether fragments come from the same or different colonies. Fossil preservation is patchy both spatially and temporally, so we can assume that even adjacent colonies are not equally represented in the fossil record. The effect of disproportionately sampling certain colonies is another interesting methodological challenge worth further investigation. Furthermore, in this study, all colonies had reached sexual maturity (Jackson and Cheetham 1990). In real fossil assemblages, colony age at death would vary according to an unknown survivorship curve. Ovicell arrangement varies with colony age and size, but we do not understand colony development well enough to model it. If we could model the colony survivorship curve and model how ovicell arrangement varies with colony age, we could account for how these factors affect our estimate of fecundity. Although explaining ovicell density is not the purpose of this study, we do find a positive correlation between colony size (number of zooids) and ovicell density (proportion of zooids with ovicells) (Supplementary Fig. 3; R ² = 0.22, p = 0.009).

Ovicell density is a valuable proxy but not equivalent to fecundity and fitness. Ovicells only represent the sexual fecundity of neocheilostome colonies. Neocheilostomes may also reproduce by asexual fragmentation, and the frequency of asexual propagation is in some cases stable for millions of years (Thomsen and Håkansson 1995), while in other cases it is rapidly evolving (O’Dea and Jackson 2009). Therefore, when using sexual fecundity as a proxy for fitness, one must also control for asexual progeny, especially in arborescent growth forms. Finally, fecundity is only one component of fitness. How to estimate offspring survivability in fossil assemblages has also long puzzled paleobiologists (Van Valen 1963), although Liow et al. (2024) have used success in inter-colony competition for space as a proxy for survivability in cheilostome bryozoans. Despite these caveats, cheilostome fossils have the benefit of large sample sizes and modern analogues on which we can experiment. We are gaining the methodological tools to make robust use of this exciting and underutilized model system, opening doors to test a number of hypotheses about microevolutionary drivers of macroevolutionary patterns.

Conclusions

Larval brood chambers in neocheilostome bryozoans, known as ovicells, are both fascinating examples of colony integration in animals and the rare system with which we can estimate fecundity in the fossil record. In Parasmittina eccentrica, we provide a proof of concept that ovicell arrangement can be quantified and propose ovicell arrangement should be explored across further taxa and environments. When ovicell arrangement is unknown, as will often be the case in fossil taxa, we demonstrate that the observed variance in ovicell density within a population is the best estimate of error. Despite nonrandom ovicell arrangement, paleobiologists can still robustly estimate fecundity with even highly fragmented specimens.