Analysing planktonic foraminiferal growth in three dimensions with foram3D: an R package for automated trait measurements from CT scans

. Foraminifera are one of the few taxa that preserve their entire ontogeny in their fossilised remains. Revealing this ontogeny through micro-computed tomography (CT) of fossil planktonic foraminifera has greatly improved our understanding of their life history and allows accurate quantiﬁcation of total shell volume, growth rates and developmental constraints throughout an individual’s life. Studies using CT scans currently mainly focus on chamber size, but the wealth of three-dimensional data generated by CT scans has the potential to reconstruct complete growth trajectories. Here we present an open-source R package to analyse growth in three-dimensional space. Using only the centroid xyz coordinates of every chamber, the functions determine the growth sequence and check that chambers are in the correct order. Once the order of growth has been veriﬁed, the functions calculate distances and angles between subsequent chambers, determine the total number of whorls and the number of chambers in the ﬁnal whorl at the time each chamber was built, and, for the ﬁrst time, quantify trochospirality. The applications of this package will enable repeatable analysis of large data sets and quantiﬁcation of key taxonomic traits and ultimately provide new insights into the effects of ontogeny on evolution.


Introduction
Foraminifera are one of the few taxa that preserve an entire ontogeny in their fossilised remains. This preservation enables reconstructions of ontogenetic trajectories in deep time, a feature of particular interest to studies investigating the influence of development on long-term evolutionary change (Brombacher et al., 2022a). Developmental plasticity can influence phenotype frequency in a population through a process called genetic accommodation (West-Eberhard, 2003). A new environmental cue will cause a plastic trait to be expressed in a novel way, and if this new phenotype has a positive effect on fitness, it will likely be selected for, increasing the frequency of both the phenotypic and genetic components (West-Eberhard, 2005. Developmental plasticity has been argued to both drive and inhibit evolutionary innovation, but the majority of research is based on theoretical models (Dewitt et al., 1998;Murren et al., 2015;Price et al., 2003) and/or modern populations (e.g. Pigli-ucci et al., 2006;Beldade et al., 2011;Moczek et al., 2011;Pfennig et al., 2010) generally limited to a handful of generations. Fossils contain information on both macroevolutionary transitions and microevolutionary change, but the lack of data on juvenile states makes reconstructing developmental trajectories difficult and so limits our ability to infer the role of developmental plasticity on macroevolution. Planktonic foraminifera could shed new light on the developmental drivers of evolutionary change. Freely distributed methodological tools would facilitate this contribution. Here we present a new, open-source R package that automatically analyses three-dimensional foraminiferal growth trajectories from micro-computed tomography (CT) scans. Our package enables fuller use of the incredible richness from x-ray CT and thus more comprehensive understanding of the role ontogeny plays in determining the size and shape of adult forms.
In the last decade, the application of micro-CT scanning to fossil planktonic foraminifera has greatly improved our un-Published by Copernicus Publications on behalf of The Micropalaeontological Society.
We present an open-source R package (available on GitHub, https://github.com/AniekeBrombacher/foram3D, last access: 31 August 2022) that automatically analyses planktonic foraminifera growth in three-dimensional space using the chamber centroid xyz coordinates (Fig. 1). The functions calculate distances and angles between centroids to arrange chambers in order of growth, calculate distances and angles between subsequent chambers, determine the total number of whorls and the number of chambers in the final whorl at the time each chamber was built, and, for the first time, quantify trochospirality. Included in the package are example data sets of Menardella limbata and Trilobatus sacculifer specimens to illustrate function usage for contrasting morphologies. The goal of this study and the intended applications of this package are envisaged to enable repeatable analysis of large data sets and ultimately provide new insights into the effects of developmental processes on the evolution of planktonic foraminifera.

Functions
The functions of the foram3D package calculate distances and angles between centre points of individual chambers. Chamber centres can be represented by centroids (the geometric centre of a three-dimensional object) or the centre of mass (the point where the entire mass of an object is concentrated). The examples in this study were reconstructed by filling in chambers on individual CT-scan image slices, cre-ating objects of uniform density (see Caromel et al., 2016Caromel et al., , 2017Kendall et al., 2020;Schmidt et al., 2013 for detailed methodology), for which centroid and centre of mass coordinates are identical (but note that these coordinates can differ for studies assigning different densities to chamber walls and cavities). No species-specific morphological assumptions are made. We illustrate the package functions using example Menardella limbata and Trilobatus sacculifer specimens.

Chamber ordering
The order of growth is essential for analyses of ontogeny. Current CT image analysis relies on the manual detection of chambers in individual slices, which are then put together to calculate chamber volume and centroid coordinates (e.g. Caromel et al., 2016Caromel et al., , 2017Schmidt et al., 2013;Kendall et al., 2020;Speijer et al., 2008;Brigulgio et al., 2011;Burke et al., 2020;Vanadzina and Schmidt, 2022). With chamber centroids plotted in three-dimensional space (see the code in the Supplement, lines 37-57 for interactive examples), foraminifera growth spirals are relatively easy to recognise, and chamber numbers can be assigned manually. However, ordering chambers this way is a time-consuming process that is prone to manual error, particularly for the smallest, earliest chambers. A repeatable automated process would speed up the process and reduce error.

Algorithm
Foraminifera grow in approximate logarithmic spirals with the distance between subsequent chambers increasing roughly exponentially (Signes et al., 1993). Starting from the earliest chambers, the algorithm finds the closest unassigned chamber to every chamber to determine the order of growth (see Fig. 2). The algorithm starts by finding the starting point of the spiral. The second chamber has always been found to be the smallest chamber (Brummer et al., 1986(Brummer et al., , 1987Huang, 1981;Sverdlove and Be, 1985), so chamber 2 is assigned first (Figs. 2a,3a). The first chamber can then be found by looking for the chamber closest to chamber 2 (Fig. 2b). The first and second chamber calcify together (Davis et al., 2020;Takagi et al., 2020) forming a flat plane as internal chamber wall, putting the centroid of the first chamber closer to the centroid of the second chamber than the third chamber. Chamber 3 is the next closest centroid to chamber 2 (Fig. 2c). From there on, when the first i chambers have been assigned, chamber i +1 is the closest remaining chamber to chamber i ( Fig. 2df). This step repeats until all chambers are assigned.
Occasionally, the proloculus is not present in fossil specimens, for example due to dissolution of the earliest chambers (Johnstone et al., 2010;Iwasaki et al., 2015). In that case, the algorithm takes the smallest chamber as the starting point. This version of the algorithm is less reliable than the version with the proloculus, as size does not always increase between subsequent chambers. For example, in Fig. 3a, if only chambers 15-22 were preserved, the algorithm would assign chamber 16 as the earliest chamber (rather than chamber 15). Generally, chamber size alone is not a reliable means to determine the order of growth. A number of species, such as the Menardella lineage, show reduced growth in the final chambers; small bullae on top of the final chamber are common, and occasionally unusually small chambers are found in random locations in the ontogeny (Duan et al., 2021). These deviations of the growth trajectory can have a profound effect on the resulting chamber order (Fig. 3b).

Usage
The function requires vectors with chamber x-, yand z-centroid coordinates and chamber volume (code in the Supplement, lines 63-89). Its default version assumes that the proloculus is present but can be changed to proloculus=FALSE if necessary. It returns a data frame with the original data ordered by the newly assigned chamber numbers. See below for an example of the function usage in R. See Sect. A1 in Appendix A for function output.

Checking the chamber order
Once chambers have been ordered (either by algorithm or manually), the R package contains a function that checks the chamber order and identifies unusual growth patterns to identify incorrectly ordered specimens. In foraminifera specimens with three or more chambers per whorl, the angle between subsequent chambers is ≥ 60 • . Therefore, sequences with angles smaller than 60 • could be an indication of incorrectly ordered chambers. Common reasons for unusual sequences include manual error, faulty data or disrupted growth. The algorithm flags up any specimens that need additional inspection.

Algorithm: angles
Angles between subsequent chambers are calculated using the chamber.angles() function (code in the Supplement, lines 97-104). It requires vectors with chamber x-, yand zcentroid coordinates and returns a vector with the angle every chamber makes with its previous two chambers (in degrees). Note that angles are generated from the third chamber onwards because the first and second chamber sit on a straight line.

Usage
The function requires vectors with chamber x-, yand zcentroid coordinates. It returns a vector with the angle each chamber makes relative to the previous two previous chambers (see the code in the Supplement, lines 97-104). Note that at least three chambers are required for angle determination; therefore, no angles can be calculated for chambers 1 and 2. Code shown below is for the exemplar M. limbata specimen. See Fig. 5 and Sect. A2 for function output. See the code in the Supplement, line 104, for the output for T. sacculifer.

Algorithm: chamber order check
To detect incorrectly ordered specimens, the algorithm calculates all angles α i that chamber i makes with the previous two chambers i − 1 and i − 2 (Fig. 4). Angle α 3 is discarded if the proloculus is known to be present: the first three chambers form a triangle, and α 3 therefore forms an angle smaller than 60 • with the proloculus and deuteroconch. In specimens where the first few chambers are absent, for example due to dissolution, all chamber angles are considered. Specimens with all angles larger than 60 • are labelled as correctly ordered (Fig. 6a). Any specimens with one or more angles smaller than 60 • are flagged to the user as potentially incorrectly ordered specimens.
Specimens with angles smaller than 60 • either represent true negatives with one or more incorrectly ordered chambers ( Fig. 6b) or specimens that experienced unusual growth resulting in sharp angles (Fig. 6c). Unusual growth can include a bulla or a smaller-than-expected chamber elsewhere in the sequence. Any manually ordered specimens flagged as incorrectly ordered can be run through the chamber-ordering algorithm to attempt to correct to sequence. If the sequence still returns sharp angles after ordering, the specimen needs to be inspected manually. In a preliminary analysis of 60 specimens, manual inspection was necessary for 5 % of the   specimens. In all cases unusual growth was the cause of the ordering issues.

Usage
The function requires vectors with chamber x-, yand z-centroid coordinates. The default option assumes that the proloculus is present, which can be changed to proloculus=FALSE if necessary. It returns a vector with either "correct order" if all angles are larger than 60 • or "check chamber order" if at least one chamber angle is smaller than 60 • . See the code in the Supplement, lines 67 and 74, for examples. See Sect. A3 for function output.  . Three examples of correct and incorrect chamber sequences using the exemplar Menardella limbata data released with the package. Specimen (a) is a correctly ordered specimen with all angles larger than 60 • . Specimen (b) is not ordered correctly and contains multiple angles much smaller than 60 • (red lines) as a result. Specimen (c) is correctly ordered but contains a sharp angle due to unusual growth of the penultimate chamber and is therefore flagged as an individual that needs closer inspection (note that the flagged angle appears larger than 60 • due to two-dimensional representation).

Number of chambers in the final whorl
The number of chambers in the final whorl is a frequently used diagnostic feature in taxonomy and can help distinguish between species and genera (e.g. Kennett and Srinivasan, 1983;Olsson et al., 1999;Pearson et al., 2006;Wade et al., 2018). Traditionally this number is determined from scanning electron microscope (SEM) or light microscope images and only describes the visible outer whorl. In contrast, CT scans also reveal all inner chambers. Using chamber centroid coordinates from CT data, we can determine the number of chambers in the final whorl at the time each chamber was built (Fig. 7). These new data could provide new insights into when ontogeny differences in whorl size between species and genera first arise. Manual determination of the number of chambers in the final whorl can be ambiguous, resulting in the inclusion of partial chambers in the final number (e.g. Globoconella inflata is described as having 3-3.5 chambers in the final whorl). These manual results are subjective and depend on visual volume determination of a partially hidden chamber. The function presented here uses distance between chambers to identify which chamber in the previous whorl is closest to the chamber under consideration. Therefore, the output is always in integers and no incomplete parts of chambers are considered.

Algorithm
To calculate the number of chambers in the final whorl at the time chamber i was built (cfw i ), the algorithm searches for the closest chamber in the whorl immediately below chamber i. This chamber marks the end of the previous whorl (see Fig. 8a). Depending on the chamber configuration, the closest chamber to chamber i can be either chamber i − 1 or the chamber in the whorl below. The algorithm selects the chamber with the shortest distance to chamber i that is not chamber i−1. If the closest non-neighbouring chamber to chamber i is chamber 1, chamber i is in the first and only whorl at time i (see Fig. 8b). In this case, cfw i = i. If the final chamber of the previous whorl is not the proloculus, cfw i = i − # final chamber of previous whorl (Fig. 8c, d).

Usage
The function requires vectors with chamber x-, yand zcentroid coordinates. It returns a two-column data frame with the closest non-neighbouring chamber to each chamber and the number of chambers in the final whorl at the time each chamber was built (see the code in the Supplement, lines 115-122). Results shown are for the exemplar M. limbata specimen. See Sect. A4 for function output. See the code in the Supplement, line 122 for the output for T. sacculifer.

Number of whorls
Taxonomic descriptions often include the number of whorls typically found in adult specimens (e.g. Kennett and Srinivasan, 1983;Olsson et al., 1999;Pearson et al., 2006;Wade et al., 2018). Three to four whorls are common in planktonic foraminifera, whereas benthic foraminifera, especially larger species, can have many more (Holbourn et al., 2013). Differences in the numbers of whorls among specimens can help species identification. A change in the number of whorls through time can point to an altered developmental history, such as paedomorphosis (reduced development) or peramorphosis (extended development).
To determine the total number of whorls, all chambers need to be visible. Determining the number of whorls manually can introduce bias similar to counting the number of chambers in the final whorl by eye. By letting an algorithm find the closest chamber in the whorl below individual chambers from CT scan data, the number of whorls can be determined objectively.

Algorithm description
The algorithm starts with the proloculus and makes its way up the logarithmic spire from there. The first whorl consists of all chambers whose closest non-neighbouring chamber is the proloculus (Fig. 9a). The second whorl starts with the first chamber whose closest chamber in the whorl below is not the proloculus (chamber 7 in Fig. 9a). Subsequently, whorl i + 1 starts with the first chamber whose closest chamber in the whorl below is in whorl i (Fig. 9b). The last whorl is often not complete (e.g. Pearson et al., 2006) (Fig. 9c). The algorithm checks if the last whorl is complete by comparing the number of chambers in the last whorl as determined algorith-mically (counted forward from the proloculus onwards) with the number of chambers in the final whorl at the time the final chamber was built (counted backwards, starting with the final chamber). If these numbers are not the same, all chambers in the last whorl are flagged as being part of an incomplete whorl. For example, the specimen in Fig. 9c had eight chambers in its final whorl at the time the final chamber was built, but only five of these were part of the actual last whorl.

Usage
The function requires vectors with chamber x-, yand zcentroid coordinates. It returns a two-column data frame with the whorl number of each chamber and whether the whorl is complete (code in the Supplement, lines 124-130). All whorls but the final whorl are complete by definition, whereas the final whorl is either complete or incomplete. See Sect. A5 for function output.

Coiling direction
Coiling direction is an important trait for biostratigraphy, with changes in the dominant coiling direction marking biostratigraphic horizons (e.g. Wade et al., 2011). For some species the coiling direction is the primary means of identification, such as in Neogloboquadrina pachyderma and N. incompta (Darling et al., 2006). Although the coiling direction of adult individuals can easily be determined by eye, the coiling of earlier, hidden ontogenetic stages is more difficult to analyse. The R package includes a function that determines the coiling direction throughout ontogeny. This could prove useful to determine the exact moment of coiling change in the rare planktonic foraminifera species that build whorl-enveloping chambers such as Orbulinoides beckmanni or more likely in streptospirally coiled benthic foraminifera species.

Algorithm
To determine the coiling direction of chamber i we compare the vertical direction of growth to the cross product of the vectors v and w between chambers i−2, i−1 and i. The cross product creates a vector perpendicular to both v and w, with the direction of v × w determined by the direction of travel from v to w (Fig. 10). Clockwise travel as viewed from the umbilical side (i.e. sinistral coiling) produces a cross-product vector pointing away from the vertical direction of growth, whereas anticlockwise travel as viewed from the umbilical side (i.e. dextral coiling) produces a vector pointing in the same direction as the direction of vertical growth (give or take minor deviations). Therefore, if the angle between v ×w and the direction of growth is smaller than 90 • , the coiling direction is anticlockwise or dextral, whereas an angle larger than 90 • indicates clockwise or sinistral coiling.

Usage
The function requires vectors with chamber x-, yand zcentroid coordinates. It returns a vector with the coiling direction of each chamber relative to the two previous chambers (see the code in the Supplement, lines 132-138). Note that at least three chambers are required for coiling direction determination; therefore, no coiling direction can be established for chambers 1 and 2. Also note that the earliest chambers in the juvenile stage grow nearly planispirally, the smallest deviations from which can result in an apparent change of coiling direction. See Sect. A6 for function output.

Trochospirality
The height of the trochospire is another trait used to distinguish among species (Fabbrini et al., 2021;Biolzi, 1991) and genera (Gradstein and Waskowska, 2021;Lipps, 1966;Pearson and Coxall, 2014). Specimens are qualitatively described as either planispiral (Fig. 11a) or having a low or high trochospire (Fig. 11b, c), but there is currently no direct way to quantify the height or angle of the trochospire. Quantifying the trochospire can help determine the exact moment of evolutionary divergence of two closely related species through statistical analyses (e.g. Pearson and Ezard, 2014). The height of the trochospire has also been noted to change through ontogeny: growth is nearly planispiral in the juvenile stage but increases as individuals mature (Caromel et al., 2016;Apthorpe, 2020;Kendall et al., 2020;Poole and Wade, 2019;Morard et al., 2019). The timing of change could be used to determine the transition from the juvenile to the neanic stage, which is currently only described qualitatively through manual inspection (Brummer et al., 1986(Brummer et al., , 1987.

Algorithm
Shells growing in approximate logarithmic spirals, such as ammonoids and gastropods, are traditionally described using "Raupian" parameters (Raup, 1966(Raup, , 1967) that quantify growth relative to the coiling axis. The translation rate or growth in the direction of the coiling axis for every ontogenetic step can theoretically be used to calculate trochospirality in foraminifera (Morard et al., 2019;Caromel et al., 2016). This method depends on correct identification of the coiling axis. Ammonoids and gastropods grow continuously and roughly isometrically, making the coiling axis easy to identify. Foraminifera, however, grow in discrete steps, and growth patterns change through ontogeny (Signes et al., 1993). Therefore, the coiling axis of a foraminifera specimen Figure 9. Determining the total number of whorls for the example Menardella limbata specimen. (a) The first whorl consists of all chambers whose closest chamber is the proloculus (black). Whorl 2 (grey) starts with the first chamber whose closest neighbour in the whorl below is not the proloculus. (b) Whorl i + 1 (grey) starts with the first chamber that has a chamber from whorl i (red) as its closest non-neighbouring chamber. (c) All three complete whorls (black, red, blue) and the fourth incomplete whorl (orange). vector w through chambers i −1 and i (blue) is clockwise as viewed from the umbilical side (i.e. sinistral coiling), the cross-product v × w (purple) points towards the proloculus in the opposite direction to the direction of growth. (b) If the direction of change from v to w is anticlockwise as viewed from the umbilical side (i.e. dextral coiling) v × w aligns with the direction of growth. Note that the orientation of the specimen does not influence the coiling direction outcome of the algorithm: if the same specimens were viewed from the spiral side, both the cross product and the direction of growth would flip by 180 • , maintaining their similar or opposite relationship necessary to determine the coiling direction. changes through ontogeny and is impossible to determine mathematically. Approximations of the coiling axis would need to be done by hand, potentially introducing bias and influencing any resulting calculations of trochospirality.
Alternatively, we can define the trochospirality of chamber i as the angle between chamber i and the plane defined by the previous three chambers (Fig. 11). This method does not rely on a coiling axis and only requires chamber centroid coordinates. Specimens following a perfect logarithmic spiral will have a constant trochospirality value throughout ontogeny.
However, natural growth is rarely perfect. The trochospirality of a specimen can be defined as the average trochospirality τ i of all chambers, and deviations in τ i can be used to assess developmental plasticity and ontogenetic constraints (Fig. 12).

Usage
The function requires vectors with chamber x-, yand zcentroid coordinates. It returns a vector with the trochospirality angle of each chamber, relative to the plane formed by the previous three chambers (see the code in the Supplement, lines 106-113). Note that at least three previous chambers are required to form a plane. Therefore, no trochospirality angle exists for the first three chambers. See Sect. A7 for function output.

Package applications and future work
For ease of use, all trait functions have been combined into a single function: foram.growth.3D(). This function requires vectors with chamber number, centroid coordinates and volumes and returns a data frame with the original data, plus angles between subsequent chambers, trochospirality, chambers in the final whorl at the time each chamber was built, the whorl number each chamber belongs to, whether the whorl is complete, coiling direction, and chamber order (see the code in the Supplement, lines 140-150). See Sect. A8 for function output.
with(Mlimbata, foram.growth.3D(n=Cham ber, x=Centroidx, y=Centroidy, z=Centroidz, proloculus=TRUE)) Plotting combinations of traits against each other can reveal ontogenetic trajectories and trait covariation patterns through ontogeny. The example data show stable chamber Figure 11. Theoretical examples of (a) planispiral growth and (b, c) low and high trochospires, respectively. Planes represent the growth plane of the entire specimen in the case of planispiral growth (a) or the plane defined by chambers i − 3, i − 2 and i − 1 for non-planispiral growth (b, c); τ i is the trochospire angle of chamber i relative to the plane defined by the previous three chambers. angles until around chamber 14, when angles increase with chamber number (Fig. 13, far left column). Trochospirality and growth rate decrease around the same time, with the number of chambers in the final whorl increasing several chambers later. The remaining columns in Fig. 13 show the correlations between all remaining trait combinations. Strong covariations between some traits are expected, such as whorl number and chamber volume (Fig. 13, bottom row, second from right), because chambers in earlier whorls are typically smaller. Others, such as the covariation between the number of chambers in the final whorl and chamber volume (Fig. 13, bottom row, third from right), are not immediately obvious and could point to a size threshold for whorl expansion. Similarly, chamber angles increase at the same time the rate of chamber growth starts to decrease, several chambers before the number of chambers in the final whorl increases, suggesting a switch away from isometric growth that leads to a reduced chamber-by-chamber growth rate and thus increasing numbers of chambers in the final whorl.

Prospects for future work
We present a new R package to automatically reconstruct foraminifera growth trajectories from chamber x, y and z centroid coordinates. The package functions arrange chambers in order of growth, calculate distances and angles between chambers, determine the total number of whorls and the number of chambers in the final whorl at the time each chamber was built, and, for the first time, quantify trochospi-rality. Applied to large numbers of specimens through multiple lineages and evolutionary transitions, the foram3D package functions can shed new light on ontogenetic trajectories and developmental constraints through time. They can be used to quantify ontogenetic stages, for example through changes in trochospirality in the juvenile and neanic stages or the increase in chambers in the final whorl for the neanic to adult transition. Additionally, they can be used to quantify species differences, for example through average trochospirality, angles between chambers in the final whorl and the total number of whorls. This will be particularly valuable for determining the exact moment of the origination of a new species, which is often difficult to identify by eye. Our package enables multivariate analyses of the incredible richness from x-ray CT data and could substantially improve our understanding of developmental constraints on planktonic foraminifera growth and form. Appendix A "NA" stands for "not available" and indicates an empty entry. "Correct" "Correct" "Correct" "Correct" "Correct" "Correct" "Correct" ## [8] "Correct" "Correct" "Correct" "Correct" "Correct" "Correct" "Correct" ## [15] "Correct" "Correct" "Correct" "Correct" "Correct" "Correct" "Correct" ## [22] "Correct" NA NA "sinistral" "sinistral" "sinistral" "dextral" ## [7] "dextral" "dextral" "dextral" "dextral" "dextral" "dextral" ## [13] "dextral" "dextral" "dextral" "dextral" "dextral" "dextral" ## [19] "dextral" "dextral" "dextral" "dextral"  Brombacher et al., 2022b), which contains all functions and data used in this paper. An example user guide to R scripts, containing information on how to download the package into R, generate interactive figures and use code is available in the Supplement.
Author contributions. AB designed the functions and wrote the code. AB, ASB, WZ and THGE tested and helped improve the functions and code. AB drafted the manuscript, and ASB, WZ and THGE provided comments on earlier drafts.

Competing interests.
The contact author has declared that none of the authors has any competing interests.
Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.