Towards a morphological metric of assemblage dynamics in the fossil record: a test case using planktonic foraminifera.

With a glance, even the novice naturalist can tell you something about the ecology of a given ecosystem. This is because the morphology of individuals reflects their evolutionary history and ecology, and imparts a distinct 'look' to communities--making it possible to immediately discern between deserts and forests, or coral reefs and abyssal plains. Once quantified, morphology can provide a common metric for characterizing communities across space and time and, if measured rapidly, serve as a powerful tool for quantifying biotic dynamics. Here, we present and test a new high-throughput approach for analysing community shape in the fossil record using semi-three-dimensional (3D) morphometrics from vertically stacked images (light microscopic or photogrammetric). We assess the potential informativeness of community morphology in a first analysis of the relationship between 3D morphology, ecology and phylogeny in 16 extant species of planktonic foraminifera--an abundant group in the marine fossil record--and in a preliminary comparison of four assemblages from the North Atlantic. In the species examined, phylogenetic relatedness was most closely correlated with ecology, with all three ecological traits examined (depth habitat, symbiont ecology and biogeography) showing significant phylogenetic signal. By contrast, morphological trees (based on 3D shape similarity) were relatively distantly related to both ecology and phylogeny. Although improvements are needed to realize the full utility of community morphometrics, our approach already provides robust volumetric measurements of assemblage size, a key ecological characteristic.

Introduction

Speciation and extinction are population-level processes with global effects on biodiversity. A species ceases to be when the last individual of the last population dies, and a new species arises when two previously connected populations become sufficiently isolated [1]. This being the case, assemblage dynamics should provide the most direct test of various regulators of biodiversity—be they environmental, biological or neutral—but this is almost never done in deep time. The vast majority of taxa simply do not have fossil records up to the task. Instead, questions of biodiversity dynamics and their drivers are typically addressed at the species level (or higher) in one of two ways [2]: (i) fitting models of diversification to modern and, rarely, fossil phylogenies (e.g. [3-7]) and (ii) assessing correlates of global diversity dynamics from fossil compilations and databases, like the Paleobiology Database (PBDB) (e.g. [8-11]). In all cases, the role of various regulators is inferred from their end-effect on phylogenetic structure or standing diversity, with varying degrees of theoretical robustness to the inference (as discussed in [6,12]). For those few taxa that do have excellent fossil records, like marine microfossils [13], assemblage-level studies of populations through time offer the exciting possibility of testing evolutionary mechanisms hypothesized from other data types (see [14,15] for a macrofossil example).

Unfortunately, even for those rare taxa with the spatial and temporal coverage needed to track population dynamics through geological time, there still remains a daunting data-collection problem. Measuring proxies of environmental and biological effects on numerous populations is extremely time intensive. As a result, most population-level work to-date has focused on environmental regulators of species abundance or range (e.g. [16]), or evolutionary dynamics within a lineage [17-19]. The tendency to investigate environmental regulators of population dynamics alone arises, in part, because environmental proxies are relatively quick and easy to measure, whereas high-throughput phenotypic methods have lagged for biological drivers. This study is aimed directly at this biological data-collection problem using a palaeontological model taxon: planktonic foraminifera.

Planktonic foraminifera are marine protists with calcium carbonate shells (technically, ‘tests’). Found throughout the global ocean today and abundantly for (roughly) the past 150 Myr, their tests rain to the sea floor on death and, in many regions, are preserved in near-continuously accumulating deposits [20,21]. Their abundance [20], the potential to measure multiple environmental proxies directly from their test geochemistry [21], and the recent completion of a phylogeny for Cenozoic macroperforate species (the major clade of planktonic foraminifera) [22] all contribute to the growing utility of planktonic foraminifera for understanding macroevolution [18,19,23-25]. Because planktonic foraminifera are a central tool in the field of palaeoceanography (the study of ancient oceans) [20,21,26], detailed records of environmental conditions already exist in many locations and time intervals over the past 66 Myr (e.g. [27-29]). Comparable biotic data are often lacking, however, and this disconnect remains a conspicuous hindrance to our ability to understand the feedbacks between biotic and abiotic processes in macroecology and macroevolution.

Here, we develop the methods for, and provide an initial view of the utility of, high-throughput semi-three-dimensional (semi-3D) geometric morphometrics as a means for rapidly tracking biotic dynamics through time. This approach rests on the assumption that the morphology of individuals, including their body size, reflects some combination of shared evolutionary history, functional ecology and individual variation. Several lines of evidence suggest that this may be the case in planktonic foraminifera. Importantly, planktonic foraminifera are well known for exhibiting iterative evolution of gross morphology. In multiple cases, complex morphologies including flat discoidal or peaked pyramidal morphologies, finger-like chambers and sharp edges (known as keels) have independently evolved from simple, globular ancestors [30-32]. Although functional morphology is poorly understood (see discussions in [33,34]), this ubiquity of convergent evolution in planktonic foraminifera and the correlation, in some cases, of morphology and life history support the inference of a functional role for gross morphology. If this is the case, then the morphological similarity of communities might provide a measure of functional similarity—an approach sometimes called ‘ecometrics’ and related to the burgeoning field of functional trait ecology [35].

That said, the relationship between diversity, functional diversity and morphological diversity is not straightforward (e.g. [36,37]), and morphological measures of community dynamics would necessarily remain just that—morphological measures—without thoughtful exploration and calibration. Even so, for deep time studies, community morphology provides a particularly promising means of assessing biotic dynamics because morphology provides a common ruler to compare across assemblages with entirely different species compositions [35,38]. For planktonic foraminifera, and many other groups, direct measures of morphology solve a second problem related to the common occurrence of morphologically intermediate individuals that lie between named taxa (for examples of morphological variation, see: [18,19,39]). Morphologically intermediate individuals can provide direct evidence of evolution in action, but their importance and implications are missed when morphologically intermediate taxa are shoehorned into named species-categories.

In planktonic foraminifera, it is clear that shared evolutionary history and factors unique to individuals influence morphology. Planktonic foraminiferal genera are often readily recognizable by shared, derived gross morphological characteristics, providing support for the inference that morphology must partially reflect the evolutionary relatedness of taxa. In some cases, speciation or evolutionary transitions have occurred across habitat types with relatively minor morphological change. Such instances include depth parapatry within a pseudo-cryptic species [40], abrupt ecological change within a gradual morphological series [41] and the occurrence of deep-water morphologies in shallow water habitats [42,43], and vice versa (as discussed in [33]). Finally, at an individual level, the availability of various resources (like light, food, temperature and oxygen) can profoundly influence adult size, shape and wall structure (e.g. thickness and porosity) [44-47].

The intent of this study is to lay the groundwork for the application of community morphology as a measure of population dynamics in planktonic foraminifera, although the underlying methods are fully applicable, and currently being used in-house, in macrofossils as well. To this end, we have:

— developed a pipeline to extract 3D data from light images (semi-3D morphometrics),

— compared morphological space represented by meshes from 16 species extracted using slow, but highly resolved, computed tomography (CT) (full-3D) and our new method (semi-3D),

— examined the relationship between morphology (full- and semi-3D), ecology and phylogeny in those 16 species,

— demonstrated the utility of our method for rapidly collecting traits like volume and surface area,

— investigated the trade-offs in data density and quality in 3D morphometrics and

— explored the potential of assemblage-wide ecometrics with four modern planktonic foraminiferal assemblages in the North Atlantic.

This work provides a new set of image-processing programs for extracting semi-3D data from light images, while highlighting the potential (and problems) of our approach, and provides a first exploration of the relationship between gross morphology (as measured with geometric morphometrics), ecology and phylogeny in 16 species of modern planktonic foraminifera.

Material and methods

Specimen sources

Two fundamentally different types of 3D data are used in this study: full-3D meshes from X-ray CT (CT scans) and semi-3D meshes from reflected light microscopy. The two data types have complementary strengths and weaknesses. CT captures the full 3D shape of planktonic foraminifera, including internal structure (although it should be noted that our subsequent analyses use only exterior shape from the CT scans), but is relatively slow to collect. By contrast, our light microscopic method (as detailed below) is very fast, but only captures exterior 3D shape from a single viewpoint (hence, semi-3D). In this study, we use the CT scans as a reference point to test the relationship between 3D morphology, ecology and phylogeny, and to test the relative information captured with the rapid semi-3D methods that we develop and introduce here.

Thirty-nine full-3D meshes from X-ray CT were obtained from Tohoku University Museum's e-Foram Stock database [48], representing 19 extant species of planktonic foraminifera (electronic supplementary material, table S1). The Tohoku University specimens were imaged with high-resolution X-ray CT scans (generated using ScanXmate-E090; Comscantecno Corporation) at 5 µm resolution [48]. We refer, hereafter, to these complete specimens as the Tohoku University specimens or as full-3D data.

