Ontogeny of highly variable ceratitid ammonoids from the Anisian (Middle Triassic)

Geological setting

The ammonoid material derives from the Fossil Hill Member of Fossil Hill in the Humboldt Range and Muller Canyon in the Augusta Mountains (Pershing County), north-western Nevada, USA (Fig. 1) and is stored in the Geosciences Collection of the University of Bremen (GSUB), Germany. The material from the Wilderness Study Area of the Augusta Mountains, Pershing County was collected with permission of the US Department of the Interior, Bureau of Land Management (BLM, Nevada State office, Winnemucca District). The Fossil Hill Member is a succession of alternating layers of mudstone with lenticular limestone and calcareous siltstone beds (see Fig. 2). The rich and diverse fossil content consists primarily of halobiid bivalves and ammonoids. Detailed geological and stratigraphic descriptions were published by Nichols & Silberling (1977), Silberling & Nichols (1982) and Monnet & Bucher (2005).

Studied specimens

The fossil material comprises 72 ammonoid specimens of the family of Ceratitidae (Mojsisovics, 1879). These represent twelve species in seven genera (Fig. 2, Table 1) that either belong to the subfamily Beyrichitinae (Spath, 1934) or Paraceratitinae (Silberling, 1962). Most of the studied species show high intraspecific variation with overlapping morphologies (see Table 1 and Figs. 3–5). Members of these genera (Gymnotoceras, Frechites and Parafrechites in particular) are sometimes hard to differentiate. They mainly differ in the ventral conch outline, ornamentation, adult ribbing and maximum growth size. The younger the individuals are, the greater the similarities. Despite their complicated taxonomy, all selected species are index fossils of the late Anisian Fossil Hill Member (see Fig. 2). It was assumed that the individual species have similar coiling rates (i.e., the individual species develop the same number of whorls in the course of their life). The total number of volutions developed by the species varies between five and a half and seven (see Table 1—Total number of volutions).

Preparation and data acquisition

We prepared high-precision cross-sections intersecting the protoconch of each specimen, following the methods by Korn (2010), Klug et al. (2015a) and Klug et al. (2015b). Subsequently, we scanned the polished surfaces in high resolution with a flat screen scanner to ensure that all pictures have the same scale. Thereafter, the scan images were digitized. CT scan images of Anisian ammonoids from Nevada do not provide sufficient contrast of the internal structures for a reliable analysis (Bischof & Lehmann, 2020).

Based on the digitized cross-sections, we performed a 2D landmark-based geometric morphometrics analysis. The landmarks were retrieved in tpsDig2 v.2.31 (Rohlf, 2010). Sixteen landmarks were digitized per half whorl (i.e., whorl stage), which resulted in 176 landmarks per specimen (16 landmarks on 11 half whorls; Fig. 6). This set of landmarks consists of two single (1, 2) and 7 pairs of landmarks (3–16), of which eight are sliding semi-landmarks. Whereas landmarks are discrete anatomical loci (i.e., point of highest curvature of venter), sliding semi-landmarks are placed along a curve (or a surface) between two landmarks in a way that best describes the curvatures of the outline. In a second step, an algorithmic approach optimizes the approximation of the outline (Zelditch, Swiderski & Sheets, 2012).

In order to omit missing values in subsequent analyses, the data set was limited to whorl stage number 5.5. From a methodological point of view, it is more practical to rotate the shells by 90° compared to conventional illustrations (cf. Stridsberg, 1990) into a lying position. Since ammonoid conchs are spiral-shaped, each whorl is cut in two parts when preparing the cross sections (see Fig. 6). The half whorls on the left side of the protoconch have odd numbers (“odd whorls”; here 0.5–5.5) and those on the right side have even numbers (“even whorls”; here 1.0–5.0”). Homologous landmarks were set in accordance to the axial plane.

Calculation of procrustes shape

