ADVANCING INDIVIDUAL FORAMINIFERA ANALYSIS BY COMBINING MOLECULAR, MORPHOMETRIC, AND TRACE ELEMENT GEOCHEMISTRY

Fossil planktic foraminifera are commonly used as proxies in paleoclimate research because their calcareous tests record physicochemical conditions in the ocean at the time of calcification (Katz et al., 2010). Modern extant species are used to validate how test isotopic and trace element geochemistry correlate with different oceanographic parameters, such as temperature, salinity, primary productivity, and carbonate chemistry; these relationships are routinely applied to the fossil record (Spero et al., 1997; Lea et al., 1999; Russell et al., 2004; Allen et al., 2016; Osborne et al., 2016; Davis et al., 2023). However, molecular phylogenies built from the small subunit (SSU) ribosomal RNA (rRNA) gene have also shown that modern foraminifera species contain cryptic genetic diversity, which is difficult to detect morphologically (Darling & Wade, 2008; Ujiié & Lipps, 2009). Indeed, molecular studies (Huber et al., 1997; de Vargas et al., 2002; Darling et al., 2006; Darling & Wade, 2008; Morard et al., 2009) have been shown to change a number of species classifications. Studies suggest these genotypes reflect distinct ecotypes with different habitats, calcification depths, and diets, which collectively determine test geochemistry (de Vargas et al., 2002; Morard et al., 2013; Sadekov et al., 2016). However, it remains to be determined how genetic diversity influences test geochemistry and the resulting paleoceanographic interpretations.

Recent advancements in analytical tools allow the analysis of individual foraminifera tests, which has the potential to reveal high-resolution population variability, but interpretation of these results is potentially complicated by intraspecific variability (Katz et al., 2010). Combining molecular and geochemical approaches could yield unique insights for elucidating mechanisms responsible for intraspecific variability. Ujiié et al. (2019) tested DNA extraction methods to assess their impact on test morphology and stable isotope geochemistry compared to non-extracted samples. They found that DNA extraction does not alter test density or stable isotope ratios (Ujiié et al., 2019). Additionally, they noticed a significant offset in carbon isotopic composition among Globigerinoides ruber (d’Orbigny, 1839) genotypes, suggesting the presence of different ecological traits or ‘vital effects’ (Ujiié et al., 2019). Although isotopes are commonly used for temperature and other environmental reconstructions, trace element proxies can yield further oceanographic insights and validate temperature interpretations (Katz et al., 2010). In this research, we advance the Ujiié et al. (2019) workflow for analyzing individual foraminifera to test how trace element geochemistry is impacted by DNA extraction. We focus our work on Globigerina bulloides (d’Orbigny, 1826), a planktic foraminifera species with cryptic diversity that is widely used for paleoceanographic reconstructions.

Globigerina bulloides is a shallow-dwelling, spinose foraminifera species found in temperate to subpolar waters with seven identified genotypes (Darling et al., 1999; Darling et al., 2000; Darling & Wade, 2008; Schiebel & Hemleben, 2017). The seven genotypes are grouped into two broader types with temperature affinities—the warmer water Type I and cooler water Type II (Darling & Wade, 2008). Previous research has already compared the geochemical signatures of G. bulloides test with the temperature affinities of the two known genotypes and found a bimodal distribution of predicted temperature based on δ¹⁸O and Mg/Ca trace element geochemistry (Sadekov et al., 2016). The researchers could not directly compare genetic and geochemical signatures of individual tests because the DNA extraction method they used dissolved the test. However, these results suggest confirming genotypic geochemical signatures would improve paleoclimate interpretations.

Specimens in this study were collected from the Southern California Bight, an upwelling region in the NE Pacific where multiple, co-occurring colder water Type II genotypes of G. bulloides have been identified in previous research (Darling et al., 2003). Type IId is prevalent in the region and has been found through upwelling and downwelling periods (Darling et al., 2003). Type IIa has been identified once in the Southern California Bight during the winter season (January) but typically prevails in the polar oceans (Darling et al., 2007). Another genotype found in the Northeast Pacific is Type IIe, which has been found to extend farther north in the North Pacific subpolar gyre (Darling & Wade, 2008).

