Introduction
The study of the development of fossilized organisms can elucidate both the state of developmental characters at ancient phylogenetic nodes and their subsequent evolution, and also permits exploration of how mechanisms controlling developmental processes have themselves evolved (Hughes et al. Reference Hughes, Hong, Hou and Fusco2017: p. 11). Extracting reliable developmental information from the fossil record requires that the data meet stringent standards (Hughes et al. Reference Hughes, Adrain, Holmes, Hong, Hopkins, Hou, Minelli, Peng, Park, Paterson, Webster, Zhang, Zhang and Fusco2020). Rarely, specimens of articulated trilobites that collectively span a wide size range are well preserved in large numbers within a short stratigraphic interval. Such cases almost invariably include specimens preserved in claystone that are flattened to various degrees, but careful selection of sites and specimens permits informative quantitative analysis of the manner of the growth within individual taxa and glimpses into the ways in which their development was regulated. Initial work in this area focused on a ~429-Myr-old Ma aulacopleurid trilobite, Aulacopleura koninckii, from one locality in the Czech Republic (Fusco et al. Reference Fusco, Hughes, Webster and Minelli2004, Reference Fusco, Hong and Hughes2014, Reference Fusco, Hong and Hughes2016). Recently, similar studies have been conducted on other species and at other sites, including some of Cambrian age (e.g., Du et al. Reference Du, Peng, Wang, Wang, Tian, Wang and Zhang2019, Reference Du, Peng, Wang, Wen and Liu2020; Holmes et al. Reference Holmes, Paterson and García-Bellido2020; Hopkins Reference Hopkins2020). Although all such cases have limitations, including the rather specific set of paleoenvironments represented, these studies do offer important insight into how ancient developmental processes were regulated in this early euarthropod clade. In parallel, and building on this descriptive and analytical work, investigations have begun on mechanistic modeling of growth and segmentation, in which ontogeny reconstruction is based on the application of estimated growth and segmentation parameters, rather than on simple data fitting (Hughes et al. Reference Hughes, Hong, Hou and Fusco2017; Hopkins Reference Hopkins2020). These works suggest the promise of trilobite studies for the discipline of “paleo-evo-devo” (Minelli and Fusco Reference Minelli and Fusco2008: p. 216).
This paper is one of a series on a particularly abundant, late early Cambrian oryctocarine trilobite, Oryctocarella duyunensis (~514 Myr old), from slope facies of the South China block at Bulin in western Hunan. This species is also similarly preserved in at least two other localities but with slightly different phenotypic features, offering the future possibility of exploring ancient microevolutionary changes in developmental patterns. Herein we document several aspects of juvenile (meraspid) growth and segmentation and their regulation in this animal, including axial growth gradients, that is, the distribution of growth rates along the main body axis. Fusco et al. (Reference Fusco, Hong and Hughes2014) identified a growth gradient in the trunk of A. koninckii and later also in the meraspid cephalon, and in the cephalon and trunk during the adult (holaspid) period (Fusco et al. Reference Fusco, Hong and Hughes2016). In the present paper, by using an analysis modified from Fusco et al. (Reference Fusco, Hong and Hughes2014), we explore the existence and the nature of growth gradients along the main body axis of O. duyunensis during the entirety of post-protaspid development in a case in which axial growth rates of tagmata (cephalon and trunk) were not constant across the ontogeny.
From a combination of these analyses, a developmental mode previously undescribed for trilobites has been revealed. Like many extant arthropods (Minelli and Fusco Reference Minelli, Fusco, Minelli, Boxshall and Fusco2013), trilobites displayed hemianamorphic development, that is, they added new trunk segments during postembryonic life until a fixed adult number was reached at a given stage (anamorphic phase), after which growth and molting continued but without any further increase in segment number (epimorphic phase) (Hughes et al. Reference Hughes, Minelli and Fusco2006; for an outline of trilobite development, see Hughes et al. Reference Hughes, Hong, Hou and Fusco2017). Trilobites are also considered to have had indeterminate growth, that is, a form of postembryonic development that in arthropods corresponds to the absence of a predetermined terminal molt (Minelli and Fusco Reference Minelli, Fusco, Minelli, Boxshall and Fusco2013), with the holaspid period characterized by multiple adult-to-adult molts (Hughes et al. Reference Hughes, Minelli and Fusco2006). Oryctocarella duyunensis ontogeny apparently included a fixed terminal molt, likely evolved through progenetic paedomorphosis, and thus qualifies as the first documented case of a trilobite with determinate growth.
Materials and Methods
Materials
The Balang Formation, which hosts the great majority of Oryctocarella duyunensis, outcrops sporadically within an area of some 15,000 km2 in eastern Guizhou and western Hunan and was part of the slope facies of the South China block that was located in peripheral equatorial Gondwana during early Cambrian time (Domeier Reference Domeier2018). Near Balang, it is some 300 m thick and constitutes a series of mudstones and siltstones differentiated most evidently by color and amounts of carbonaceous and carbonate material (Lei Reference Lei2016; Zhao et al. Reference Zhao, Peng, Yuan, Esteve, Yang, Wu and Chen2019; Du et al. Reference Du, Peng, Wang, Wen and Liu2020). The total range of O. duyunensis within the Balang Formation is as much as 290 m (Zhao et al. Reference Zhao, Peng, Yuan, Esteve, Yang, Wu and Chen2019: fig. 2; Du et al. Reference Du, Peng, Wang, Wen and Liu2020: fig. 1), and it also extends into the overlying Qingxudong Formation, but it is most common within an interval of about 150 m in the middle to lower part of the Balang Formation. At Bulin, the Balang Formation is about 200 m thick, and the great majority of our collections were recovered from interbedded argillaceous, arenaceous, and calcareous mudstones in an interval only 4 m thick in the lower part of Balang Formation at the Bulin section, 6.3 km northwest of Jiwei village, Huayuan County, Hunan Province, South China, GPS coordinates 28.355°N, 109.384°E. Further locality details are provided in Dai et al. (Reference Dai, Hughes, Zhang and Pengin press). Biostratigraphically, these fossils occur in Cambrian Series 2, Stage 4, and possibly also in Stage 3, depending on how the boundary between those stages is ultimately defined (Zhao et al. Reference Zhao, Peng, Yuan, Esteve, Yang, Wu and Chen2019: fig. 2).
A total of more than 1700 specimens of O. duyunensis, including 1276 complete specimens, were recovered during fieldwork in the interval 2012–2014. Of these complete specimens, three subsets in this dataset were identified for analysis. Dai et al.'s (Reference Dai, Hughes, Zhang and Pengin press) sketch of the ontogeny of O. duyunensis was based on a dataset of 968 specimens that permitted measurement of their sagittal dorsal length and the number of thoracic segments (dataset 1) and on a subset of this dataset composed of 643 specimens, for which segment number in the pygidia could also be counted with confidence (dataset 2).
Here, based on 337 specimens drawn from dataset 2, for which a complete set of landmark-based measurements were taken (dataset 3), we quantitatively investigate the absolute (stage-based) postembryonic development of this trilobite. The same dataset 3 will be used in a forthcoming study on the relative growth (allometry) based on geometric morphometrics.
Data Acquisition
Body length (BOL, dataset 1), was taken with a caliper (precision 0.05 mm) directly on the specimens. The number of thoracic segments (NTH, dataset 1) and the number of pygidial segments (NPY, dataset 2) were counted directly on the specimens as well.
The fossils of dataset 3 were photographed with a Canon 5Ds Digital SLR camera equipped with a Canon EF-S 60mm 1:2.8 macro lens, in lower-angle lighting from the northwest direction and higher-angle lighting from the northeast direction, or, for the specimens smaller than 5 mm in length, with a Leica M205C stereomicroscope with a Plan-Apo 1.0× lens, and the associated Leica Application Suite v. 4.10 software. The resulting images were digitized using the ImageJ software package (Abràmoff et al. Reference Abràmoff, Magalhaes and Ram2004) with the x- and y-coordinates recorded for each of a series of landmark points on each specimen. A scale in 0.5 mm divisions was included in each image. Because the sagittal position of the posterior of the axial rings is evident in O. duyunensis, axial landmarks were positioned with the aid of a line drawn along the sagittal axis on the image of each specimen.
Morphometric Characters of Dataset 3
Morphometric characters of dataset 3 were extracted from landmark data (Fig. 1). One, the cephalic centroid size (CCS, based on a suite of 30 unreflected cephalic landmarks), is a standard geometric morphometric measure used as a main reference in most size analyses (e.g., Fusco et al. Reference Fusco, Hughes, Webster and Minelli2004). The others are linear distance measures along the main body axis. The length of the cephalon along the main axis (CEL) was partitioned into a frontal piece (FRP) and six cephalic segments (CS1–CS6) from the anterior backward. The whole trunk axial length (TRL) was partitioned into the lengths of individual thoracic segments (TS1–TSx, the actual number varying among specimens), pygidial segments (PS1-PSx, the actual number varying among specimens), and a terminal piece (TEP). Thoracic length (TXL) is the sum of all thoracic segment lengths, and pygidial length (PYL) is the sum of all pygidial segment lengths plus the terminal piece, thus TXL + PYL = TRL. Finally, CEL + TRL = ATL, that is, axial total length. Because of a notch in the terminal margin of the pygidium (Fig. 1), ATL (measured on images, dataset 3) does not correspond exactly to body length (measured on the material, dataset 1), the former being slightly shorter. Relative position of the posterior boundary of a given segment (either cephalic or trunk) with respect to the body region to which it belongs (cephalon or trunk, respectively) or to the whole body was computed as ratios of the relevant distance measurements. Details on computations are given in the Supplementary Material.
Figure 1. Landmark configuration used in this study of Oryctocarella duyunensis. This specimen (NWU-DYXJT 1282) has 11 thoracic segments and 3 pygidial segments. Diamonds and the two squares represent axial landmarks; circles are cephalic exsagittal landmarks. The most anterior square marks the axial junction of the cephalon and the trunk. It and all other cephalic landmarks were used in the calculation of cephalic centroid size. The most posterior square marks the axial junction of the thorax and the pygidium. Arrowheads mark the pleural thoracic–pygidial boundary. Colors highlight the three main axial regions: cephalon (green), thorax (pink), and pygidium (orange).
Morphometric Analyses
Adhering to Hughes et al. (Reference Hughes, Adrain, Holmes, Hong, Hopkins, Hou, Minelli, Peng, Park, Paterson, Webster, Zhang, Zhang and Fusco2020), we refer to different segmental conditions of the specimens in the sample by the term morph. Morphs are labeled as mX, where X is the number of thoracic segments (for meraspids corresponding to the traditional concept of degree, mX ≡ dX), or as mXY, where Y is the number of pygidial segments. Morphs (as well as degrees), established on a strictly descriptive basis, should not be confused with stages, which are developmental steps in the ontogenetic progression separated by molting events (Hughes et al. Reference Hughes, Adrain, Holmes, Hong, Hopkins, Hou, Minelli, Peng, Park, Paterson, Webster, Zhang, Zhang and Fusco2020). Assignment of a specimen to a stage requires an interpretative undertaking based on explicit criteria (see next section). Stages are labeled as sX, where X is a progressive ordinal number from zero (considering s0 to be the first post-protaspid stage).
Morphometric analyses were carried out in R (R Core Team 2020) and STATGRAPHICS Centurion (v. 18.1.11). Cephalic centroid size computation and accessory calculations were performed in Microsoft Excel 2010.
Preliminary Assessment for Staging
In arthropods, absolute growth dynamics and segmentation are developmental descriptions based on the ontogenetic succession of stages (or instars). Specimen stage assignment is therefore a critical step in any ontogenetic analysis of absolute growth, and the choice of a given staging criterion needs to be supported by evidence (Hughes et al. Reference Hughes, Adrain, Holmes, Hong, Hopkins, Hou, Minelli, Peng, Park, Paterson, Webster, Zhang, Zhang and Fusco2020).
We took into consideration alternative staging criteria, based on distinct segmentation dynamics, and evaluated them on the basis of observed data. The criterion finally adopted, which assumes the steady release into the thorax of one pygidial segment per stage for a significant section of the meraspid ontogeny, was used as the basis for data ordination in all subsequent analyses.
Absolute Growth Analyses
The average per-stage growth rate (AGR) for a given size variable was calculated as the antilogarithm of the arithmetic mean of the size increments from one stage to the next, after natural logarithmic transformation of the variable (lnGR) (Fusco et al. Reference Fusco, Garland, Hunt and Hughes2012).
Deviation of the ontogeny from the growth progression at a constant rate (Dyar's rule; Dyar Reference Dyar1890) was quantified with the index of conformity to Dyar's rule (IDC; Fusco et al. Reference Fusco, Garland, Hunt and Hughes2012). This metric takes values between 0 (maximal divergence from Dyar's rule, the whole growth realized in a single stage) and 1 (perfectly constant growth rate).
Growth Gradient Detection
Growth gradients are revealed as a differential in the allometric coefficients of serially arranged body parts with respect to a more inclusive body region. Allometric coefficients were calculated as the ordinary least-squares (OLS) regression slope of mean values of the relevant variables for each stage, that is, cephalic segment lengths with respect to CEL and thoracic segment lengths with respect to TRL. The use of means, rather than specimen data directly, is justified by the need to control for the bias of static (within-stage) allometry, which can be quite substantial for more posterior thoracic segments, which were present only during a few stages, because they were released into the thorax relatively late during ontogeny.
Hypothesis Testing for the Implementation of Growth Gradients
Upon the detection of a clear gradient in the trunk, we considered two primary hypotheses for the developmental implementation of the growth gradient in this body region throughout the interval of ontogeny studied that encompass different mechanisms of growth regulation and control. These are the same hypotheses discussed in Fusco et al. (Reference Fusco, Hong and Hughes2014), but in this case allow the trunk growth rate to vary during ontogeny, either as gradual change in growth dynamics or as the result of biased sampling at different stages.
Under the segmental gradient (SG) hypothesis, the thoracic segments represent individual growth fields specified early in ontogeny and each segment, once released into the thorax, grew at a rate with a fixed relationship to the overall trunk growth rate, depending on its position in the sequence. The SG hypothesis is further articulated in (1) SG-R, in which segments grew at a rate that was a constant ratio of the trunk growth rate, and (2) SG-A, in which segments grew maintaining a constant allometric relationship with the trunk. The different values of growth rates among the segments at each stage were determined by the segmental gradient g(i), a discrete function in the domain of natural numbers from 1 to 11 inclusive (the maximum number of thoracic segments observed in the sample), where i is the ordinal position of a segment counting from the anterior. At some point in ontogeny, either at the time of the segment release into the thorax or before, the subsequent growth relationship of each segment was fixed, and trunk growth and segmental size composition of the thorax were derived from the autonomous growth rates of the various segments. The SG hypothesis reflects the standard assumption that ontogenetic allometry results from differential growth rates of distinct body parts (Nijhout Reference Nijhout2011).
Under the alternative, trunk gradient (TG) hypothesis, the whole trunk was a growth field that exhibited a continuous growth gradient. Growth patterns of segments thus derived from the global growth pattern of the trunk. The different growth rates at each relative position along the trunk determined the scaling trunk gradient g(x), a continuous function in the domain of real numbers from 0 to 1 inclusive, where x is the relative position of a point along the trunk, from the anterior to posterior. During ontogeny, individual segments changed their relative position (an interval of x-values) along the trunk as a result of the differential growth of the various sections of the trunk. Accordingly, as they were subjected to different values of the gradient, their growth rates changed as well. Segment growth under the TG hypothesis involves a form of positional specification along the main body axis, that is, the regulation of tissues’ activity (in this case, axial growth) according to their positional values within a developing field (Wolpert Reference Wolpert2011).
The competing growth implementation hypotheses have different structures, and none are nested within any other; thus, they have been compared with each other using the corrected Akaike information criterion (AICc), based on the predictions of each hypothesis with respect to observed relative size of trunk segments across ontogeny. An initial testing was based on thoracic segment data only, while a second testing was extended to all trunk segments, irrespective of their thoracic or pygidial identity. Further investigations of the gradients, including the possibility of a single global body gradient for the whole body axis, were carried out on other relative-size data derived from the same dataset. A nonlinear least-squares regression procedure was performed using Marquardt's algorithm as an estimation method (implemented in STATGRAPHICS). Details on fitting function derivations and fitting procedure implementations are given in the Supplementary Material.
Investigation of the Holaspid Period
Alternative hypotheses on the number of holaspid stages were evaluated on the basis of cephalic size distributions and morph frequency of occurrence.
Results
Preliminary Assessment for Staging
We used lnCCS as the main logsize variable, because it is considered to be relatively unaffected by anamorphosis (the postembryonic addition of body segments; Fusco and Minelli Reference Fusco and Minelli2021), as cephalic segment number remained constant throughout postembryonic ontogeny.
Preliminary data inspection showed that lnCCS does not exhibit the kind of multimodal distribution that would allow reliable stage assignments based on size clustering, because, especially in later ontogeny, within-stage size distributions are apparently largely overlapping (Hunt and Chapman Reference Hunt and Chapman2001; Webster Reference Webster2015).
The criterion that equates stage with the segmental condition defined by the number of thoracic segments and in suborder by the number of pygidial segments (i.e., with the morph; e.g., Dai et al. Reference Dai, Zhang, Peng and Yao2017) was not taken into consideration, because it would result in a decrease in the number of trunk segments between some consecutive stages (e.g., stage identified by morph m57 would have more trunk segments [12] than the following putative stage m65 [with 11]), something unknown in living arthropods.
Accordingly, we contrastively tested two staging criteria, both based on segmentation dynamics that are regular in some way and on the assumption of a constant per-stage growth rate of a representative size variable or at least a regular change in growth rate across ontogeny that allows a linear approximation in a small subset of stages. These are:
1. Steady anamorphosis (SA). The criterion of stage assignment is the number of trunk segments (sX ≡ NTR). It assumes the steady addition of one trunk segment per stage. The pace of thoracic segment release can vary across ontogeny and among individuals.
2. Steady release (SR). The criterion of stage assignment is the number of thoracic segments (sX ≡ NTH). It assumes the steady release into the thorax of one pygidial segment per stage. Segment addition can vary across ontogeny and among individuals. This corresponds to the conventional trilobite degree-based staging, often used by default.
By using lnCCS as a logsize variable, we applied AIC to compare the two staging criteria on the basis of their performance in producing a regular ontogenetic size progression on the whole ontogenetic series (n = 337; Fig. 2).
Figure 2. Ontogenetic progression of cephalic centroid size (lnCCS) under two alternative staging criteria. A, Steady anamorphosis (SA, stage ≡ number of trunk segments, NTR); and B, steady thoracic segment release (SR, stage ≡ number of thoracic segments, NTH). Boxes represent the interquartile interval, with median (transverse line) and mean (diamond); vertical lines are ranges of variation, to the exclusion of outliers (dots). SR is the better supported criterion (see text).
SR has a normalized probability p > 0.9999 of being the right model in comparison to SA. When restricting the comparison to morphs m4–m9 (n = 277), the normalized probability does not change. Also, the ontogenetic progression of within-stage variation is more regular under the SR criteria.
Based on these results, we adopted the SR criterion for all subsequent elaborations, but with the following caveats: (1) the number of thoracic segments might not correspond to the stage during the holaspid period, when multiple stages maintain the same NTH (although this caveat proved unnecessary in the particular case of Oryctocarella duyunensis; see “Holaspid Period”); (2) the CCS growth progression (see “Per-Stage Growth Rates”) suggests that the level of confidence in degree-based staging varies during ontogeny. In particular, morph m0 has very small sample size associated with large variation in body size (Fig. 3A, dataset 1) and number of pygidial segments (2 to 5, dataset 2), thus hinting that more than one stage may have been spent with 0 thoracic segments. Accordingly, the more sensitive analyses were restricted to that section of ontogeny in which confidence in the adopted staging was highest.
Figure 3. Ontogenetic progression of body length (A, lnBOL), cephalic centroid size (B, lnCCS), cephalic length (C, lnCEL), and trunk length (D, lnTRL). lnBOL is based on dataset 1, the others on dataset 3. Dots are means and bars are standard errors (often not visible because they are smaller than the dot symbol).
Per-Stage Growth Rates
BOL (dataset 1), CCS, CEL, and TRL (dataset 3) grew at average rates of 1.23, 1.18, 1.17, and 1.39, respectively. However, these four size variables did not grow at a constant rate during meraspid ontogeny. They exhibit a first phase of higher growth increment between s0 and s3, little growth between s3 and s4, fairly constant growth but lower than the first phase from s4 to s9, and finally a smaller increment still from s9 to s11 (Fig. 3, Table 1). This apparently unparsimonious description of growth data, is justified by (1) small standard errors in mean size of all stages, such that the minimal growth between s3 and s4 cannot be treated as a sampling artifact; (2) negligible improvement in fitting using more complex polynomials able to capture a progressive change in growth rates with ontogeny; and (3) connections to the segmentation process (see “Segmentation”).
Table 1. Average growth rates for body length (BOL), cephalon centroid size (CCS), cephalic length (CEL), trunk length (TRL), and pygidial length (PYL). BOL based on dataset 1, the others on dataset 3. Estimates are followed by standard errors (not calculated for s0–s3 for dataset 3 as there is only one s0 specimen).
For the cephalon, CCS growth between s3 and s4 was small but significant (p = 0.0144, two-tailed Student's t-test), and growth from s9 to s11 was significantly smaller than previous stages (p = 0.0050, one-tailed Student's t-test). Overall (12 stages) index of Dyar's conformity was quite low (IDC = 0.86) in comparison to other trilobites (Fusco et al. Reference Fusco, Garland, Hunt and Hughes2012), although higher in the stage interval s4–s9 (IDC = 0.92).
For the trunk, TRL, the growth rate between s0 and s3 was high and declined slightly across these stages. Overall IDC was very low (0.70), but higher in the stage interval s4–s9 (IDC = 0.95).
The meraspid pygidium, which exhibited within-stage variation in the number of segments (see “Segmentation”) grew in length for a significant part of the ontogeny, somehow reproducing the complex pattern seen the other size variables despite its changing complement of segments (Fig. 4A, Table 1). However, its growth was not significant in the stage intervals s3–s4 and s9–s11 (two-tailed Student's t-tests). Pygidial relative length with respect to trunk length declined steadily across ontogeny, approximately exponentially (OLS regression, r 2 = 0.997) with an allometric coefficient with respect to TRL of 0.524 (OLS regression, r 2 = 0.835) (Fig. 4B).
Figure 4. Ontogenetic progression of pygidium length (A, lnPYL) and relative pygidial length with respect to trunk length (B, PYL/TRL). Dots are means and bars are standard errors (often not visible because they are smaller than the dot symbol).
Segmentation
Based on dataset 2, within-stage variation in NPY is observable in all stages with a significant sample size (Fig. 5).
Figure 5. Ontogenetic variation in the number of pygidial segments (NPY; based on dataset 2, n = 643). Point size is proportional to the frequency of specific values.
The process of segment addition seems to have been triphasic. A first phase, until s2, consisted of the addition of more than one segment per stage, with individual variation in the rate of addition. This produced the increase in modal NPY and set the within-stage variation in NPY observed from stage s2 onward. Because of the small sample size of stages s0 and s1, it is not possible to assess whether NPY variation first originated during the protaspid period. A second phase, from s2 to s7, saw the segment addition of one segment per stage (according to this interpretation expected morph m77 is absent from the observed dataset). This, in association with segment release, which was occurring at the same pace, maintained both the modal number of NPY and kept variation in NPY constant. A third phase, from stage s7 onward, saw the end of segment addition and a decline in NPY at a pace of one segment per stage. The anamorphic phase apparently ended at stage s6, and thus the developmental mode of O. duyunensis qualifies as protomeric (more precisely, hypoprotomeric; Hughes et al. Reference Hughes, Minelli and Fusco2006).
Holaspid Period
The holaspid period encompasses development after completion of the thorax (Raw Reference Raw1925). In trilobites, the holaspid period is generally characterized by one morph with the highest number of thoracic segments observed in the species that shows a disproportionately wide size range compared with those of preceding (meraspid) stages. This is because multiple stages were spent in the holaspid period, in which, by definition, the number of thoracic segments is fixed. This growth pattern, with no signs of a predetermined terminal molt is called indeterminate growth (Minelli and Fusco Reference Minelli, Fusco, Minelli, Boxshall and Fusco2013), and in extant arthropods is common among myriapods and several crustacean taxa.
This pattern is not observed in O. duyunensis from Bulin as described herein (and not from other localities either; see “Discussion”). The hypothesis that the sample of morph m9 specimens (n = 95) corresponds to a single developmental stage (s9) can be supported by different lines of evidence based on dataset 1. First, morph m9 does not exhibit a significant increase in lnCCS size variance with respect to m8 (Fisher's test, F = 1.51, p = 0.08; Fig. 2A). Second, the frequency of occurrence of m9 (23% of the total sample) is similar that of morph m8 (18%). Third, assuming (1) a constant logsize increment from stage s8 onward (average lnGR = 0.185, estimated from s6 to s8), (2) a constant within-stage lnCCS variance equal to that observed in stage s8 (var(lnCCS) = 0.016), (3) an equal frequency of occurrence of specimens from stage s9 on, it is possible to apply a likelihood-ratio test to see whether the observed lnCCS distribution of m9 is more probable under the hypothesis of a single stage (s9) or the hypotheses of two or more stages (s9 + s10, s9 + s10 + s11, …). There is a modest evidential ratio in favor of the hypothesis of one stage versus the hypothesis of two stages (Λ = 1.62) and very strong evidence for the hypothesis of one stage versus any hypothesis of more than two stages (Λ > 1011). The likelihood ratio of the hypothesis of one stage versus the hypothesis of two stages is very sensitive to the parameters used for estimating the size distributions under the two hypotheses, therefore this test, while allowing us to reject with confidence the hypotheses with more than two m9 stages, does not decisively exclude the possibility that m9 is a mix of two stages.
Taking these lines of evidence together, the most plausible scenario is that O. duyunensis exhibited a single holaspid stage (instar) with a polymorphism in the number of mature thoracic segments (9–11), produced by one or two extra molts associated with growth and segment release (but not to segment addition) after stage s9 (Fig. 6). Considering that in dataset 1 (n = 968), 222 specimens have NTH = 9, 31 have NTH = 10 and 3 have NTH = 11, and that there is no marked variation in the frequency of occurrence between morphs m8 and m9, we can estimate that about 14% (31/222) of individuals underwent one extra molt and achieved a final number of 10 thoracic segments and that about 1% (3/222) did that but also underwent a further molt to achieve a total of 11 thoracic segments. Accordingly, about 15% (34/222) of the individuals with NTH = 9 and 10% (3/31) of the individuals with NTH = 10 were actually meraspids.
Figure 6. Schematics of the ontogenetic progression of the number of trunk segments (NTR, upper broken line) and the number of thoracic segments (NTH, lower straight line) for an individual with a modal number of pygidial segments (NPY = NTR − NTH). Only a minority of individuals entered stages s10 (medium-thickness line) and s11 (thin line).
Irrespective of the possibility of one or two stages with nine thoracic segments, the results presented herein suggest that all O. duyunensis from Bulin had determinate growth, that is, a development mode that comprised a predetermined terminal molt for the individual, associated with some modest variation in the population for the stage (s9–s11) and the number of thoracic segments (9–11) at which molting ceased.
These calculations are supported by taphonomic considerations. Up to a stratigraphic resolution within 2 cm, specimens with 9, 10, and 11 thoracic segments occur together and may have been contemporaries. The different frequencies of occurrence of the different morphs are not easily related either to stratigraphy or to taphonomy.
Growth Regulation
Variance in lnCCS shows a complex ontogenetic pattern (Fig. 7).
Figure 7. Ontogenetic progression in lnCCS variance. Bars are standard errors. Variance not computable for stage s0 with a sample size of n = 1; standard error not computable for stage s1 with a sample size of n = 2.
Restricting the test for the homogeneity of variance to stages s2–s10, which have a sample size ≥12, differences are significant (Levene's test, p < 0.0001). Variance decreases significantly (p < 0.0001) from stage s2 to s3, increases significantly (p < 0.0001) from stage s3 to s9, and then declines, although not significantly (p = 0.0713), from stage s9 to s10.
There is no evidence of steady compensatory growth across the whole ontogeny, although a possible indication of compensation might occur within the short stage interval s5–s8 (Levene's test, p = 0.69).
Growth Gradients
Segmental Growth Gradient Detection
Cephalic segmental allometric coefficients with respect to CEL show oscillating values (Fig. 8), resulting in a nonsignificant gradient (Spearman's rank correlation test, r = −0.14, p = 0.75, n = 6). The thorax exhibits a significant gradient (Fig. 8), with declining segmental allometric coefficients with respect to TRL from posterior to anterior (Spearman's rank correlation test, r = 0.9879, p = 0.0030, n = 10). The gradient is significant even excluding thoracic segment 10, where standard error cannot be calculated (r = 0.9833, p = 0.0054, n = 9).
Figure 8. Growth gradient in the cephalon and trunk. Cephalic segment ontogenetic allometric coefficients with respect to cephalic length (left) and thoracic segment ontogenetic allometric coefficients with respect to trunk length (right), calculated as ordinary least-squares (OLS) regressions (with n = 12 for the cephalon, and with n decreasing from n = 11 for TS1 to n = 2 for TS10). Bars are standard errors (not calculable for TS10). Declining values from posterior to anterior are significant only in the trunk (Spearman's rank correlation test).
Trunk Growth Gradient Hypothesis Testing
Stage s0 has no thoracic segments and stage s1 is an outlier under any hypothesis; thus gradient testing was restricted to the developmental interval s2–s11. In this interval, the TG hypothesis has a normalized probability p > 0.9999 of being the correct model when compared with SG-R, and p = 0.8063 when compared with SG-A, with an evidential ratio (ER) on the order of 1035 in the first case and of ER = 4.16 in the second case (Table 2). From here on, the comparison was carried out between the TG hypothesis and the best-performing of the two SG hypotheses, SG-A.
Table 2. Results of Akaike information criterion (AIC) comparison onto relative length of thoracic segments (RLS) dataset for the three competing hypotheses for the developmental interval s2–s11. No. par., number of estimated parameters; n, number of observations (estimated RLS values); ΔAICc, difference in corrected Akaike score with respect to the model with the minimum score; wAICc, Akaike weight, the probability of being the correct model among the set of candidate models.
Nonlinear least-squares regression based on thoracic segment data returned multiple best-fit curves for each of the two hypotheses, depending on the trunk growth rate (Fig. 9).
Figure 9. Estimated growth gradients under the two hypotheses at three different values of observed trunk per-stage growth rates (TRG) in s2–s11: minimum (TRG = 1.10, light color), mean (TRG = 1.27, intermediate color), maximum (TRG = 1.56, dark color). A, Segmental gradient (SG-A) of ontogenetic allometric coefficients of thoracic segment with respect to trunk length. B, Trunk gradient (TG) of the growth rates at each relative position along the trunk as a fraction of the trunk growth rate.
The TG model explains 98.59% of the observed variance in the relative length of thoracic segments (RLS) in the stage interval s2–s11 (n = 54), while SG-R explains 98.45%. The modest difference in the evidential support of the two competing hypotheses depends on the subdued curvature of the growth gradient, either seen as a segmental (barely distinguishable from a straight line; Fig. 9A) or a trunk gradient (Fig. 9B), which results in similar segment growth dynamics in the thorax under the two growth mechanisms (Supplementary Fig. S1). Because of this modestly curved gradient, the profile of the relative lengths of all the segments at any given stage did not change much during ontogeny, with the segmental position of the longest segment in the series moving from the second segment at stage s2 to the third segment at stage s11.
The comparison of the two hypotheses based on trunk segment data (in the same stage interval s2–s11, n = 108) sees a significant increase of the support of TG with respect to SG-A, with a normalized probability p > 0.9999 and an ER = 89697. However, this requires that in the SG-A model, the allometric coefficient of a given segment is established at the moment of its first expression in the pygidium, rather than at the moment of its release into the thorax.
Further Investigations of the Growth Gradient
Assuming TG to be the correct hypothesis, we can fit it to a dataset that is more suited to the dynamics it entails, namely, to the set of relative positions of the boundaries between segments (simply considered as trunk landmarks) irrespective of their relation to the corresponding segment. This produces more accurate parameter estimates for the trunk gradient. The TG model explains 99.97% of the observed variance in the relative position of the posterior boundary of thoracic segments (s2–s11, n = 54) (Fig. 10) and 99.94% of the variance in the relative position of the posterior boundary of trunk segments (s2–s11, n = 108).
Figure 10. Growth gradients under TG hypothesis estimated by fitting of the relative position of the posterior boundary of thoracic segment dataset at three different values of observed trunk per-stage growth rates (TRG) in s2–s11: minimum (TRG = 1.10, light color), mean (TRG = 1.27, intermediate color), maximum (TRG = 1.56, dark color). A, Growth rates at each relative position along the trunk as a fraction of the trunk growth rate. B, Absolute growth rates at each relative position along the trunk.
Again, under the assumption that TG is the correct model for the trunk and that the cephalon might have a very flat (and thus undetected) gradient, we tested whether the growth dynamics within the cephalon can be combined into a single whole-axis gradient, together with the trunk. The dataset consists of the relative positions of the posterior boundary of all segments (cephalic, thoracic, pygidial) rescaled with respect to the length of the whole axial length (AXL). The fit is very good (s2–s11, r 2 = 0.9996, n = 171); however, when using AICc to compare the predictions for the thorax data (s2–s11, n = 54) in the context of the trunk with the thorax data in the context of the whole body, the fit within the trunk results significantly better (ER = 6.5 × 106, p > 0.9999).
Discussion
Average per-stage growth rates for the main cephalic size variables (CCS and CEL) were very close to trilobite median values (1.18 and 1.17 vs. 1.18 and 1.16, respectively; Fusco et al. Reference Fusco, Garland, Hunt and Hughes2012), while the index of conformity Dyar's rule for CCS was relatively low with respect to the median of the group (0.86 vs. 0.94; Fusco et al. Reference Fusco, Garland, Hunt and Hughes2012). Departure from growth at a constant rate, little of which can be attributed to sampling error, included a distinctive “pause” in growth between s3 and s4 that does not correspond to any other developmental event detected. A phase of minimal growth, possibly from s2 to s4, is apparent also in the collection of the same species (under the name Arthricocephalus chauveaui) from the Lazizhai section in Guizhou Province (Du et al. Reference Du, Peng, Wang, Wen and Liu2020; see also Dai et al. Reference Dai, Hughes, Zhang and Pengin press), although sample size is small in that case. Owing to the strong correlation between the main size variables, trunk length and pygidial ontogenetic progressions approximately reproduced the pattern of cephalic size variables, despite the trunk's different growth regime. Pygidial length ontogenetic progression depended on the combination of the incremental effect of axial growth in the posterior of the trunk and the decremental effect of segment release into the thorax. These two processes had virtually independent dynamics across ontogeny, but their combined effect produced a steady decline across ontogeny in pygidial relative length with respect to trunk length that contrasts with the variable pace of the two processes.
An axial growth gradient could be detected with confidence only in the trunk. Comparing the gradient in Oryctocarella duyunensis under the TG hypothesis with the meraspid trunk gradient in the well-studied Aulacopleura koninckii (Fusco et al. Reference Fusco, Hong and Hughes2014, Reference Fusco, Hong and Hughes2016), several differences emerge. First, the O. duyunensis axial trunk gradient is characterized by a marked dependence on trunk growth rate. This was not observed in A. koninckii, for which the trunk growth rate was fairly constant in the investigated ontogenetic interval. Second, when the shapes of the gradients of the two species are compared, normalizing for the difference in average trunk growth rate (Supplementary Fig. S2), that of O. duyunensis is less markedly curved. The disparity between local growth rates at the two ends of the trunk is similar in the two species at their respective average trunk growth rate (TRG = 1.27 for O. duyunensis; TRG = 1.12 for A. koninckii) but decisively smaller in O. duyunensis when the two species are compared at the same growth rate of TRG = 1.12 (Supplementary Fig. S2). These features, associated with a shorter relative length of the thorax with respect to the trunk in O. duyunensis (i.e., the portion of the trunk where the expectations of the two hypotheses can be compared), make the predictions of the two alternative gradient implementation hypotheses (SG and TG) very similar (Supplementary Fig. S1), which results in an overall modest support for the TG hypothesis in O. duyunensis in contrast to the stronger support in A. koninckii. Third, while A. koninckii had separate meraspid cephalic and trunk gradients of opposite anteroposterior polarity, in O. duyunensis, meraspid cephalic and trunk distribution of growth values can fit a single monotonic gradient from the posterior, with no head/trunk divide. Again, this may be an artifact of the almost flat gradient in the anterior of the trunk combined with the possibly nonexistent gradient in the cephalon. In the only other published work reporting on trilobite continuous growth gradients (Hopkins Reference Hopkins2020), the shape of the meraspid trunk gradient of the Cambrian aulacopleurid Elrathia kingii, normalizing for the average trunk growth rate, is very similar to that in A. koninckii.
The developmental mode of O. duyunensis qualifies as protomeric (Hughes et al. Reference Hughes, Minelli and Fusco2006), the most common developmental mode in trilobites (62%; Fusco et al. Reference Fusco, Garland, Hunt and Hughes2012). However, this segmentation pattern is associated with a peculiar feature of this species. Taken at face value, the results presented herein suggest that O. duyunensis from Bulin had determinate growth, along with variation in the number of thoracic segments at which further molting ceased. If this is the case, it is, to our knowledge, the first time that this growth pattern has been detected in trilobites. Given the potential importance of this result, it is appropriate to consider whether it could be a preservational artifact, associated with a selective bias against preserving larger specimens, or some feature of the life cycle of the species that might hide the characteristic pattern of holaspid growth.
Given the depositional setting and common preservation of specimens in exuvial configuration, it is most unlikely that a process of mechanical sorting could have physically removed larger specimens from the sample. One might conjecture a lifestyle change associated with the onset of maturity that resulted in larger individuals migrating away from the habitat of the younger stages, but this is also unlikely for a number of reasons. Physical sedimentological evidence suggests a relatively deep-water, shelf-edge environment, which is consistent with the broadly “olenomorphic” morphotype of the species that characterizes such settings through much of trilobite evolutionary history (Fortey and Owens Reference Fortey, Owens and McNamara1990). A change from a benthic to a nektonic habit would be unlikely to result in the complete absence of larger individuals from the sample in such a setting, and an out-of-habitat migration is both implausible given the setting and inconsistent with other occurrences of this taxon: we know of no locality bearing O. duyunensis only of large size or from a different paleoenvironment. Rather, where O. duyunensis is known from elsewhere (e.g., Du et al. Reference Du, Peng, Wang, Wen and Liu2020), we see not only the same style of preservation of articulated meraspid exoskeletons in claystones, but also a pattern of growth that is likely consistent with the determinate pattern described herein, albeit with locally differing numbers of maximum observed thoracic segment numbers. For example, in the Lazizhai section near Balang, the maximum number of segments observed in an individual belonging to this species is 9 (Du et al. Reference Du, Peng, Wang, Wen and Liu2020), whereas at the Balang type section, a specimen with 11 thoracic segments is preserved with an additional pygidial segment partially released (Peng et al. Reference Peng, Babcock, Zhu, Lei and Dai2017: figs. 9.14, 9.15). These subtle variations in growth mode attest to the fidelity of the regional stratigraphic record in preserving the microevolutionary fabric of body plan evolution, even in rocks more than 500-Myr-old.
The assemblage studied includes some specimens that are almost certainly exuviae (such as those showing the rotated and inverted free cheeks; see Dai et al. Reference Dai, Hughes, Zhang and Pengin press: figs. 3.9, 4.1) but also others that are likely carcasses, including many of those in which the hypostome is in place. This issue will be discussed in more detail in a forthcoming study, but cases in which specimens can be assigned with confidence as exuviae or carcasses form only a small proportion of the total sample. The merger of carcasses and exuviae in our study limits our ability to access the significance of morph frequency data (see Hartnoll and Bryant Reference Hartnoll and Bryant1990), but the gradual rise in morph frequency to morph m9 and then abrupt drop (Dai et al. Reference Dai, Hughes, Zhang and Pengin press: fig. 9) are unlikely to be an artifact of a mixed sample. This pattern contrasts markedly with the gradual decline in frequency of the largest individual from comparable paleoenvironmental settings (Sheldon Reference Sheldon1988) that is consistent with indeterminate growth. Rather, the pattern seen among Bulin O. duyunensis supports the interpretation of determinate growth.
No detailed phylogeny yet exists to locate O. duyenensis among its closest relatives, but the species belongs to a subfamily of spineless oryctocephalid trilobites, the Oryctocarinae (Dai et al. Reference Dai, Hughes, Zhang and Pengin press). An origin via progenetic paedomorphosis (the heterochronic change in which development stops prematurely with respect to the ancestral condition; McNamara Reference McNamara1986a) has commonly been suggested for this group (Hupé Reference Hupé1953; McNamara Reference McNamara1986b) and receives some support from phylogenetic analyses (Sundberg Reference Sundberg2014). Features consistent with a progenetic origin are the reduced holaspid segment count and small size, both at the onset of trunk maturity and of the largest holaspid specimens (Dai et al. Reference Dai, Hughes, Zhang and Pengin press). Here, we show that growth in O. duyenensis is determinate, such that specimens reach their largest size coincident with the onset of trunk maturity, consistent with the progenetic description.
Determinate growth is considered to be a derived growth mode for arthropods (and possibly for metazoans) that evolved independently in different groups (Hariharan et al. Reference Hariharan, Wake and Wake2016). This growth mode seems to be a relatively labile character, stable in some lineages, but varying among closely related taxa in others (Hartnoll Reference Hartnoll and Bliss1982). Among extant arthropods, determinate growth can be found in many spiders, among millipedes (in the form of teloanamorphosis), in several “crustacean” (non-Hexapoda Pancrustacea) lineages, and in most Pterygota among the insects. Depending on the taxon, sexual maturity can either occur before the final instar (e.g., in Cumacea crustaceans) or with the final instar (e.g., in Polydesmida millipedes) (Minelli and Fusco Reference Minelli, Fusco, Minelli, Boxshall and Fusco2013). Within the Brachyura (crabs), one of the groups in which this aspect of growth has been investigated in some detail (Hartnoll Reference Hartnoll, Castro, Davie, Guinot, Schram and von Vaupel Klein2015), determinate growth is interpreted as a truncated indeterminate growth pattern, with no compensation for the smaller number of mature molts, as it tends to be associated with smaller body size and shorter life span (McLay Reference McLay, Castro, Davie, Guinot, Schram and von Vaupel Klein2015). However, in Brachyura, which are not anamorphic, determinate-growth taxa exhibit the same number of trunk segments typical of the whole clade (20).
We have no way to establish whether in O. duyenensis reproductive maturity was obtained with the last molt or beforehand, or whether intraspecific variation in the total number of stages was associated with parallel variation in the stage of the so-called pubertal molt (the molt to the sexually mature stage), and thus in the number of reproductively mature stages. However, we note that, similar to extant taxa, (1) these features might have evolved in this species and not in other closely related species, (2) determinate growth might have evolved by truncation of the ontogeny (progenetic paedomorphosis), and (3) the stage of the final molt could have been delayed by one or two stages in some individuals, as observed in some determinate-growth crabs (with no variation in the number of trunk segments; McLay Reference McLay, Castro, Davie, Guinot, Schram and von Vaupel Klein2015) and millipedes (with the addition of extra trunk segments with respect to the typical number for the species; Minelli Reference Minelli2020).
In O. duyunensis from Bulin we recognize three varieties of holaspid, with 9, 10, and 11 thoracic segments, respectively. The reduced growth rates between s9 and s11 compared with those of immediately preceding stages may suggest that those individuals molting again after s9 were drawn from the smaller end of the s9 size-range distribution, and similarly for those that ultimately attained 11 segments from the s10 distribution. This speculation is supported by a common form of regulation of the final size in extant arthropods, that is, by altering the number of stages (Minelli and Fusco Reference Minelli, Fusco, Minelli, Boxshall and Fusco2013). At one or more “size checkpoints” along ontogeny, individuals that have not reached a predetermined critical threshold size undergo one or more cycles of growth and molting before entering the next phase of the life cycle (in holometabolous insects, typically the pupa stage; Nijhout and Callier Reference Nijhout and Callier2015). This form of adult size regulation could perhaps compensate for the limited (if any) molt-by-molt compensatory growth observed during the meraspid period. However, why these extra molts should have been associated with an extra thoracic segment release in O. duyunensis remains to be explained.
Individuals in their terminal stage, with 9, 10, or 11 thoracic segments, were holaspids, because the development of their thoraces had been completed (Raw Reference Raw1925). However, because of the pattern of determinate growth, they do not show the notable range of holaspid size that characterizes almost all trilobite species. The ontogeny of O. duyunensis thus provides another demonstration of the variety of developmental modes present among Cambrian arthropods (Hughes Reference Hughes2007), to which determinate growth can now be added.
Data Availability Statement
Data available in the Dryad Digital Data Repository: https://doi.org/10.6086/D1HX1D.
Acknowledgments
We thank J. D. Holmes for fruitful discussion on gradient modeling. M. J. Hopkins, J. D. Holmes, and one anonymous reviewer provided insightful comments on a previous version of the article. T.D. and X.Z.'s contributions were supported by the Natural Science Foundation of China (grant nos. 41621003, 41890840, and 41890844), the 111 Project (grant no. D17013), and the Strategic Priority Research Program of Chinese Academy of Sciences (grant no. XDB26000000). This paper is a contribution toward IGCP 668 project “The Stratigraphic and Magmatic History of Early Paleozoic Equatorial Gondwana and Its Associated Evolutionary Dynamics.” N.C.H.'s contribution was supported by U.S. National Science Foundation grant EAR-1849963 and by the Fulbright Academic and Professional Excellence Award 2019 APE-R/107 held at the Indian Statistical Institute, Kolkata.
Open access