All geometric morphometric analyses were carried out using the R software v 3.6.3. (R Core Team, 2020) packages Morpho v2.8 (Schlager, 2017), geomorph v3.3.1. (Adams et al., 2020) and RRPP v0.6.0 (Collyer & Adams, 2018; Collyer & Adams, 2020). Plots were drawn with the R package ggplot2 (Wickham, 2016). Using the Morpho::procSym function, the 2D landmark coordinates were subjected to a full generalized Procrustes alignment (GPA). The semilandmarks were slid minimizing Procrustes distance. The full Procrustes fit standardizes size, orientation and position, leaving only the Procrustes shape coordinates (Bookstein, 1991, chap. 7.1, p. 258–270; Hammer & Harper, 2005; Zelditch, Swiderski & Sheets, 2012). Since the “odd” and “even” whorls cannot be made congruent by any of these operations (i.e., alignment, translation, rotation), all “even” whorls were manually mirrored before the GPA.

The individual whorls were regarded as different structures of the ammonoid conch. Therefore, the GPA was performed separately for every whorl. The procSym function performs Procrustes superimposition including sliding of semi-landmarks on curves and accounts for the symmetry of the object. Subsequently, the R function geomorph::combine.subsets was used to normalize the configurations of all whorl stages to unit centroid size or with a customized weighting (see “Developmental morphospaces”). The centroid size (CS) is regarded as a proxy for the size of the whorls and equals the square root of the summed squared distances of each landmark from the centroid of the landmark configuration before the GPA (Zelditch, Swiderski & Sheets, 2012). The function geomorph::combine.subsets was originally introduced to combine different parts of a body (e.g., heads and tails; Collyer, Davis & Adams, 2020).

To visualize the multivariate data in two-dimensional morphospaces, we ran a Principal Component Analysis (PCA) on the aligned Procrustes shape coordinates using the R function stats::prcomp. Thereby we used two different types of visualization: Ontogenetic trajectory spaces and developmental morphospaces.

Ontogenetic trajectory spaces

It is well-known that ammonoids have a very characteristic but also complex ontogenetic development (e.g., Klug, 2001). To visualize the ontogenetic development of ammonoids, there are different types of morphospaces. Ontogenetic trajectory spaces (originally called ontogenetic morphospaces), as defined by Bischof & Lehmann (2020), p.2), illustrate the differences in total ontogenetic development of individuals. They show the data in an artificial state of combined morphologies of different ontogenetic stages. To calculate an ontogenetic trajectory space, all Procrustes shapes (i.e., whorls) of an individual are re-assembled before running the Principal Component Analysis. This means that, in an ontogenetic trajectory space, the ontogenetic trajectory of every individual is reduced to a single data point. Ontogenetic trajectory spaces are a tool to examine if the ontogenetic pathways of individuals differ, but they do not show how the trajectories vary. To test whether ontogenetic trajectories statistically differ between species, a multivariate analysis of variance (MANOVA) using the R function stats::manova was applied.

Developmental morphospaces

Developmental morphospaces as defined by Eble (2003, p. 40) are morphospaces that directly contain developmental information. In terms of this study, this means that every individual dot in the morphospace reflects a specific ontogenetic stage (i.e., half whorl) of an individual. By connecting all points of an individual, its ontogenetic trajectory can be obtained. In contrast to ontogenetic trajectory spaces, developmental morphospaces show how individual whorls differ from each other.

General Procrustes Analysis (GPA) removes all information about size from a given set of data leaving only the pure shape coordinates. However, as can be seen in Fig. 6, size differences between different whorl stages are tremendous. If normalized to unit centroid size (i.e., non-weighted morphospace), the earliest whorls of ammonoids therefore get enormously enlarged and the last whorls scaled down. In general, deviations (measurement uncertainties as well as actual morphological variation) are increased for the initial whorls and reduced for older whorl stages. Therefore, a second morphospace with weighted Procrustes shape coordinates was calculated. Thereby, the logarithmic centroid size (log₁₀CS) of all configurations of a whorl stage were normalized to the proportional centroid size of the respective stage to the sum of all whorl stages (log₁₀ CS_{whorl i}/∑ log₁₀ CS_whorls). The principal components of the PCA on the weighted shapes were called wPC (weighted principal components).