Northeast Pacific G. bulloides genotypes are morphologically indistinguishable from each other, meaning they cannot be distinguished by microscopy. However, some studies have suggested the colder-water type IIa might be more heavily encrusted than type IId, which may account for differences in geochemistry and, notably, area-normalized test weight, a proxy for carbonate ion concentration (Sautter & Thunell, 1991; Osborne et al., 2016). Osborne et al. (2016) used size-normalized test weight to identify the thicker encrusting genotype, since the differences were not visually discernable. No molecular data have been paired yet with morphological or geochemical data to corroborate that type IIa is more heavily encrusted.

Here, we advance and test a workflow for analyzing individual foraminifera for molecular, morphometric, and trace element geochemical analysis on G. bulloides. We specifically make use of X-ray micro-computed tomography (XMCT), a non-destructive technique which offers new opportunities to discern morphological variability. XMCT scanning provides high-resolution information about test thickness, density, and morphology, including the ability to project 2D areas and 3D images. The mean computed tomography (CT) number is the average intensity value for each test determined by the X-ray attenuation coefficient and reflects minor changes in test density (Iwasaki et al., 2015). The mean CT number has been proposed as a proxy for test calcification intensity and with ambient seawater carbonate chemistry in G. bulloides specimens (Iwasaki et al., 2015, 2019a, 2019b).

In this study, we tested whether DNA extraction influences trace element geochemistry by dividing 60 adult G. bulloides from a single net tow in the California Bight (Fig. 1) into an extraction and control group. Trace element ratios were analyzed by Laser Ablation Inductively Coupled Plasma Mass Spectrometry (LA-ICP-MS), which is minimally destructive and leaves potential for additional analyses on the same test. Following Ujiié et al. (2019), we also collected information about test morphology and density using XMCT. This workflow allows for highly-resolved genetic, trace element and morphometric data from a single foraminifera test, and thus provides tools to distinguish trace element variability and morphometric features in a genetically diverse morphospecies like G. bulloides.

Foraminifera samples were collected from a single net tow between 0–55 meters deep using a one-meter, 200-µm mesh net off the University of South California’s Wrigley Marine Science Center at Catalina Island, CA at 33°28.9′N, 118°29.9′W (Fig. 1). The sea surface temperature, measured at the time of collection, was 21.0°C (Table 1). Sixty live individual G. bulloides with colorful cytoplasm were picked from the net contents, rinsed with 0.2-µm filtered seawater and photographed (Fig. 2). Eight additional G. bulloides specimens, described below, were collected for genetic analysis on different dates and in different regions of the Northeast Pacific with the aim of capturing genotypic variability (Table 1).

The control group specimens (∼30 specimens) were rinsed in sodium hydroxide buffered deionized water to remove organic material and air dried. Three specimens were damaged or lost in subsequent analysis. The other half of the specimens was placed in individual vials containing GITC* buffer for DNA extraction and stored frozen until further processing. DNA extraction followed the standard protocol (Morard, 2010; Weiner et al., 2016). After extraction, 27 of the 30 foraminifer tests were intact; these specimens were rinsed in deionized water and air dried. Two specimens did not sequence successfully. The following cleaning, imaging, and geochemical analysis treatments were the same for control and extraction tests.

A small number of additional G. bulloides specimens was collected on opportunistic cruises in the North Pacific and included in the DNA extraction group to see if we could identify other genotypes in the region (Fig. 1; Table 1). Five specimens were collected from the same study region on a different date (CAT 02-04, 06, 08; Table 1). Two specimens (NCC Y01 & Y02) were collected at ∼42°N and one specimen (NCC SPR18) was collected at 44°N during an upwelling event (Table 1). Due to the small sample size, these specimens were excluded from statistical and geochemical analyses but discussed qualitatively in terms of morphological differences.

Polymerase chain reaction (PCR) was carried out using the Phire™ or Phusion™ HotStart polymerase enzyme (Thermo Fisher Scientific) with spinose-specific primer S14p (5′ – AAGGGCACCACAAGMGCG – 3′) and universal primer SBf (5′ – TGATCCATCRGCAGGTTCACCTAG – 3′) from Weiner et al. (2016).

The PCR cycling parameters for the Phire polymerase enzyme, with slight variations for Phusion denoted in parentheses, included an initial denaturation at 98°C for 30 seconds, followed by 35 cycles of denaturation at 98°C for 5 seconds (10 seconds), annealing at 65°C for 5 seconds (30 seconds), and extension at 72°C for 20 seconds (30 seconds). The final extension at 72°C lasted for 1 minutes (5 minutes). The PCR products were visualized on a 2% agarose gel. At least two successful PCRs from each sample were pooled for sequencing. PCR cleanup used the Qiagen QIAquick® PCR Purification kit. DNA sequencing was completed using ABI Prism 3730 Genetic Analyzer with BigDye Terminator v. 3.1. Cycle Sequencing Kit at the Center for Quantitative Life Sciences at Oregon State University.