We generated new, semi-3D meshes of 281 individuals identified to species level, representing 24 species of extant planktonic foraminifera. The complete list of all specimens used in this study (including database, catalogue number and species) is presented in electronic supplementary material, table S1. Five coretop locations in the North Atlantic were used as focal sites for this study to maximize taxonomic diversity/disparity at a community level and to explore the size distributions obtained from two-dimensional (2D) versus 3D imaging (figure 1). An additional 116 planktonic foraminifera from these sites, not identified to the species level, were included for community comparison purposes. The sites are: KC 78 (5°16′01″ N and 44°07′59″ W), CH 82-21 (43°29′17″ N and 29°49′48″ W), EW 93-03-04 (64°43′ N and 28°55′ W), VM 20-248 (33°30′ N and 64°24′ W) and AII 42-2-2 (18°01′59″ N and 24°27′ W). An example slide image for each site is available via the Yale University Peabody Museum of Natural History (YPM) online database (KE EMu) via the listed catalogue numbers and the Division of Invertebrate Paleontology Portal (http://peabody.yale.edu/collections/search-collections?ip). Full raw images are available upon request, as no public repository exists, free of charge, for large image files. YPM catalogue numbers for the five focal sites are: KC 78 (IP.307630–IP.307634), CH 82-21 (IP.307625–IP.307628), EW 93-03-04 (IP.307636–IP.307640), VM 20-248 (IP.307715) and AII 42-2-2 (IP.307647–IP.307651). To improve species coverage, additional exemplar individuals were included from 12 species: Truncorotalia crassaformis (IP.307720), Pulleniatina obliquiloculata (IP.307727), Truncorotalia truncatulinoides (IP.307733), Neogloboquadrina dutertrei (IP.307862), Globigerinella siphonifera (IP.307863), Globorotalia tumida (IP.307749), Menardella menardii (IP.307754), Menardii fimbriata (IP.307754), Sphaeroidinella dehiscens (IP.307757), Hirsutella hirsuta (IP.307760), Globigerinoides conglobatus (IP.307763, IP.307764), Globoconella inflata (IP.307767), Candeina nitida (IP.307772, IP.307773) and Globigerinella calida (IP.307865). Note that although only 281 specimens were identified to species level for morphological analyses, a total of 9681 individual planktonic foraminifera were used to compare body size distributions across four focal sites (no. of individuals by site: 1768 from KC 78; 2879 from CH 82-21; 3034 from EW 93-03-04; and 2000 from AII 42-2-2).

Figure 1.

Map of the five Atlantic coretops sites used in this study.

Specimen completeness and species identification were conducted by eye by PMH and LEE. Species were identified following the naming scheme of Aze et al. [22], and taxonomic concepts of Kennett & Srinivasan where applicable [49], to facilitate direct comparisons with the macroperforate phylogeny. Exceptions to the Aze et al. [22] species-naming scheme were as follows:

(i)Truncorotalia: All Truncorotalia were identified as either T. truncatulinoides or T. crassaformis. Aze et al. [22] recognizes six extant Truncorotalia. The first two (T. crassaformis and T. oceanica) are allied with the T. crassaformis complex and the remaining four (T. cavernula, T.excelsa, T. pachytheca and T. truncatulinoides) with the T. truncatulinoides complex, with morphological and genetic definitions varying among authors (e.g. [22,50-53]).

(ii)Pulleniatina: All Pulleniatina were identified as P. obliquiloculata, ignoring the possibility of P. finalis. This was a practical decision: from the common imaging angle (umbilical), these taxa are not readily distinguished.

(iii)Globigerinoides triloba: All G. triloba (morphospecies concept) were classified as Trilobatus sacculifer following genetic evidence for a single modern species in this morphologically variable taxa [54,55].

(iv)Neogloboquadrina: We recognized N. incompta, in addition to the N. pachyderma and N. dutertrei of Aze et al. [22], given widespread support for the genetic separation of this readily identified taxa [56,57].

In short, we generally favoured a ‘clumped’ taxonomic approach in naming (excepting N. incompta), allowing morphological variation to highlight differences among closely related lineages. Planktonic foraminifera morphospecies commonly harbour a few (pseudo-)cryptic genetic species [58,59], but the degree of splitting still varies among authors and is still being resolved [57].

We analyse the relationship between morphology, phylogeny and ecology, in the 15 species of macroperforate foraminifera in common between the full- and semi-3D datasets. These 15 species are listed in table 1, along with their ecological characteristics. The 15 macroperforate species and one microperforate species, C. nitida, were also used to assess the relative morphological information contained in full- and semi-3D meshes.

Table 1.

Ecological traits of overlapping macroperforate planktonic foraminifer species.

scientific name	symbiont type	habitat depth	geographical range	refs
Globigerina bulloides	none	mixed layer	mid-latitudes	[60–66]
Globigerinella siphonifera	chrysophytes	mixed layer/thermocline	low–mid latitudes	[40,60,63,67–71]
Globigerinoides conglobatus	dinoflagellates	mixed layer	low latitudes	[60,64,68,69]
Globigerinoides ruber	dinoflagellates	mixed layer	low latitudes	[21,60,64,68,69,71–73]
Globoconella inflata	chrysophytes	thermocline	low–high latitudes	[21,22,60,68,69,74]
Globorotalia tumida	none	thermocline/sub-thermocline	low latitudes	[21,22,60,63,68,69,73,75,76]
Hirsutella hirsuta	none	thermocline/sub-thermocline	low–mid latitudes	[21,60,68,69,74]
Menardella menardii	chrysophytes	thermocline	low latitudes	[21,60,63,68,71,74]
Neogloboquadrina dutertrei	chrysophytes	mixed layer/thermocline	low latitudes	[22,60,64,68,69,71,74,77,78]
Neogloboquadrina pachyderma	none	mixed layer/thermocline	low–high latitudes	[21,60,68–70,74,79]
Pulleniatina obliquiloculata	chrysophytes	mixed layer/thermocline	low latitudes	[21,60,68,69,71,74,79]
Sphaeroidinella dehiscens	dinoflagellates	thermocline	low latitudes	[22,60,68,80]
Trilobatus sacculifer	dinoflagellates	mixed layer	low latitudes	[21,60,64,68,69,71–73]
Truncorotalia crassaformis	none	sub-thermocline	low–mid latitudes	[22,68,69,71,79]
Truncorotalia truncatulinoides	none	sub-thermocline	low latitudes	[21,22,60,62,68,69,71,74,79]

Slide preparation and imaging for semi-3D data

Semi-3D data were collected from five sites (figure 1) and our species exemplar slide collection. For each of the five sites, a micropalaeontological split of the greater than 150 µm fraction was taken with a target sample size of 5000 individuals per site. Splits were scattered, oriented to an umbilical view, and glued to plain black micropalaeontogical slides. Approximately four slides were used per site to accommodate the roughly 5000 individuals in the split. Slides were imaged on a Leica Microsystems DM6000M compound microscope with transmitted light, a 5× objective and a 5-megapixal Leica DFC450 digital camera, using a 64-bit beta version of the controlling Surveyor Software. Each slide was scanned in a series of tiles in the x- and y-dimensions to cover the full length and width of the slide (using the automated stage). For each x–y tile, the automatic drive focus took a series of images at different heights (with a prescribed z-step) to capture the full-depth dimension of the fossils.

Each slide scan was saved as 32 bigTiff files, with the first 31 files capturing a single depth slice of the scan (x- and y-tiles composited). The 32nd bigTiff file is a 2D extended depth of focus image of the slide. A z-step size of 31.1 µm was used to account for the depth of focus of the 5× objective. This step size sets the limit of our ability to resolve shape in the z-dimension. Pixel size in the x- and y-dimensions was 0.975 µm. Microscope settings were the same for the exemplar slide scanning, with the only difference being the number of individuals scanned per slide (several individuals rather than thousands of individuals).

The reproducibility of the 3D-mesh extraction pipeline to variation in the orientation of mounted specimens and imagining problems (e.g. image tiling, glare) was tested using a single representative of each of the following eight species: Trilobatus sacculifer (IP.307747), G. tumida (IP.307749), P. obliquiloculata (IP.307751), N. dutertrei (IP.307753), Orbulina universa (IP.307756), Globigerinoides ruber (IP.307758), T. truncatulinoides (IP.307762) and G. siphonifera (IP.307765). Each individual was (re-)mounted five times and (re-)imaged in five different settings (per mount) with varying field of views and white balance/shade correction settings. The five imaging tests conducted per mount included: (i) a vertical image seam along specimen (e.g. specimen at the joint of a left and right image tile); (ii) a horizontal image seam along specimen (e.g. specimen at the joint of a upper and lower image tile); (iii) a centred specimen (single image tile) with the same white balance and shade correction for tests #1–3; (iv) a centred specimen (single image tile) with a readjusted white balance/shade correction; and (v) a centred specimen (single image tile) with a readjusted white balance/shade correction. In total, 25 images (five mounts per specimen and five imaging tests per mount) of the same individual were included for each species listed above, resulting in a final set of 200 reproducibility test ‘individuals’.

Preliminary image processing for semi-3D data

The segment (v. 1.10) and focus (v. 1.10) functions from the image-processing module of the AutoMorph software package (current software available on GitHub at https://github.com/HullLab) were used to extract individual objects from the scanned slides. segment chops slide scans up into individual objects and focus generates an extended depth of focus image for each object.

For segment, a black/white thresholding value of 0.18 was used for all slides except for CH 82-21 (threshold = 0.14), KC 78 (threshold = 0.20) and EW 93-03-04 (threshold = 0.20). The absolute thresholding value is unimportant for the analyses that follow, and are varied by slide to optimally segment out all foraminifera (common errors include segmentation of background glare, edge clipping of transparent or darkened individuals, etc.). The size range filter for a valid object was set to 125–2000 µm (width).

focus was run using the Zerene Stacker software [81] to generate a best extended depth of focus image (EDF) per object. Zerene Stacker settings included brightness correction between frames, automatic order of images (e.g. focus stacking begins on image with the narrowest field of view), default values of the estimation radius (10) and smoothing radius (5), a contrast threshold of 25%, and a grit suppression algorithm to reduce noise and pixellation during image stacking.

run2dmorph (v. 1.07, also available on GitHub), another module of the AutoMorph package, was run on the focused images to extract 2D morphology (e.g. outlines and corresponding coordinates) and shape parameters (e.g. major and minor axes length, enclosed area, eccentricity, rugosity, perimeter and aspect ratio) for body size analyses. All three processes (segment, focus, run2dmorph) were performed on the Department of Geology and Geophysics Tide server at Yale University.

Extraction of semi-3D mesh

A new 3D mesh extraction module, run3dmorph, was written by AYH as part of the AutoMorph software package to extract semi-3D meshes and object volumes from the objects (in this case, planktonic foraminifera) after preliminary image processing via segment and focus. A beta version of run3dmorph is available on GitHub at this time, with early adopters encouraged to check for updates (https://github.com/HullLab).

The generalized pipeline for the 3D mesh extraction is shown in figure 2. First, a greyscale EDF image and height map is generated for each individual using the Stack Focuser plugin [82] for ImageJ [83,84] and FIJI [85]; an 11 × 11 pixel kernel was used for the height map generation (figure 2b). The 3D mesh is then extracted from a cleaned height map (figure 2c,d) using a series of custom MATLAB scripts, which require MATLAB version 2015b or above [86].

Figure 2.

Visual pipeline illustrating the steps involved in 3D mesh extraction using the run3dmorph software. (a) Z-stacks of each individual object (taken at varying focal planes of known height above the object) are processed using the Stack Focuser plug-in for ImageJ/FIJI, resulting in a focused image of the object and a height map (built using an 11 × 11 pixel kernel size). (b) The focused image and height map are rescaled such that each pixel has a height and width of 1 µm, and the 2D outline of the focused image is extracted using the run2dmorph software (see https://github.com/Hull-Lab). Each pixel of the binary 2D outline image is then multiplied against the corresponding pixel in the height map (element-wise multiplication). This effectively deletes background noise and results in a cleaned-up height map (c). The greyscale value of each pixel in the height map is then used, in conjunction with the distance between each z-stack slice, to back-calculate the real-world height of each pixel and generate an unfiltered 3D mesh (d). High and low outlier noise is then filtered from the 3D mesh using a custom neighbourhood pixel-averaging algorithm (using a user-defined n × n pixel kernel, where n is a positive odd integer). Vertex and face coordinates are then extracted from the cleaned 3D mesh and outputted in standard 3D ASCII formats (OBJ and OFF).

More specifically, the semi-3D mesh extraction pipeline is as follows: First, all images are rescaled such that pixels are squared and 1 pixel = 1 µm (though any base unit can be used), using the conversion factor generated by the microscope calibration (0.975 µm/pixel for both the x- and y-dimensions in all cases here). The 2D outline from the run2dmorph module of the AutoMorph software package [55] is then used to exclude the background of the height map (figure 2b). This effectively deletes all background noise before the 3D mesh extraction step (figure 2c). The presence of a pronounced aperture (the opening in the final chamber) in some foraminiferal species resulted in errors in 3D mesh extraction owing to Stack Focuser's inability to capture aperture depth with fidelity. In most cases, pronounced apertures resulted in large noisy spikes in the final 3D mesh (figure 3a,b). To mask apertures, run3dmorph calls and runs an adjusted version of run2dmorph that skips the hole-filling step (figure 3c), thus allowing the aperture to be identified and excluded along with the background (figure 3d–f).

Figure 3.

Correction of artefacts arising from foraminifer apertures during 3D mesh extraction. (a) Example of aperture artefact, which manifests as a large peak of noise. (b) Lateral view of the aperture artefact. (c) Binary outline of the object with aperture excluded, as outputted by an adjusted version of the run2dmorph software. This binary image is then used, via element-wise multiplication, to remove the aperture during initial mesh extraction (d). The aperture height is then artificially set to the lowest height value in the mesh, as illustrated in top (e) and bottom view (f).

The cleaned height map (figure 2c) is then used to extract the height of each pixel based on the value of each pixel (between 0 and 255) in the height map, the number of slices in each z-stack (in our case, 31 z-slices) and the distance between each z-stack slice (31.1 µm), where the extracted height of each pixel was equal to:where H is the greyscale value of a pixel in the height map, s is the number of slices in the z-stack for the object in question, and Z is the distance between each slice in µm. This results in a matrix A with dimensions x × y, where x is equal to the width of the object and y is equal to the height of the object, and where each entry a corresponds to the actual height, in µm, of the pixel located at i,j. The accuracy of this approach in extracting heights from planktonic foraminifera was checked using the spherical species Orbulina universa. For the seven specimens examined, the extracted height was within 7.67% of the major and minor axes length.

A semi-3D mesh is then extracted from the scaled height map whereby every non-zero pixel of the height map is given a point in x, y, z coordinate space (i.e. x = horizontal, y = vertical, z = height). The extracted 3D mesh is then passed through a custom sliding neighborhood filter to remove outlier noise (figure 2d). This filter processes a kernel of size n × n pixels (where n is an odd positive integer). For each pixel, the filter calculates the upper and lower quartile ranges of all the pixel values encompassed in the neighbourhood (i.e. the n × n kernel). If the focal pixel falls outside the inner quartile range (i.e. below 25% or above 75%), it is considered an outlier and replaced with the mean value of all the pixels in the neighbourhood. This filter has the effect of removing both high and low outliers that result from noise created during height map generation, and also of smoothing the surface of the final mesh. For our specimens, we used n = 45 after testing several values of n to optimize processing time, noise deletion effectiveness and smoothing. Larger values of n are required to effectively delete larger patches of noise, but result in significantly increased computational resource requirements, and may also result in over-smoothing of the final mesh. The optimal balance between the kernel size used for height map generation and the kernel size used for outlier filtering varies between objects/sample sets, and must be determined through testing for each particular dataset.

Once the height map is filtered, the number of pixels present on each z-level (i.e. height) is counted. If the number of pixels in a given z-level is smaller than 1% of the total number of pixels in the object, that z-level is removed. This step has the effect of deleting any background noise along the edge of the object that may have been missed by the previous outlier filtering step. For our foraminifera, we also removed the bottom-most z-level for every object, as the majority of the meshes retained a rim of background around the object, thus obscuring the outline of the shape. Finally, all objects with apertures are then given an aperture depth equal to the lowest height in the overall object (figure 3f). Because the semi-3D meshes consist only of the visible upper portion of the foraminifera, this method results in apertures terminating approximately in the centre of the foraminifera.

Once fully processed, the semi-3D mesh is then extracted as a series of vertices and faces using the pointCloud2mesh function from the geom3d package [85] and saved in both Wavefront OBJ and Object File Format (OFF) format. The 3D mesh can be downsampled when saving (to minimize file sizes) but, for our purposes, the full mesh was retained (i.e. no downsampling was conducted). A CSV file containing the raw x-, y- and z-coordinates is also saved at this step. For quality control, run3dmorph can also output 3D PDFs of the extracted meshes, and uses the u3d_pre [87] function, a modified version of save_idtf, IDTFConverter (all from the mesh2pdf package; [88,89]) and the LaTeX package media9 (v.0.60; [90]) to do so.

In addition to semi-3D meshes, run3dmorph also estimates the volume and surface area of all objects, and saves an additional CSV file of these values. The volume and surface area of the extracted semi-3D hull are calculated exactly (by summing up the heights (or areas) of every pixel), while the volume and surface area of the bottom, un-imaged half are estimated. Three estimations of the bottom half shape are made in order to bound the volumetric and surface area uncertainty arising from the lack of direct measurement. They include base shapes of an irregular cone (figure 4a), an irregular cylinder (figure 4b) or a spheroidal dome (figure 4c). The surface area and volume of the irregular cylinder are calculated as P2D × H + A2D and A2D × H, respectively, where P2D = the length of the perimeter of the 2D outline, H = height and A2D = the area enclosed in the 2D outline. Height (H) above the background is calculated as the distance between the lowest image plane (e.g. the slide background) and the lowest z-level of the extracted semi-3D mesh. The 2D parameters, perimeter length and 2D area are outputs of run2dmorph.

Figure 4.

Illustration of idealized base shapes, used by run3dmorph for estimating surface area and volume of the complete object. (a) Irregular cone base; (b) irregular cylinder base; (c) spheroidal dome base.

The volume of the spheroidal dome is calculated as 1/2 * 3/4 * πs, or one-half of an ellipsoid where s = semi-axis X, s = semi-axis Y and s = semi-axis Z. For the dome, semi-axis X is equal to the length of the major axis of the 2D outline and semi-axis Y is equal to the length the minor axis of the 2D outline (figure 4c). Semi-axis Z is equal to the height H. The surface area of the spheroidal dome is estimated using Thomsen's formula, where k = 1.6:

Although the volume of the irregular cone can be calculated as one-third of the volume of the irregular cylinder with the same base, the surface area of the irregular cone cannot be calculated exactly. We estimate the surface area of the irregular conical surface using 100 perimeter coordinates extracted from the 2D outline, the height (H) and centroid of the cone, and the angle of inclination from each perimeter point and the centroid. The centroid of the 2D outline is determined using the MATLAB regionprops function, and then the angle of inclination (φ) between the horizontal and the line formed between the centroid and each perimeter coordinate is calculated, such that φ increases monotonically along the perimeter from 0 to 2π. The corresponding Euclidean distance between each perimeter coordinate and the centroid is also calculated. Then, using the midpoint integration rule, the generalized conical surface area is estimated by summing the areas of the trapeziums formed between each adjacent perimeter coordinate and a point of height H directly above the centroid.

The bottom estimates (conical, cylindrical and spheroidal) are used because they represent the full range of possible background shapes—from perfectly domed in the spherical Orbulina universa, to conical in H. hirsuta, to flat or filling-in taxa like T. truncatulinoides and N. dutertrei. They also span the theoretical maximum (cylindrical) and minimum (conical) possible volumes (or surface areas) of the back-half and thus allow for estimating the range of possible volumetric uncertainty introduced by assuming the shape of the back-half of the object. This uncertainty is calculated aswhere Econ = estimate of volume (or surface area) of the object with a conical back, Ecyl = estimate of volume (or surface area) of the object with a cylindrical back and Edom = estimate of volume (or surface area) of the object with a spheroidal back.

Mesh alignment and landmark placement for semi- and full-3D morphometrics

In this study, we use 3D semi-landmark geometric morphometrics to assess the morphology of species and communities. To use these analytical approaches, we first had to convert the semi- and full-3D meshes (from our specimens and the Tohoku University specimens) to landmarks. For a visual comparison of these two data types, see figure 5. Landmarks were placed using Boyer et al.'s [91] automated alignment and shape comparison algorithm, as implemented in the R (v. 3.0.2; [92]) package auto3dgm [93] and the parallelized cluster version of the algorithm implemented in MATLAB as PuenteAlignment [94]. Analyses using auto3dgm were conducted on the Tide server, whereas analyses using PuenteAlignment were conducted on the Yale High Performance Computing Omega cluster.

Figure 5.

Visual comparison of data types: Tohoku University 3D specimens from CT and semi-3D half-hulls extracted using the run3dmorph software on stacked microscopic images. Three examples specimens are shown: Trilobatus sacculifer, Truncorotalia truncatulinoides and Neogloboquadrina dutertrei. The Tohoku University specimens were digitized using X-ray CT at 5 µm resolution. Run3dmorph-extracted 3D-meshes are shown next to their corresponding focused 2D-image.

Two batches of landmark placement analyses were run for the semi- and full-3D specimens respectively. For the semi-3D specimens, landmarks were optimized using the parallelized PuenteAlignment for 597 objects that comprise the 281 species-identified individuals from the five Atlantic coretops and selected exemplar species, the 200 reproducibility tests ‘individuals’, and an additional 116 complete and well-imaged individuals (without species-level identification) from the five Atlantic coretops, examined by eye to ensure proper mesh extraction. For the full-3D analysis, landmarks were optimized using auto3dgm for the 39 Tohoku University specimens. For both the semi- and full-3D analyses, 256 landmarks per object were placed and used in the subsequent construction of morphospaces.

Morphospace construction

Four morphospaces were constructed to consider the key questions of the study: What is the relationship between morphology, ecology and evolutionary history in extant planktonic foraminifera? And, can semi-3D morphometrics capture community dynamics? To get at these issues, we first constructed two morphospaces for the complete set of semi- and full-3D data, respectively (comprising a total of 587 semi-3D individuals and 39 full-3D individuals). In order to directly compare morphospaces between the semi- and full-3D approaches, we also constructed two morphospaces, one for semi-3D and one for full-3D, including only the 16 overlapping species between the two datasets (e.g. the 15 macroperforate species listed in table 1 and the microperforate C. nitida). The pruned semi- and full-3D datasets contained 421 and 34 individuals, respectively.

In each case, morphospace was constructed using the Geomorph R package [95]. We first used the gpagen function to conduct a Generalized Procrustes Alignment (GPA) of the auto3dgm-/PuenteAlignment-extracted landmarks, with all landmarks set as sliding surface semi-landmarks. The reason for using semi-landmarks is twofold: first, as the Boyer et al. algorithm does not assume homology between landmarks [91], the resulting landmarks are not true geometric morphometric landmarks and should be allowed to move to minimize the potential distance between objects in morphospace. Second, Gonzalez et al. [96] found curvature in shape space when using the landmarks outputted by the Boyer et al. [91] algorithm as traditional landmarks, a pattern that we also observed in the planktonic foraminiferal shape spaces. With semi-landmarks, shape space is uncurved, supporting the use of semi-landmarks for a more accurate characterization of shape space. After a GPA is conducted on the semi-landmarks, a principal component analysis (PCA) was then conducted using the plotTangentSpace function in the Geomorph package. The resulting principal components capture the major axes of variation in the full- and semi-3D morphospaces for planktonic foraminifera.

Hierarchical clustering: morphology and ecology

Relationships among planktonic foraminifera in shape space (and subsequently, ecological space) were considered using hierarchical clustering with the hclust function in R. Hierarchical clustering was carried out on the complete full- and semi-3D morphological datasets (containing 39 and 587 individuals, respectively), and on the pruned datasets containing only the 16 species present across both datasets (table 1; containing 34 and 421 individuals, respectively).

For the morphological clustering, principal component (PC) scores were first averaged within species along each PC. A Euclidean distance matrix was then calculated for the species-averaged PC scores. This distance matrix was then used for hierarchical clustering using the Ward, single, complete, average and McQuitty linkage agglomeration methods in hclust. A 50% majority-rule consensus tree was then built from the resulting dendrograms using the consensus function from the ape (v. 3.0-10; [97]) R package in order to identify stable clusters. The final consensus dendrogram depicts consensus linkages between each species in the various analyses.

For the ecological clustering, ecological characteristics (symbiont type, habitat depth and geographical range) were taken from Ezard et al. [98] for 15 of the 16 overlapping species (C. nitida, the only non-macroperforate species, was excluded). Jaccard distances were then calculated between each species using the vegdist method from the vegan (v. 2.2-1; [99]) R package. Using the Jaccard distance matrix, hierarchical clustering and consensus dendrogram generation were conducted as described above.

Linear discriminant analysis for semi-3D data

One of our key findings (discussed in detail below) is that species overlap to a much greater degree in semi-3D morphospace than in full-3D morphospace. Reasoning that the species, identified by PMH based on morphology, are, in fact, morphologically distinct, we conducted a linear discriminant analysis (LDA) on the principal component coordinates prior to conducting hierarchical clustering for the semi-3D dataset. This was done to bring the variance that is useful in distinguishing species clusters to the forefront of the semi-3D analysis, thereby minimizing the potential noise introduced by the current limitations of the semi-3D approach. We constructed the PCA–LDA model using the function lda from the MASS package [100], with PCs 1 through 457 (the maximum number of PCs that could be included without resulting in collinearity) as the continuous explanatory variables and species identity as the dependent categorical variable. The predict function was then used to apply the linear function derived from the LDA for hierarchical clustering and, later, for the examination of morphospace occupation by the five North Atlantic sites. Clustering on LDA output followed the procedure described above for the PCA data.

Assessing topological similarity: morphology, ecology and phylogeny

An analysis of topological similarity was used to assess the overall similarity of morphospace as captured by the semi- and full-3D analyses, and to examine the relationship between planktonic foraminiferal morphology, ecology and phylogeny. These analyses were conducted using the consensus dendrograms for the semi- and full-3D morphological clustering, the ecological clustering, and a pruned version of Aze et al.’s [22] phylogeny of macroperforate foraminifera pruned to the 15 macroperforate species common to all datasets (table 1). Topological similarity was assessed by calculating the path difference (PD) [101] between unrooted topologies using the treedist function from the phangorn (v. 1.99-1) [102] R package. Because the path difference is only well defined for fully bifurcating trees, we randomly resolved multifurcations using the multi2di function from the ape R package. To account for differences arising in tree topology as a result of randomized uncertainty resolution, we conducted 2000 replicates for each dendrogram that required multifurcation resolution and report the average path difference of all replicates.

Assessing phylogenetic signal in ecology and morphology

Phylogenetic niche conservatism describes the pattern of closely related species exhibiting similar ecological traits [103]. To assess the degree of niche conservatism present in our ecological dataset, we calculated Pagel's λ using the time-calibrated Aze et al. [22] phylogeny pruned to include only the 15 macroperforate species that overlap in the Tohoku University and semi-3D datasets (electronic supplementary material, figure S1). These pruned datasets contained 33 full-3D individuals and 414 semi-3D individuals. The Aze et al. [22] time-calibrated phylogeny was built using the data for a fully bifurcating morphospecies tree in the R package paleoPhylo (v. 1.0-108) [104] and pruned using the drop.tip function from the ape R package. Pagel's λ was then calculated for each of the three ecological characters (i.e. symbiont type, habitat depth and geographical range) using the phylosig function from the phytools R package (v. 0.4-31) [105].

Phylogenetic signal was also assessed using Pagel's λ for each individual principal component as a measure of the degree of phylogenetic signal present in the morphological data, and as a proxy for identifying which PCs are most informative in both the full- and the semi-3D datasets. For every PC (39 total for the full-3D Tohoku University dataset and 597 for the semi-3D dataset), Pagel's λ was calculated using the average PC values for each of the 15 focal macroperforate species as described above (results in the electronic supplementary material, table S2). The morphological PCs with high phylogenetic signal and strong support (i.e. λ > 0.5 and p < 0.05) were identified for each dataset (electronic supplementary material, table 2a,d). These high phylogenetic signal PCs are interpreted to capture the phylogenetically informative aspects of morphology. To understand more specifically what aspect of shape variation these phylogenetically informative orthogonal axes were capturing, the five maximum and minimum individuals along the top three high-λ PC axes (PCs 2, 4 and 12 for the Tohoku University dataset, and PCs 2, 17 and 96 for the semi-3D dataset) were identified and considered in turn (electronic supplementary material, table S2b,e).

Table 2.

Comparing topologies of morphological, ecological and phylogenetic clusterings using the path difference pairwise distance metric. Lower distances correspond to more similar topologies.

	path difference
	Tohoku University	coretop/exemplar	ecology	phylogeny
Tohoku University	—
coretop/exemplar	32.882	—
ecology	35.642	24.570	—
phylogeny	32.019	24.429	23.367	—

For each morphological dataset, hierarchical clustering was then conducted for: (i) all of the PCs with high phylogenetic signal and strong support and (ii) the PC with the highest λ value and a substantial amount of morphological variance captured (cutoff of 3% or more) (i.e. PC 4 for the Tohoku University dataset and PC 2 for the semi-3D dataset). The path difference between the high-λ PC dendrogram, the single highest PC dendrogram, and the Aze et al. [22] phylogeny was then calculated for each dataset (electronic supplementary material, table S2c,f). Hierarchical clustering and tree distance calculation was carried out using the same methods described in §2g,i.

Results

Pipeline for semi-3D morphometrics: towards high-throughput community dynamics

With the completion of the beta version of run3dmorph, we have a complete pipeline for the extraction of semi-3D morphometric data (including surface area and volume estimates) from light microscopic and photogrammetric images. The first two components of the pipeline (segment and focus) are written in Python, a free programming language that runs across platforms. run2dmorph and run3dmorph currently execute in MATLAB (version 2015b or above), a proprietary software, but will be ported into Python in future versions. All the software is available, with frequent updates, from GitHub (see https://github.com/Hull-Lab). Key features of run3dmorph, discussed in detail in Material and methods, are aggressive noise reduction routines (figures 2 and 3) and multiple options for estimating surface area and volume given unknown backs to objects (figure 4). From the angle imaged, semi-3D meshes visually reflect the 3D morphology captured by more traditional, CT scanning methods (figure 5).

Semi- and full-3D morphospace in modern planktonic foraminifera

Tohoku University specimens and full-3D morphospace

We constructed the first 3D morphospace for modern planktonic foraminifera using the CT scans of 39 Tohoku University specimens spanning 19 extant taxa (figure 6). Across all of the datasets that we examined, we found that the morphological variation captured by single principal component axes was always quite low: in the case of the full-3D morphospace, the first three PCs cumulatively account for just 21.24% of the total morphological variance, with PC1 accounting for 8.61%, PC2 for 7.32% and PC3 for 5.31%. However, we do not find this low variance captured inherently problematic. In allowing the semi-landmarks to slide, Gonzalez et al. [96] found that the amount of variance capture in the first two principal components fell by about half for rodent molars, rodent brains and primate brains, with sliding 3D-semi-landmark PC1 values ranging from 18–40% of total morphological variance. We likewise found a similar decline in variance capture between non-sliding and sliding semi-landmarks of about half. Even so, the variance captured in these 3D measures of planktonic foraminiferal shape seem low considering that 2D outline methods can capture approximately 75% of morphological variance in the first three principal components (i.e. [19]). Although we have yet to compare 2D and 3D morphometrics in planktonic foraminifera, we suspect that the drop in variance captured is owing to the fact that much of the shape variation in the third dimension is independent of that in the other two dimensions. This remains to be tested in future work.

Figure 6.

Morphospace (PC1 versus PC2) generated from the 39 specimens that comprise the Tohoku University dataset, aligned using the Boyer et al. automatic 3D geometric morphometrics algorithm [91] and sliding semi-landmarks. Coloured lines define convex hulls for each species (for species that include only two individuals, this is manifested as a single connecting line). H. pelagica, Hastigerina pelagica; H. scitula, Hirsutella scitula; G. bulloides, Globigerina bulloides; G. glutinata, Globigerinita glutinata.

The low morphological variance captured by the first few PCs does not affect the robustness of the exploratory morphological analyses, as these are conducted using all non-collinear PCs (i.e. excluding highly correlated PCs). The non-collinear PCs together capture more than 95% of the total variance.

For the Tohoku University specimens, overlap in the relative position of species in morphospace generally occurred between taxa with close evolutionary affinities. For instance, in PC1/PC2 space (figure 6), the two Neogloboquadrina species overlap (N. dutertrei and N. pachyderma) and the two Truncorotalia species overlap (T. crassaformis and T. truncatulinoides) and all the disc-shaped taxa (technically ‘globorotaliform’) cluster in the lower left quadrant of morphospace (e.g. Menardella menardii, H. hirsuta, Hirsutella scitula and Globorotalia tumida). The morphological clustering of some taxa by taxonomic affinity is also clearly apparent in the consensus dendrogram (figure 7a). Compact to spherical forms, spanning a range of taxonomic groups, occupy the centre of PC1/PC2 space and include the closely related Globigerinoides ruber and Globigerinoides conglobatus, along with other taxa (e.g. S. dehiscens, Globoconella inflata, P. obliquiloculata and the Neogloboquardiniids). More lobulate forms (e.g. Globigerina bulloides and Globigerinita glutinata) have the most positive PC1 scores. The largest amount of morphological variation along PC1–PC2 is encompassed by S. dehiscens, a species noted for its unusual crust with supplementary apertures. The next largest amount of variation is encompassed by T. crassaformis, which exhibits a cone-like axial form that varies greatly within the pseudeo-cryptic species complex [53].

Figure 7.

Majority-rule consensus cluster dendrograms for (a) the Tohoku University dataset and (b) our novel coretop exemplar dataset. Tip label colours correspond with species colours in the morphospaces depicted in figures 6 and 8a. The consensus dendrogram is constructed from five cluster dendrograms, each built from the same mean distance matrix using a different clustering algorithm (see text).

Figure 8.

Morphospace generated from the 597 specimens that comprise our novel coretop/exemplar dataset, aligned using the Boyer et al. automatic 3D geometric morphometrics algorithm [91]. (a) The 281 species-identified coretop/exemplar individuals, grouped by genus to aid visualization (note: Trilobatus sacculifer grouped with Globigerinoides); (b) the 116 individuals from four Atlantic coretops; (c) the 200 reproducibility test individuals. Coloured lines define convex hulls for each genus/coretop/species. All three morphospaces are plotted on the same axes (PC1 versus PC2) and scale.

When the phylogenetic signal of each individual PC in the full-3D dataset is assessed using Pagel's λ, the PCs with the strongest phylogenetic signal are PC2, 4, 12 and 39 (electronic supplementary material, table S2a). Although PCs 2 and 4 have a strong phylogenetic signal (λ of 0.66 and 0.82), the species representing the extremes of the PC axes do not sort out according to species or morphotype (electronic supplementary material, table S2). For instance, S. dehiscens appears as a representative species on both ends of the PC2 axis.

Semi-3D morphospace

The 281 individuals from 24 species, the 200 replicate ‘individuals’, and the 116 coretop objects were used to construct the first semi-3D morphospace for modern planktonic foraminifera. As with the full-3D morphospace, the variance captured by the first three PCs was low (14.66% total), with PC1 accounting for 7.6%, PC2 for 4.6% and PC3 for 2.5% of the total variance. A striking difference between the semi- and full-3D morphospaces is the amount of overlap between different species (e.g. figure 8a) along all PCs. The low number of individuals included in the full-3D analyses precludes a quantitative analysis of the relative overlap between species in the full- and semi-3D approaches, but it is clear that while the full-3D analysis resulted in sensible separation of major taxonomic groupings along the principal PCs, the semi-3D analyses did not.

This complete overlap in semi-3D space is also apparent in the complete morphological overlap between North Atlantic sites with very different species compositions (figure 8b), and in the wide variance observed when the same individual is imaged under the multiple reproducibility test conditions, as described in §2b (figure 8c). This latter analysis (figure 8c: repeat measurements on single individuals) reveals the general tendencies of PC1 and PC2 in semi-3D morphospace, with a progression from relatively low (e.g. flat) to high (e.g. domed) umbilical profiles on PC1 and from smooth-edged to lobulated-edged along PC2. Based on Pagel's λ, PC2 is identified to be one of the most taxonomically informative PCs for the semi-3D dataset, along with PC 17, 96, 457, 555 and 588. As with the full-3D data, the taxa loading on the extremes of these PCs are typically mixed (e.g. Globigerinoides ruber appears at both extremes of PC17), with the exception of PC2. For PC2, five individuals of G. ruber have the most positive PC2 scores, and five individuals of Menardella menardii and Globorotalia tumida have the lowest scores.

The high overlap between species along all the PCs directly examined led us to perform a PCA–LDA on the semi-3D data before clustering. The PCA–LDA allowed us to examine the relationships among taxa using the variance relevant to distinguishing among species (cf. figure 7b (PCA cluster) and figure 9b (PCA–LDA cluster)). To consider the relationship between the semi- and full-3D datasets, we performed a second PCA–LDA that included just those 16 species overlapping between the semi- and full-3D datasets (figure 10). This second PCA–LDA emphasizes the difference in morphospace between the full- and semi-3D analyses. There are no species pairs in common between the two analyses (figure 10). Both dendrograms contain a major cluster of nine species, but only four species are in common to both clusters—a result which might be expected by chance.

Figure 9.

(a) Plot of the linear discriminant (LD) space (LD1 versus LD2) for the 24 planktonic foraminifer species present in the semi-3D dataset. Zoomed-in box shows the enclosed species on different axes to show relative cluster positions more clearly. (b) Corresponding majority-rule consensus cluster dendrogram for the LD space shown in (a). (c) LD spaces for coretop individuals, with species identity predicted using the LD model estimated using the species-identified individuals in (a) as the training set. Each coretop is plotted with identical axes to show relative occupation.

Figure 10.

Morphospaces and majority-rule consensus cluster dendrograms for the 16 species that overlap between the Tohoku University dataset (a,c) and our novel coretop/exemplar dataset (b,d) (post-LDA; see text). Coloured lines define convex hulls for each species. Dendrogram tip label colours correspond to morphospace species colours.

Semi-3D morphospace: reproducibility test

To explore the question of ‘other sources of variance’, we consider the results of the reproducibility tests in more detail here. For the reproducibility tests, eight individuals (representing eight different species) were imaged under 25 different conditions each, for a total of 200 ‘individuals’ imaged. The goal of this test was to determine how much of the semi-3D morphological variability could reflect various sources of error introduced during slide preparation and imaging, including variable angles of imaging, image compositing and glare. Ideally, the reproducibility test individuals would cluster closely together in morphospace—while this is somewhat true (figure 8c), together these eight individuals span most of PC1 and PC2 morphospace. Several individuals exhibit particularly high levels of variation: N. dutertrei, Globigerinoides ruber and Orbulina universa (figure 8c). We considered the variability of the perfectly spherical O. universa in greater detail (figure 11a), and it was apparent that the large variation exhibited by O. universa is due in large part to two outlier individuals (marked in red). Examination of the meshes of these two outlier individuals reveals obvious errors in the integrity of the extracted mesh (figure 11b), which appear to have resulted from imaging artefacts during z-stack focusing (figure 11c). These artefacts appear as a smeared, unfocused area on the object surface in the focused image. The error is then perpetuated through the mesh extraction pipeline via a poorly constructed height map (figure 11d). The pathological nature of the two outliers is clearly distinguishable via comparison with non-pathological individuals (figure 11e–g).

Figure 11.

Exploration of the reproducibility tests of the Orbulina universa individual (a) The PC1 versus PC2 coordinates for O. universa from figure 8c. The points highlighted in red are outliers for which the run3dmorph-extracted semi-3D half-hulls exhibited pathologies (see text). (b–d) Extracted mesh, focused image and height map showcasing pathologies arising during z-stack focusing. In (c), the smeared, unfocused portions of the object are outlined in white. (e–g) The corresponding mesh, focused image and height map for a properly extracted O. orbulina individual.

Without the two poorly extracted individuals, the variation encompassed by the O. universa is comparable to that of the relatively low variation individuals of T. truncatulinoides, Trilobatus sacculifer and Globigerinella siphonifera (figure 8c). All of these low variation individuals have relatively flat spiral surfaces and, in two cases, small (T. truncatulinoides) or edge-facing (G. siphonifera) apertures—both factors that increase the consistency of the image perspective. Globigerinoides ruber, N. dutertrei and Globorotalia tumida have some combination of large apertures and domed-spiral sides, both of which dramatically affect the shape of the semi-3D mesh extracted from slightly different viewpoints.

Full- versus semi-3D morphospace

To facilitate the direct comparison of the full- and semi-3D morphospaces, we trimmed the two datasets down to include just the 16 species in common to both datasets. These 16 species versions were clustered, as before, on the PCA data for the full-3D dataset and on the PCA–LDA data for the semi-3D dataset (figure 10). The most readily apparent feature of the morphological consensus for both datasets is how little the clustering appears to relate to taxonomic affinities. There are, of course, exceptions, with N. dutertrei and Neogloboquadrina pachyderma paired in the full-3D analysis, and Menardella menardii and Globorotalia tumida paired in the semi-3D analysis. Beyond this, both trees have structure that can be interpreted as reasonable given the morphology of the taxa in related clusters, but overall the groupings were unexpected given our morphological understanding of these species.

Relationship between gross morphology, ecology and phylogeny in 15 extant species

The PD tree distance metric was used to quantitatively assess the similarity between the two morphological dendrograms (semi- and full-3D) and to examine the relationship between 3D morphology and the phylogenetic and ecological relatedness of taxa (table 2). As a Euclidean distance-based metric of tree similarity, each tree is represented as a vector of pairwise edge distances between all terminal taxa pairs and is thus compared. This comparison is purely topological (i.e. edge lengths are not accounted for) and provides a direct metric for assessing the topological similarities between the morphological, ecological and phylogenic dendrograms.

The two most similar dendrograms were the ecological and phylogenetic trees (table 2). This result suggested a strong signal of phylogenetic niche conservatism, a possibility we tested further using Pagel's λ on the three ecological characters (symbiont type, habitat depth and geological range) individually. In the subset of 15 macroperforate species of planktonic foraminifera examined morphologically, we found phylogenetic signal in all ecological traits: symbiont type (λ = 1.049; p = 0.002), habitat depth (λ = 1.106; p = 0.002) and geographic range (λ = 0.565; p = 0.046).

Both morphological dendrograms (i.e. full- and semi-3D) are more similar to the phylogeny (i.e. exhibits the shortest PD between the two trees) than they are to one another or to the ecological dendrogram. These results suggest a stronger influence of evolutionary history on morphology than of ecology (figure 12c). It also implies that the morphological structures captured by the full and semi-3D approaches differ, given the relative dissimilarity between these two morphological dendrograms. To our surprise, the path difference suggests a greater similarity between the semi-3D dendrogram and phylogeny (PDsemi-3D = 23.791) than between the full-3D dendrogram and phylogeny (PDfull-3D = 31.081). This same relative ordering is also true of the ecological similarity: there is greater similarity between the ecological dendrogram and the semi-3D morphological dendrogram (PDsemi-3D = 24.207) than with the full-3D morphological dendrogram (PDfull-3D = 34.728).

Figure 12.

Ecological and phylogenetic clustering of the 15 macroperforate species that overlap between the Tohoku University dataset and the coretop/exemplar dataset. (a) Consensus cluster dendrogram of the full-3D Tohoku University specimen morphospace (same as figure 10c). (b) Ecological cluster dendrogram built using Jaccard distances calculated from three ecological traits (table 1). (c) The phylogenetic relationships between the 15 macroperforate species, as pruned and redrawn from Aze et al.'s [22] stratophenetic phylogeny. Dendrogram tip label colours correspond to morphospace species colours from figure 9.

The higher congruence between the semi-3D morphological dendrogram and the phylogeny, as compared to the full-3D data, also exists in the morphological dendrograms based on high phylogenetic signal PCs only (λ > 0.5 and p < 0.05). When the high phylogenetic signal dendrograms are considered (PCs 2, 4, 12 and 39 for the full-3D dataset; PCs 2, 17, 96, 457, 555 and 588 for the semi-3D dataset), the path difference between the semi-3D dendrogram and the phylogeny is lower than that between the full-3D dendrogram and the phylogeny (PDsemi-3D = 24.576 versus PDfull-3D = 26.609) (electronic supplementary material, table S2c,f).

Assemblage-wide ecometrics: community structure in four coretop locations

In addition to generating 2D-outlines and semi-3D meshes for downstream morphometric analyses, the run2dmorph and run3dmorph software can also automatically generate estimates of 2D size (e.g. major and minor axes length, enclosed area, perimeter length) and 3D size (volume and surface area). As a key macroecological trait, body size is often measured in fossils either by species exemplars [106,107] or by 2D metrics like major axis length or 2D surface area [108]. For the 9681 individuals included in the community structure analysis from the four coretop locations, the potential uncertainty owing to the unknown back-morphology is estimated as %volumetric difference between the low- and high-end estimates. We find that the high-end estimate (cylindrical) is on average 23.55% (σ = 11.89%) larger volumetrically and 3.19% (σ = 1.40%) smaller in surface area than the low-end estimate (conical), normalized to the dome estimate. For species like planktonic foraminifera, which range from nearly completely flat to completely spherical, we consider the importance of estimating volume (over 2D area) with a direct comparison of 2D and 3D metrics of community size distributions (figure 13); these comparisons are discussed in the following sections.

Figure 13.

Kernel density plots for (a) 2D outline area and (b) semi-3D volume (assuming spheroidal dome base) for the four Atlantic coretops: CH 82-21 (n = 2879), KC 78 (n = 1768), EW 93-03-04 (n = 3034) and AII 42-2-2 (n = 2000).

Discussion

In recent years, advances in imaging and image analysis have driven forward our understanding of diversity dynamics in the history of life by unlocking previously inaccessible aspects of the fossil record [109]. There have been three main foci of efforts in this field of virtual palaeontology: (i) the automatic recognition of taxa [110-112], with arguably the most successful application of this approach in calcareous nanofossils [113,114]; (ii) the collection of refined morphometric data [115], with technologically advanced approaches like X-ray and synchrotron CT at the forefront [116,117] and (iii) the rapid assessment of assemblage body size distributions [118], typically in 2D silhouettes [119,120]. High-throughput morphometric approaches have generally lagged because in morphometrics, as in most things, the devil is in the details. Small differences in specimen orientation can introduce more variation into a dataset than evolution itself, and many taxa are differentiated at the species level by fine features missed by conventional methods like light microscopy.

In full recognition of the ultimate need for precision, we set out to develop a high-throughput approach (semi-3D morphometrics). We reasoned that image processing and machine learning will ultimately allow us to extract the ecological and evolutionary signals from the noise and that rapid approaches are critically needed to address questions of assemblage dynamics (see call for population studies in [109,121]). Towards this end, we have developed and tested a method for the extraction of 3D hulls from light microscopic and photogrammetric images, examined the potential of morphology as a tracer of assemblage dynamics in planktonic foraminifera, compared rapid and refined 3D morphometric approaches (e.g. semi- versus full-3D meshes), and explored a preliminary application of these methods in the North Atlantic. We discuss each in turn.

Towards high-throughput community dynamics

Coupled with upstream functions in the AutoMorph software suite, the code presented here (run3dmorph) can process a slide scan containing thousands of individual fossils and extract individual semi-3D meshes and volume and surface area estimates for all unique objects (figures 2–5). This code is flexible enough to work with any light-coloured object imaged on a black background, and we are currently using it in-house with photographs of limpet assemblages. All code is freely available on GitHub (https://github.com/HullLab) and is being actively updated to fix known bugs, increase functionality and speed, and remove dependencies on proprietary software. Our reproducibility tests highlight the robustness of the semi-3D mesh extraction code—the two major problems encountered (artefacts during z-stack image focusing and/or tilt variation during sample preparation prior to microscope imaging) do not result from the 3D mesh extraction code itself. Rather, the 3D mesh extraction behaves as expected, with the quality of output determined by the quality of the input (an example of mesh extraction resulting from optimal versus suboptimal inputs is shown in figure 11). However, regardless of where the errors occur in the high-throughput approach, they need robust solutions to enable high-throughput study of community dynamics. We discuss ongoing work to address these problems below. At present, in light of the reproducibility test results, users of the run3dmorph software are advised to exercise care during initial set-up and positioning of specimens prior to imaging in order to reduce the amount of variability due to tilt, and to inspect 3D PDFs of height map meshes to ensure proper mesh extraction.

Already, run3dmorph provides a robust means of measuring assemblage body volume distributions. Even without measuring the backs of objects, we are able to measure the volume of fossils with an uncertainty of approximately 20% introduced by estimating back-half shape. In other words, an individual with a fully filled (i.e. cylindrical) back is approximately 40% more volumetric than the cone-backed equivalent. Although still somewhat uncertain, this is a dramatic improvement on the existing 2D silhouette volume-estimation approaches for assessing assemblage body size distributions, given the great variation among foraminiferal species in the depth dimension. In modern tropical assemblages, some of the largest individuals, as measured by 2D length and width, are also the flattest (i.e. Menardella menardii, Globorotalia tumida, H. scitula, to even relatively flattened taxa like Trilobatus sacculifer and Globigerinella siphonifera). By assuming that 2D silhouette captures the true body volume, we calculate that 2D methods effectively overestimate the body mass in flat individuals by roughly 375% (as calculated for Menardella menardii assuming spherical equivalent). As such, our approach reduces the volumetric uncertainty by more than ninefold, yielding body size measurements that are more accurate and better suited to many of the ecological and morphological questions targeted with body size data. Through its effect on metabolic rates and intraspecific interactions, body mass is considered to be an ecological trait of first-order importance [122,123]. We therefore see the volume (and surface area) estimates of run3dmorph as one of the most useful, and certainly most readily applicable, outputs of the code.

Two major advances are planned for the AutoMorph software in order to enable accurate 3D reconstructions of assemblages morphometrically by targeting the two major problems in the current analyses: noise and perspective. First, an additional noise reduction step, based on inferred surface rugosity, will be used as a final check for both pitting (as in the Orbulina example in figure 11) and glare peaks. Second, we aim to use the semi-3D meshes as hulls for matching, and warping, full-3D shapes to fit. In this way, the semi-3D information can be used to fit and interpret the full 3D-shape of the geometrically formed planktonic foraminifera.

Determinants of gross morphology in planktonic foraminifera

Planktonic foraminiferal taxa differ in gross and fine morphological features—only the first of which is accounted for here. Fine morphological features beyond our ability to measure include differences in foraminiferal wall construction, including wall porosity, thickness, texture and the presence or the absence of spines. Many of these fine morphological features have hypothesized (and in some cases, well-tested) functions, such as the use of pores for metabolic gas exchange, wall ornamentation for cytoplasmic anchoring and spines (or other surface ornamentation) for exposing symbionts to sunlight (in photosymbiotic taxa) and prey capture [60]. In future work, we aim to consider the importance of this omission through the inclusion of categorical species-specific fine morphological associations, but at present we have focused on what we can measure directly: gross morphology.

Gross morphology can vary dramatically among planktonic foraminifera, with species ranging from completely spherical (see Orbulina universa in figure 11), to pyramidal (see T. truncatulinoides in figure 5), to globular (see N. dutertrei in figure 5), discoidal (i.e. Menardella menardii), or digitate (i.e. finger-like projections in taxa like Hastigerinella digitata). In contrast to fine morphological features, surprisingly little is known about the functional importance of these gross morphological features. Earlier studies had often hypothesized hydrographic functions for these morphologies (as summarized in [34]), but these were recently elegantly laid to rest [33]. Given the importance of interspecific interactions in plankton communities [59], we suspect that interspecific associations (mutualistic or commensal), predation and oxygen limitation are likely the dominant forces shaping the evolution (and re-evolution) of gross morphology in planktonic foraminiferal lineages.

In this context, the analysis of the gross morphological similarities among 19 taxa (out of the roughly 40 morphological species globally) of modern planktonic foraminifera is particularly interesting for highlighting morphological affinities among taxa (figures 6 and 7a) and for exploring morphological, ecological and phylogenetic relationships (figure 12). Along the first two morphological axes (PC1 and PC2), the globorotaliform species cluster distinctly in the lower left quadrant, from the very flat Menardella menardii (most negative values along both axes) through the relatively inflated hirsutellids and Globorotalia tumida, to the umbilically domed truncorotaliids (figure 12). The relatively spherical taxa in this analysis (Globigerinoides conglobatus, N. dutertrei, C. nitida, Globoconella inflata and P. obliquiloculata) all cluster in the middle of PC1/PC2 space, with the relatively lobulate taxa along the edges with more positive PC1 (Globigerina bulloides, Globigerinita glutinata) or PC2 (Globigerinoides sacculifur, Globigerinoides ruber, Globigerinella siphonifera). Aside for the clustering among the taxonomically related globorotaliform taxa (along PC1 and PC2; figure 6), the phylogenetically mixed nature of morphospace occupation along the first two PCs is reinforced by the consensus clustering (figure 7a). Consensus clustering used all non-collinear PCs, which together account for more than 95% of the morphological variance. In the morphological consensus dendrogram, aside from two close taxonomic pairs (i.e. Hirsutella and Neogloboquadrina), the most striking feature of the dendrogram is the lack of apparent taxonomic or ecological cluster of the data, or obvious morphological affinities of the dominant clades—an observation supported by the distant relationship of full-3D morphology to both ecology or phylogeny alone (table 2). There appears to be more similarity between the ecological consensus dendrogram and the phylogeny than the morphological dendrogram and the phylogeny.

Rather than decoupling morphology from these two important influences (i.e. phylogeny and ecology), we suspect that this result emphasizes four aspects of planktonic foraminiferal evolution and morphology. First, major ecological features like depth habitat, biogeography and the presence/absence of symbionts appear to be relatively conservative. The relative similarity of the ecological and phylogenetic trees, and the tests of Pagel's λ on the three ecological characters (table 1) included in this study suggest that there is very high phylogenetic signal in these ecological traits, particularly for depth habitat and symbiont type. Second, in spite of widespread convergence in gross morphology [30,124], it is still unknown whether this convergence is functional or due to morphological constraints. It is also possible the most relevant ecological factors controlling morphology (like inter- and intraspecific interactions and metabolic demands) are the most poorly understood and therefore not considered. Recent studies (e.g. [125,126]) emphasize the importance of considering intraspecific variation in characterizing community ecology in the light of phenomena such as trait distribution overlap and niche packing. Third, our analysis is just a first attempt to explore this problem in planktonic foraminifera. For instance, we considered the relationships between morphology, ecology and phylogeny in just the 19 taxa with freely available 3D data via the Tohoku University Museum. Including all extant taxa (roughly 40–50 morphospecies) would provide a comprehensive test of the problem, as would including other time periods to explore the issue of morphological and ecological convergence. Fourth and finally, our method is designed to extract gross morphology and, as previously discussed, the traits used to segregate foraminifer species also include fine morphological features such as wall texture and spine presence/absence. These fine features are either impossible to extract by our method, or, if extractable, are not well resolved and therefore account for very little of the total morphological variation. By contrast, the gross morphological features that dominate the 3D shapes—including the most obvious first-order ones such as relative inflatedness/flatness—are often not highly informative taxonomically.

The tendency of gross morphology to capture many aspects of shape variation unrelated to ecology and phylogeny is unsurprising and would be a default expectation in many taxonomic groups. There are several reasons for this. First, even when functional morphology is well understood, there is no reason to expect a direct one-to-one mapping of morphology and ecology. In many cases, it is the interaction of multiple morphological features that underpins functionality, so a wide range of morphologies could have the same function. Second, gross morphology may poorly capture the underlying features that directly relate to ecological function or to phylogenetic relationships. For instance, ontogeny can greatly affect the first-order gross morphology, but the morphological differences between a juvenile and an adult (e.g. relative proportions, roundedness of features, number of segments and limbs, etc.) are not the features that are used to distinguish among species. For example, we suspect that a gross exterior morphological analysis of canids (e.g. dogs, wolves, foxes, coyotes, jackals, etc.) in varying ontological stages would likely find a large amount of variance related to the difference between adult and juvenile forms owing to the conspicuous differences between adult and juvenile builds and proportions (e.g. snout length, limb proportion, etc.). The importance of ontogeny to morphology is also well known in foraminifera, with juveniles showing a high degree of morphological conservatism when compared with adults [127-129]. Widespread morphological variation is well known intraspecifically as well in planktonic foraminifera—with, for instance, food availability dramatically influencing adult morphology in T. sacculifer [130]. In other words, the default expectation in a gross morphological analysis such as ours is that multiple factors are likely to influence shape, only some of which (and perhaps the least of which) are taxonomically or ecologically informative.

A practical implication of this phenomenon for principal component analyses of 3D morphology would be that the PCs capturing most of the variation may not be taxonomically informative. A first glance at the morphological phylogenies (figures 9, 10 and 12) would certainly seem to support this. To explore this idea further, we assessed the phylogenetic signal present in each individual PC in the semi- and full-3D datasets using Pagel's λ. We found four taxonomically informative PCs in the full 3D analysis (PCs 2, 4, 12 and 39) and six in the semi-3D analysis (PCs 2, 17, 96, 457, 555 and 588). In neither case was PC1 found to exhibit high phylogenetic signal, despite PC1 encompassing the highest amount of variation. However, PC2 was informative for both datasets, with a clear separation of the 3-chambered, globular Globigerinoides ruber and the bi-covex, multichambered Globorotalia tumida and Menardella menardii along the semi-3D PC2 (other morphological distinctions along phylogenetically informative PCs were difficult to interpret; see the electronic supplementary material, table S2). We repeated our clustering and tree distance analyses using only the most informative PCs for each dataset. For the full 3D-dataset, this resulted in greater congruence between the morphological and phylogenetic clusterings (PD = 32.019 versus 26.609). For the semi-3D dataset, there was little difference (PD = 24.429 versus 24.576).

Efficacy of semi-3D morphometrics: semi- and full-3D morphospaces

Semi- and full-3D approaches provide a very different perspective on gross morphological affinities between modern planktonic foraminifera (figure 10), with only the second being based upon the full-3D shape of individuals. For some questions, particularly those of size variation and morphospace occupation, the large differences in morphological space quantified by semi- and full-3D approaches may not matter. However, we generally think that methodological advances, similar to those described in §4a (towards high-throughput community dynamics), are needed to extract full-3D-like data from the semi-3D hulls. There are two reasons for this: first, some taxa have rather similar 3D shapes, like the pyramidal form of both T. truncatulinoides and H. hirsuta, which are reversed with regards to the umbilical/spiral axis. By ignoring half of an object's shape in semi-3D morphometrics, this inherent similarity is missed (compare figure 10c and d with regards to the inferred similarity of T. truncatulinoides and H. hirsuta). Second, the critical problem of viewing angle on inferred morphology in the semi-3D analyses is most easily addressed by projecting semi-3D morphologies into full 3D-shape space, rather than attempting to orient each half-hulled individual to the same perspective.

Given those caveats with regard to the current semi-3D morphological space, it is worth noting that the semi-3D data (after LDA) exhibits more structural similarity to the phylogenetic and ecological clustering than the full-3D shape space (table 2). This result is unexpected given the high resolution and completeness of the Tohoku University specimens compared to the semi-3D meshes extracted using our novel method. One possible reason for this result is the effective difference in sampling density between the full- and the semi-3D datasets. There are two aspects of sampling density that differ between the full- and semi-3D analyses. First, we examine hundreds as opposed to tens of individuals in the semi- versus full-3D analyses. By sampling more individuals in the semi-3D analysis, we likely captured more intra- and interspecific variation. Our results would thus suggest that it may be more important to optimize the number of individuals sampled over the quality of the 3D data, when such a trade-off must be made. Second, by using the same number of semi-landmarks (i.e. 256) for both the full- and the semi-3D meshes, we may have effectively under-sampled morphology in the full-3D dataset, as the 256 semi-landmarks are spread across the entire foraminifer rather than across only the top half of the foraminifer. In this sense, although the raw mesh resolution of the full-3D dataset is higher than that for the semi-3D dataset, the effective mesh resolution used for morphospace construction may actually be lower for the former than the latter. Further work is required in order to parse out these possible effects.

Regardless, it seems likely that for some questions, semi-3D shape spaces, like those briefly explored in figure 9c, might suffice for tracking community dynamics, although much future work is needed to explore this potential. By contrast, assemblage volume, discussed next, is a readily applicable trait, whose measurement is already enabled by the run3dmorph software.

Assemblage-wide ecometrics: quantification of community volume

Body size is a key determinant of interspecific interactions and alone can provide an important metric of biodiversity dynamics (and biotic drivers) through time [131]. For instance, assemblage dwarfing is one of the most coherent cross-clade responses to mass extinction and other abrupt perturbations (e.g. [132,133]). Similarly, over the course of the Cenozoic, there is a massive change in body size of well-fossilized marine plankton [131], with a dramatic increase in planktonic foraminiferal body size since the Miocene [108] and a Cenozoic-long decrease in lith size of coccolithophores [131,134]. Within extant assemblages, foraminiferal body size is correlated with latitude, with high latitudes corresponding to smaller average assemblage sizes and low latitudes corresponding to larger average assemblage sizes [135]. Even in the largest, most diverse modern assemblages in subtropical to tropical regions, adult body size separates otherwise co-occurring taxa [136], such that increasing assemblage diversity coincides with increasing size disparity.

Although many patterns and relationships have been observed, one limitation of most population-scale body size studies to-date is that the 2D silhouette is used as an implicit proxy for total size (i.e. volume or organic mass). In the case of foraminifera, and most shelly fauna, this is only approximately true. In planktonic foraminifera, for instance, relatively flattened morphologies are most common in the taxa with the largest 2D silhouettes. These taxa include most of the globorotaliform species and a number of other relatively un-inflated forms like Globigerinella siphonifera and Trilobatus sacculifer, all with higher surface area to volume ratios than are observed in relatively inflated taxa. Because of variation in the depth dimension, 3D measures of body volume are needed to check patterns inferred from 2D measures, a problem also discussed by Caromel et al. [127].

A first exploration of the relationship between 2D silhouette area and semi-3D volume in four North Atlantic assemblages (figure 13) reveals important differences between size distributions as inferred from 2D and 3D data. In both cases, the northern-most sites (CH 82-21 and EW 93-03-04) have the smallest maximum body sizes and a lower relative proportion of individuals in intermediate-large body sizes, with the northern-most site (EW 93-03-04) exhibiting smaller maximum body sizes than the next most northern site (CH 82-21). In the upper end of the size distributions, it is notable that the use of 3D body size estimates reduces the overall effect of this latitudinal trend on body mass. In both the 2D and the 3D data, the site with the largest body size (KC 78) is approximately fivefold larger than the site with the smallest (EW 93-03-04). If body volume increased equally in the depth dimension, however, the volumetric comparison (figure 13b) should have an 11-fold increase and not the roughly fivefold difference observed. In other words, tropical to subtropical assemblages are large, but they are very flat by comparison to high-latitude assemblages. This ‘flattening’ effect perhaps reflects the effects of metabolic constraints on morphology. While more data are certainly needed to explore the relationship between 2D silhouette size and body volume, this first analysis suggests that semi-3D volume data will be important for understanding the determinants of body size evolution through time, and for quantifying constraining inferences with regards to body mass distributions of populations in time and space. Excitingly, this semi-3D method can be seen as additive to the existing (much faster) 2D approaches, expanding rather than replacing the types of questions that can be addressed.

Conclusion

The main impetus for this study was to develop a high-throughput method for quantifying biotic dynamics in the abundant shelly-fossil record of life. To this end, we developed and presented the run3dmorph code for extracting 3D meshes and estimating surface area and volumes from light and microscopic images, tested the utility of planktonic foraminiferal shape as a proxy of community dynamics, and explored the importance of 3D measures of body size. Although this case study focused exclusively on planktonic foraminifera, run3dmorph is applicable for any object imaged with a series of z-stacks of known height. This broad applicability of run3dmorph, in combination with automated landmark-generation methods such as the Boyer et al. [91] auto3dgm/PuenteAlignment algorithm, is a step forward towards laying the groundwork for a ‘Next-Generation’ approach to morphometrics. Our semi-3D approach sacrifices data resolution and completeness relative to full-3D approaches like CT scanning, but has the advantages of accelerated processing speed, higher data density and reduced expense. Together, these improvements will allow a wider range of researchers to conduct robust analyses on population-level problems addressing the dynamics of 3D community size in time and space. At this point, much future work is still needed to make this semi-3D approach as robust and reliable as full-3D approaches to morphometrics, but the path to this ultimate goal is already clear.

For planktonic foraminifera, our early analyses of morphological patterns and body size variation are intriguing, as they suggest remarkably little similarity between overall gross morphology (in the Tohoku University specimens), phylogeny and certain aspects of ecology. This finding, combined with the pronounced ‘flattening’ of taxa towards tropical regions, emphasizes an exciting set of unknowns with regard to planktonic foraminiferal functional morphology, whose importance is underscored by the use of planktonic foraminifera as a palaeontological model species of macroevolution.