If the relative log₁₀CS is used to normalize the centroid size of the configurations, this approach is extremely similar to a Relative Warp Principal Component Analysis (RW-PCA; size-shape space) after Mitteroecker et al. (2004). To calculate a RW-PCA the shape matrix of a configuration is augmented by an additional column containing information about the log₁₀CS of the configurations. Whereas the R function geomorph::combine.subsets scales every configuration accordingly, the size information in the RW-shape matrices are stored in the additional variable. The resulting RW size-shape space can be analyzed with an ordinary PCA. Typically, RW size-shapes spaces are strongly dominated by the log₁₀CS and PC1 therefore often accounts for more than 90% of the variation. If proportional log₁₀CS-values (log₁₀CS_{configuration i}/∑log₁₀CS_{configurations}) are used, the analysis is less dominated by size, but the eigenvalues are very similar to the ones of the weighted PCA (wPCA). For simplicity the R function was used here.

Because weighting does not change the shapes itself, weighted and non-weighted developmental morphospaces look very similar. The main difference is the placement of the individual configurations within the morphospace. Whether weighted or non-weighted shape coordinates should be used, depends on what the analysis is intended to show. A summary of the three different approaches how to combine landmark configurations can be found in Table 2.

To model the shapes at the maximum and minimum PC-values, the R function GeometricMorphometricsMix::reversePCA (Fruciano, 2019) was used. The function is designed to recalculate artificial Procrustes shape variables from the extreme PC-values in a morphospace. The thin-plate spline deformation grids were calculated using the R function geomorph::plotRefToTarget.

Trajectory analysis

In morphometric studies, ontogenetic trajectories represent a series of measurement values of different ontogenetic stages of an individual or a group, called longitudinal data (Klingenberg, 1998). To quantify the differences of the ontogenetic trajectories of the individual species, the R function RRPP::trajectory.analysis with 999 iterations was used. The function calculates a linear model with at least one categorical interaction variable (here: Shape ∼ Species * WhorlStage) and assesses differences in path distance (magnitude differences, length of trajectories), trajectory shape and the angle between the individual trajectories (trajectory correlation) (Collyer & Adams, 2013).

If weighted shape-coordinates were used, the artificial size-shape relationship could overlay true differences between the trajectories. Therefore, only non-weighted shape coordinates were analyzed in the trajectory analysis.

Ontogenetic trajectory spaces

Ontogenetic trajectory spaces are means to visualize whether or not the ontogenetic development of two or more individuals differs. The first three components (PCs) of the Principal Component Analysis (PCA) on the shape coordinates with combined ontogenetic stages of an individual account for 58.3% (PC₁ = 38.3%, PC₂ = 13.4%, PC₃ = 6.6%) of the total variation. Especially considering that there are a total of 352 primary components (x and y coordinates of 176 landmarks), but only 72 specimens, this result can be regarded as satisfactory. The convex hulls of the ontogenetic trajectory space of most species reveal a large overlap (Fig. 7). Nevertheless, ontogenetic trajectory spaces differ significantly between species (MANOVA: Pillai’s trace = 9.8962; ∼ F (1, 11) = 1.6437, p < 0.001).

Since PC1 accounts for 38.3% of the total variation, the most important characteristic is the position of the individuals on the x-axis. In fact, there are certain species that primarily have negative PC1 values (B. vogdesi, F. nevadanus, F. occidentalis, M. spinifer, P. meeki) and some that are more restricted to positive PC1 values (G. blakei, G. mimetus, G. rotelliformis, G. weitschati, P. dunni). B. cordeyi and D. lawsoni, both cover a wider range of different PC1 values, but are generally restricted to negative PC2 and positive PC3 values.

Non-weighted developmental morphospace occupation

The first three components of the PCA on the non-weighted shape space account for 93.8% (PC₁ = 78.5%, PC₂ = 11.5%, PC₃ = 3.1%) of the total variation. The PCA plot (Fig. 8) of PC1 and PC2 shows that the whorls of early ontogenetic stages cover the lower left quadrant of the morphospace (negative PC1 and PC2 values), which characterizes extremely depressed, broad whorls with a flat venter (Fig. 9A). The center of the morphospace (PC1 equals 0 and PC2 is positive) is occupied by intermediate growth stages (juveniles), which have a more quadratic outline with an only slightly triangular venter (Fig. 9B). The lower right quadrant (high PC1 and low PC2 values) is associated with the latest ontogenetic stages (adults). Towards maturity, the whorls increase mainly in height and have a clearly triangular venter and sometimes a keel (Figs. 9C and 10). Overall, there are two extreme adult shapes: Type (A) describes rather depressed, stout conchs with only a slight overlap with the preceding whorl and are associated with much shorter ontogenetic trajectories. Type (B) describes compressed conchs with a clearly triangular venter and a higher degree of overlap and are associated with longer ontogenetic trajectories. For the assignment of the species to the two types see Table 3.