Nucleotide sequences were visually inspected and trimmed using Geneious (Geneious Prime 2022.0.2, 2022). To identify genotype, sequences were queried using the National Center for Biotechnology Information Nucleotide Basic Local Alignment Search Tool or BLAST (Altschul et al., 1990). Sequences were submitted to GenBank (PP249933-PP249968).

To explore the genetic similarities of the extracted samples, genetic analysis included nucleotide sequences from all extracted samples (G##), other Catalina Island samples (CAT##), other G. bulloides Northeast Pacific samples (NCC), and reference G. bulloides sequences pulled from the Planktonic Foraminifera Ribosomal Reference (PFR²) database (Morard et al., 2015). Sequences were aligned using the MUSCLE multiple alignment algorithm with 10 iterations (Edgar, 2022) across a ∼900 bp region. A neighbor-joining phylogenetic tree was built using the Geneious Tree Builder algorithm using a Tamura-Nei model (Tamura & Nei, 1993) and a reference Type Ia G. bulloides sequence as the outgroup. To estimate the significance of the branches, the tree was resampled 1000 times using bootstrapping with a random seed. Only bootstrap values with a support threshold >50% were included.

Both control and extraction specimens were oxidatively cleaned to remove any residual organic matter, following previously established protocols (Bonnin et al., 2019). Samples were immersed for 10 minutes in 1:1 mixture of NaOH and H₂O₂ in a water bath at 60°C, sonicated briefly, then triple-rinsed in deionized water and air dried. Samples were weighed three times using a copper weigh boat on a Sartorius Ultra Microbalance SE2 with a tolerance of 0.7 μg (Dale Weighing Systems, Wood Dale, IL). Sample weights represent the mean of the three repeated measurements.

To acquire high-resolution XMCT scans, all samples and a NBS 19 calcite standard were mounted on a peg with carbon tape. Specimens were scanned on a Hamamatsu L10711-19 (Hamamatsu Phototonics Corporation, Shizuoka, Japan) at the School of Chemical, Biological and Environmental Engineering at Oregon State University. Each scan took approximately 4.5 hours and had a spatial resolution of 1.9637 μm per pixel. Additional scanning electron microscopy (SEM) imaging was conducted on a subset of control and extraction group samples using a benchtop Hitachi SEM TM-4000 (Hitachi, Ltd., Tokyo, Japan) at the Advanced Technology and Manufacturing Institute at Oregon State University.

XMCT image processing was completed using the imaging software Fiji (Schindelin et al., 2012) and followed an automated macro (Fritz-Endres, 2022). The scans were concatenated into a 3D image block, cropped, and the background was subtracted by applying a rolling radius of 142 μm (Fritz-Endres, 2022). A 3-D Object Counter plugin was utilized to identify and map a projected area for each specimen and generate a 2D surface area measurement (μm²). Then each foraminifera was individually processed to generate a greyscale histogram and a 3D surface. Specimen 3D images are available at https://github.com/Foraminarium/Globigerina-bulloides/tree/main.

To examine trace element geochemistry, all samples were analyzed on a laser ablation system comprising a Photon Machines 193nm ArF laser with an ANU HelEx dual-volume laser ablation cell coupled to an iCAP quadrupole ICP-MS at the College of Earth, Ocean and Atmospheric Sciences at Oregon State University. Each of the foraminiferal chambers in the final whorl was ablated at least once, and larger chambers (typically F0, F1) were ablated twice using a laser spot size of 38 μm, a repetition rate of seven Hz, and a fluence of ∼1.0–1.2 J/cm², following recommended guidelines (Fehrenbacher et al., 2015). Trace element (TE) analytes included ²⁴Mg, ²⁵Mg, ⁸⁸Sr, ¹¹B, ¹³⁸Ba, ⁶⁶Zn, ⁴³Ca, and ²⁷Al. Because the tests were ablated from the outside to the inside, intratest TE/Ca variability is not explored in this manuscript. Fully resolving intratest variability typically requires ablating the test from the inside to the outside, starting from the flat inner test surface. Data were processed using LA-tools data reduction software, including screening for outliers, correcting background values, and normalizing to known concentrations with trace-element glass standards NIST-610 and NIST-612 (Branson et al., 2019). We applied a maximum aluminum (²⁷Al/⁴³Ca) threshold of 200 mmol to remove contaminant phases (Branson et al., 2019).

All statistics were completed using R (R Core Team, 2023). To test for normality between groups, a Shapiro-Wilk test was performed between the control and extraction group at each analysis step and between trace element analytes. A Welch’s two sample t-test assessed differences between groups if data were normally distributed. A Mann-Whitney U test was used when data were nonparametric.

Twenty-eight of the 30 extraction samples successfully amplified and sequenced, in addition to five other G. bulloides samples from Catalina Island and three samples from the Northeast Pacific. All extracted samples from Catalina Island were the same genotype (type IId), including a few samples where tests were lost or broken in the extraction process (Table 2). All type IId samples, including reference sequences and two Northeast Pacific sequences, were nearly identical (>98% matching identity; Fig. 3). The type IId G. bulloides samples separated from other types at >99% confidence level (Fig. 3). Therefore, subsequent analysis assumed extraction and control groups were genotype IId. This finding was expected because type IId is the most common G. bulloides genotype in the region and may be endemic to the area (Darling & Wade, 2008).

Phylogenetic analysis revealed that the G. bulloides sequences included in this study, including reference sequences, were similar (>78% matching identity; Fig. 3). The rarer regional subpolar genotype IIa was not found in this study (Darling et al., 2003). This result is not surprising because we sampled during the summer season and temperatures were warmer than average due to a marine heatwave (Harvey et al., 2023). Of the other North Pacific samples, two were also type IId and genetically similar to the Catalina Island samples (Fig. 3). We sampled only one different genotype—a single G. bulloides specimen collected at 45°N during an upwelling period—which was identified as the subpolar type IIe (Table 1; Fig. 3). This sample, SPR18, had a 100% matching identity with a reference type IIe sequence (Fig. 3). This type is the only genotype identified in the North Pacific subpolar gyre to date, so we would expect to encounter it at this latitude (Darling et al., 2007). In future studies in this region, we aim to expand the number of other G. bulloides genotypes sampled.

Microscope and SEM images revealed variable morphologies within the same genotype (type IId) between Catalina Island samples (Figs. 2, 4). Representative control and extraction group samples show variations in coiling direction and chamber number, but all present the expected G. bulloides morphotype with four globular chambers in the outer whirl and a wide aperture (Schiebel & Hemleben, 2017). Three additional Northeast Pacific samples were included for comparison, including another sample from Catalina Island, CAT02, collected a few days prior to the extraction experiment; and two samples collected from the Northern California Current; a type IId G. bulloides sample, NCC Y1 and genotype IIe sample, SPR18 (Fig. 4). The CAT02 sample was identified as G. bulloides genotype IId, but its morphology is not consistent with the morphotype. Environmental or laboratory contamination is possible. The test contains a fifth chamber which is opening slightly to the side (Figs. 4), morphologically more similar to Globigerinella calida (Parker, 1962). The pairing of morphometric and molecular analysis on the same test can thus improve our understanding of species morphologies. 3D images for each test provide further training tools and information for fellow researchers on morphological variability within and between genotypes (Fig. 5).

Previous research has shown that DNA extraction does not alter test physical properties (Ujiié et al., 2019). In our study, test mass was slightly lower in the extraction group, with specimens in the control group weighing 5.5 ± 1.4 μg and specimens in the extraction group weighing 4.6 ± 1.4 μg (t-test, p = 0.04; Table 2). Although this difference is statistically significant, we note that during DNA extraction, cleaning, and subsequent processing, some tests lost thin, final chambers. Remnants of broken chambers could be seen in the SEM images and 3D tomography (see C04, G31, and G39; Fig. 6). Furthermore, samples were collected live and had variability in the thickness of the final chamber, although all were from the same size fraction (>250 μm). Other studies have found that younger foraminifera have different final chamber thicknesses which will impact various test metrics compared to fully adult, gametogenic foraminifera specimens, such as those collected in sediment traps or core top sediments (Osborne et al., 2016).

Multiple studies suggest normalizing test weights by test size or area to account for differences in calcite thickness between similar size fraction test, particularly for G. bulloides (Barker & Elderfield, 2002; Marshall et al., 2013; Osborne et al., 2016). Area density, which accounts for test size and weight, was similar between the two groups (t-test, p = 0.1; Fig. 6a). Additionally, our study’s area density (0.12–0.21 × 10³ μg/μm²) was within the same range as other G. bulloides studies using XMCT estimates (Iwasaki et al., 2015, 2019a). Calcification intensity, as measured by the mean CT number, was similar between the control and extraction groups (573.5 ± 26 versus 564.5 ± 27, t-test, p = 0.2; Table 2; Fig. 6b).

Of the six elemental ratios, the two most widely used ratios (i.e., Mg/Ca, Mann-Whitney U test, p = 0.1, Table 3, Fig. 7a; and B/Ca, Mann-Whitney U test, p = 0.3, Table 3, Fig. 7b) were not influenced by GITC* DNA extraction. However, some trace element ratios were significantly higher after DNA extraction, notably Sr/Ca, Ba/Ca, and Zn/Ca (Mann-Whitney U test, p > 0.05; Table 3, Figs. 7d–f). Elevated Al/Ca ratios in the extraction group, though still below the threshold filter, indicate that some contamination remains on extracted samples (Mann-Whitney U test, p > 0.5; Fig. 7c).

We propose contaminant phases may have adhered to the test. The single oxidative cleaning treatment typically used for live-caught foraminifera tests was not sufficient for all analytes. We opted for the oxidative cleaning process because previous cleaning studies have recommended this protocol for planktonic foraminifera and found that further cleaning steps can remove the trace element elevated inner calcite (Fritz-Endres & Fehrenbacher, 2021). Although the reagents for GITC* buffer do not include trace metals, they are molecular grade and not trace metal free. One of the buffer reagents, Tris buffer, has been shown to form metal complexes, so trace element contamination could have been introduced during extraction and was not removed with a brief, oxidative clean (Fischer et al., 1979). In the next steps for this workflow, we propose including a quick, surface acid leaching step to remove any contaminants using a 0.001N HNO₃ solution (Boyle & Rosenthal, 1996). Further study would need to validate that this leaching step removes contaminants, but not inner calcite (Fritz-Endres & Fehrenbacher, 2021). However, for the widely used Mg/Ca and B/Ca proxies, this workflow is validated and promising.

Furthermore, this workflow could include stable isotope analysis, which could elucidate the known stable isotope variability in G. bulloides.Sadekov et al. (2016) proposed that the genotypic variability they observed in the Arabian Sea could be the cause of large bimodal variability in their oxygen isotopes. Specifically, their δ¹⁸O values suggest habitat temperatures between 10° and 28°C, where 10°C would imply a calcification depth of ∼700 meters, much lower than the known habitat range of G. bulloides (cf., near-surface conditions associated with the chlorophyll maximum in the Southern California Bight; Field, 2004). While Sadekov et al. (2016) collected genetic information from different specimens than the specimens they analyzed for geochemical information, our workflow would allow genetic and geochemical analysis to be applied to the same test and could help to resolve this conundrum. Although we did not collect such data, our workflow would allow tests to be analyzed for stable isotopes after LA-ICP-MS (Fig. 5). Previous research has shown DNA extraction does not influence stable isotope geochemistry (Ujiié et al., 2019).

We advanced and validated a workflow for highly resolved individual foraminiferal analysis including trace element geochemistry. Results yield genetic, morphometric, and geochemical data on G. bulloides, a species commonly used in paleoceanographic studies with high cryptic diversity. While Mg/Ca and B/Ca ratios are not impacted by DNA extraction, we found elevated Sr/Ca, Ba/Ca, and Zn/Ca ratios after extraction. More research is needed to explore the cause of this observation, but our study bears promise for elucidating the morphological and geochemical features of foraminiferal morphospecies with known or suspected cryptic diversity in a given region. While we only determined the genotypes of our specimens, the remaining DNA material could be used for further analyses, including the study of the microbial community associated with each specimen (i.e., its microbiome). This workflow could help illuminate the life history, population variability, and ecology of planktic species, and their respective impact on test geochemistry.

We would like to acknowledge the Catalina 2022 Culturing Crew, including Ingrid Izaguirre, Yoon Kim, Elise Poniatowski, Cate Fehrenbacher, and the staff of the USC Dornsife Wrigley Marine Science Center. We thank Theresa Fritz-Endres and Dorthe Wildenschild for assistance with processing microCT data. We thank Raphael Morard for help with molecular data. We acknowledge our funding sources, including the Cushman Foundation for Foraminiferal Research Loeblich and Tappan Student Research Award (MKL), the ARCS Foundation (MKL), the NSF Graduate Research Fellowship Program (MKL), NSF OCE21-03461 (BH, LH), and NSF OCE17-37165 (JSF, MKL).