The developmental morphospace of beyrichitine and paraceratitine ammonoids comprises three basic shape stages, which are not separated by sharp borders (Fig. 9): (1) Earliest whorls: broad and very flat; (2) Juveniles: more rounded and depressed; (3) Adults: mostly high and compressed whorls. Since type A species stop their development at more rounded and depressed whorls, their adult whorls resemble the juvenile stages of type B (Fig. 10).

Ontogenetic trajectories in the non-weighted developmental morphospace

The ontogenetic trajectories of species in the non-weighted developmental morphospace share many similarities: They all have similar directions of propagation and a slight parabolic shape (Fig. 8). The variation detected by the trajectory analysis (R function RRPP::trajectory.analysis) revealed significant differences in trajectory length (path distance), trajectory shape and trajectory slope between most species (File S1, summarized results below).

Members of the type A ontogeny have smaller magnitudes of shape change and different trajectory shapes than members from type B ontogeny (path distances A: 0.2005–0.0225; B: 0.2276–0.2778). For the assignment of the species to the two types see Table 3. Only pairwise differences of the path distance and the trajectory shape between type A and type B species are statistically significant (magnitude of shape change: 17/66 possible pairs; trajectory shape: 16/66 possible pairs).

Most species have statistically significant pairwise differences in trajectory slope (57/66 possible pairs). Species of all pairs with non-significant pairwise p-values are in the same ontogenetic group (i.e., both belong to type A or B). The trajectories of the nine pairs with non-significant differences in slope, have non-significant magnitudes and shapes also. However, none of the species that share a common slope have overlapping biostratigraphic ranges.

Weighted developmental morphospace occupation

The first three components of the PCA on the weighted shape space account for 94.8% (wPC₁ = 85.8%, wPC₂ = 7%, wPC₃ = 3.0%) of the total variation. In comparison to the regular PCA, the wPCA is—by definition—more strongly controlled by the centroid size of the configurations, which is mainly expressed by the domination of PC1.

Similar to the regular PCA plot (Fig. 8), the wPCA morphospace (Fig. 11) can be divided into three main parts: (1) The extremely depressed earliest whorls cover the lower left quadrant (low PC1 and PC2 values); (2) the center of the plot (PC1 equals 0, PC2 positive) is occupied by the more depressed whorls of juveniles and (3) adult whorls are associated with positive PC1 values. In contrast to the PCA, the wPCA reveals a more distinct separation of the type A and type B groups of adult whorls (see Table 3). Representatives of the more depressed type B clearly occupy the lower right quadrant (positive PC1 and negative PC2 values). This division into type A and B can also be seen in the mean shapes of the whorl 5.5 of the respective species (Fig. 10).

Ontogenetic patterns in Ceratitidae

The ontogenetic trajectories of the studied species comprise the biphasic development from strongly depressed to weakly depressed to compressed whorl profiles (Figs. 8–11). It is commonly accepted that sudden changes in ontogenetic allometry often mark the onset of sexual maturity (i.e., Kullmann & Scheuch, 1970; Klug, 2001; Klug et al., 2015a; Klug et al., 2015b).

The studied species can be divided into two main ontogenetic groups: Type (A) Truncated trajectories that are associated with depressed adult whorls; type (B) longer, complete trajectories that lead to a compressed adult whorl shape. The process of lengthening and shortening of the trajectories (i.e., related to changes in rate and timing of the development) account for the ontogenetic differentiation of the species in focus. This contrasts a previous traditional morphometric analysis by Bischof & Lehmann (2020) of ptychitids, which revealed that the spherocone-cadicone morphospace is much more distinct. The highly ontogenetically differentiated genus Ptychites directly differed through characteristic ontogenetic trajectories.

While precise temporal growth rates of ammonoids are unknown (Lécuyer & Bucher, 2006; Knauss & Yacobucci, 2014), a basic assumption herein was that the individual species have similar coiling rates (i.e., the individual species develop the same number of whorls in the course of their life). Modified rate/timing of shape change from any ancestor to any descendent within an evolutionary framework is called heterochrony (Zelditch, Swiderski & Sheets, 2012, p. 317). Between type A and type B species, interspecific variation of the species in focus arises from an acceleration, a special case of peramorphosis; (for discussion of this term, see Alberch et al., 1979; McNamara, 2012) that allows type B species to occupy an extended portion of the morphospace characterized by more compressed whorls. Therefore, the studied ceratitids do not primarily differ in shape, but rather in the timing of the development of individual shapes. Heterochrony as a mechanism in macroevolution is known to be a key driving factor in phenotypic diversification (e.g., Gould, 1977; Alberch et al., 1979; McKinney & McNamara, 1991; Gerber, Neige & Eble, 2007; Gerber, 2011; Korn et al., 2013; Knauss & Yacobucci, 2014). The quantification of patterns of morphologic disparity and the relationship between size and shape (i.e., heterochrony) will be the subject of future studies.

Anisian ammonoid diversity

It is widely agreed that ammonoid diversity reached its maximum during the Triassic period (House, 1993; Brayard et al., 2009; Whiteside & Ward, 2011). Thereby, the late Anisian ammonoid diversity peak was dominated by members of the family Ceratitidae (Brayard et al., 2009; supporting material Fig. S2). However, there are a growing number of studies critically questioning diversity peaks by arguing that—to some extent—the high diversity might be artificially inflated by taxonomic over-splitting (Forey et al., 2004; De Baets, Klug & Monnet, 2013; Knauss & Yacobucci, 2014).

The results obtained here do justice to the general opinion that ontogenetic trajectories can be a powerful tool to describe (e.g., Korn & Klug, 2007) and discriminate ammonoid species (e.g., Rieber, 1962; Bischof & Lehmann, 2020): The newly introduced methods succeeded in statistically discriminating the ontogenetic pathways of the pre-defined ceratitid species. Based on this analysis, the high diversity of the Anisian ammonoid assemblages of Nevada appears not to be artificially inflated and the alpha taxonomy is regarded to be adequate. However, the high morphological resemblances of the investigated species cannot be denied. Therefore, this study supports the main idea of McGowan (2004), McGowan (2005) and Brosse et al. (2013) that taxonomic diversity and morphological disparity need not necessarily be closely linked.

It is important to be aware of the fact that GMM carry no direct biological information. They help to understand if and how configurations differ, but not what the underlying mechanisms for their morphological development are. In the complex discoidal morphospace landmark-based approaches have proven to be useful to evaluate a priori defined taxonomic groups. Nevertheless, geometric morphometric methods cannot be considered as being a phylogenetic or taxonomic tool per se. But they certainly represent an improvement and valuable supplement to traditional methods.

Why it is worth the effort

There is no doubt that preparation and analysis of ontogenetic cross sections involves a lot of work (Korn, 2012). However, geometric morphometric methods (GMMs) open the door to a new world of objectified, statistically quantifiable descriptions. For example, in the case of the fauna described herein, conventional descriptions and traditional morphometric methods did not succeed to differentiate species adequately. Landmarks and semi-landmarks, however, make it possible to statistically quantify shape variations of entire morphologies (Neige, 1999) and allow the analysis of shape and size separately (Hammer & Harper, 2005).

The high resolution of the ontogenetic trajectories of the herein studied material was achieved owing to the accretionary planispiral growth of ammonoids with conservation of previous growth stages (Korn, 2012), which adds an intuitive, relative time-component to the ontogenetic analysis. Even though it is likely that small-scale ontogenetic changes are overlooked at a measurement density of one measurement per 180 degrees, it can be assumed that no major developmental steps were skipped (Tajika & Klug, 2020). Leaving out complete ontogenetic stages would most likely prevent the recognition of ontogenetic processes such as heterochrony. If, for example, only the earliest and latest stages of the ceratitid development were analyzed, representatives of type A and type B would differ fundamentally. The accretionary growth of many ammonoid conchs therefore not only adds an individual time component to the analysis, but more importantly ensures that no major developmental steps have been overlooked. This reinforces the general opinion that ontogenetic trajectories of ammonoids are a powerful tool to study evolutionary processes.