Oceanic-Scale Species Diversity of Living Benthic Foraminifera: Insights into Neogene Diversity, Community Structure, Species Duration, and Biogeography

Benthic foraminifera are abundant and speciose members of the meiofauna in all marine environments from marshes and bays to abyssal depths and are important in marine ecosystem functioning. Moreover, they have been so for millions of years. Consequently, they are ideally suited to record diversity patterns of modern oceans as well as those of the past.

When an individual of the living benthic foraminiferal population dies or reproduces the empty test is often preserved in the sediment and becomes part of the dead population. Because the dead population is more abundant, the total population (live + dead) usually resembles the dead population. Over time the dead population becomes the fossil population. Researchers, of course, hope that the transition from living to dead to fossil population faithfully records the structure and composition of formerly living foraminiferal communities.

Hessler & Sanders (1967) demonstrated that within-habitat (α) diversity of the macrofauna in the deep-sea was as high as in the shallower depths of the tropics. Buzas & Gibson (1969) also found high within-habitat (α) diversity for the total population of meiofaunal foraminifera at abyssal depths along the Gay Head to Bermuda transect in the North Atlantic. A related finding is that the latitudinal diversity gradient, the trend of increasing diversity with decreasing latitude (Fisher, 1960), is apparent not only on the shelf (e.g., Culver & Buzas, 1998; Dorst & Schönfeld, 2013; Jablonski et al., 2017) but also in certain deep-sea benthic macrofaunal and meiofaunal groups (e.g., Rex et al., 1993, 1997) including the benthic foraminifera (Culver & Buzas, 2000; Dorst & Schönfeld, 2013).

The patterns noted above for benthic foraminifera were based on dead or total populations in the Atlantic Ocean basin (Buzas & Gibson, 1969; Culver & Buzas, 2000). The present study examines benthic foraminiferal diversity gradients in the Atlantic Ocean basin, from shelf, slope, and abyssal depths, using data sets of living populations (rose Bengal-stained; Walton, 1952) from an extensive compilation by Murray (2015) that was used by Jones & Murray (2017) in their analyses of standing crop (density) of benthic foraminifera on an oceanic scale. If the same patterns are found for the living population as formerly found for the total population, then we can be confident of our diversity assessment throughout the Neogene (Miocene, Pliocene) and Quaternary to the present (e.g., Thomas & Gooday, 1996; Culver & Buzas, 2000).

In this article, we: 1) investigate whether a latitudinal species diversity gradient exists at shelf, slope, and abyssal depths for within-habitat living populations of benthic foraminifera; 2) analyze differences in within-habitat diversity over shelf, slope, and abyssal depths for the living population; 3) compare the within-habitat diversity patterns exhibited by living and dead or total populations; and 4) integrate these data with previous studies of benthic foraminiferal community structure, species durations, and biogeography.

This study uses a subset of the data used most recently by Jones & Murray (2017) in their analysis of standing crop (density) values from the Atlantic basin. Jones & Murray (2017) extracted density data from a larger dataset compiled by Murray. This larger dataset was published in full in the electronic supplement to a book (Murray, 2006) and summarized in Murray (2015). This dataset included 2423 samples grouped by study (Murray, 2015) of which 1167 included counts of live specimens for each species encountered. The data from the Murray (2006) supplement (presented as multiple Excel workbooks) were collated into a single species-by-sample matrix (of count data) using the R programming environment (R Core Team, 2020). Each study encountered different species and in some cases used differing species names, either as a result of taxonomic revision or use of grouped species names (e.g., Ammonia group) or various taxonomic qualifiers (e.g., aff., spp., ?). As such, a manual quality control of the species names was undertaken by Murray to ensure that each species was represented by a single name in the final matrix. The cleaned dataset included a total of 1227 distinct “species”. Grouped species counts (e.g., unidentified agglutinated), if used, were counted as a single species in analysis, potentially leading to some underestimation of species totals.

The metric chosen for analysis of diversity is the Shannon (1948) information function, because this function includes not only species richness, but also species proportions. Single samples of sediment, each normalized to 10 ml, were analyzed, and each is considered to represent a foraminiferal habitat. Consequently, this study is about within-habitat or α diversity and not of regional (or γ diversity) or between-habitat (or β diversity; Whitaker, 1972). For a summary of the geologic, paleoceanographic, and paleoclimatic utility of benthic foraminifera, their important role in marine ecosystem functioning, see Gooday et al. (2008), and for detail on the dataset analyzed as part of this paper, see Murray (2015) and Jones & Murray (2017).

where p_i is the proportion of the ith species. Reasonable estimates of species richness (S) as well as species proportions (p_i) are required for the calculation of H. In temperate areas about 200 to 500 individuals are required for the species effort curve (plot of accumulated S vs accumulated N) to become asymptotic (Hayek & Buzas, 2010). In tropical shelf areas, the species effort curve shows no sign of abatement even when thousands of individuals are accumulated (Buzas et al., 1977), but a representative estimate of species proportions (p_i) is obtained by using 200 to 400 individuals (Hayek & Buzas, 2010). Consequently, we chose N = 200 as a minimum number of specimens counted in a sample as the criterion for inclusion in this study. Of the 1167 samples in the compiled matrix, 411 met this criterion (Fig. 1 ).

To compare means we used the Analysis of Variance (ANOVA). We choose to reject the null hypothesis when p < 0.05. Levene's test for homogeneity of variance was applied in each case where relevant and each was non-significant.

Murray (2007) estimated that the number of hard-shelled modern species of benthic foraminifera to be ∼3,200 to 4,200. However, this number ignored the many rare species. The WoRMS database (https://www.marinespecies.org) currently lists 8,953 recognized and named Recent species (Hayward et al., 2020). Further, Gooday (2019) noted that Murray's estimates do not include the many undescribed, single-chambered, soft-bodied (monothalamous) forms nor the “huge diversity” of unknown phylotypes (Lecroq et al., 2011). Delicate, loosely agglutinated tests are underrepresented in typical samples of benthic foraminifera owing to destruction during sampling and processing (Gooday et al., 1998) and, of course, below the calcium carbonate compensation depth (CCD), hard-shelled populations are dominated by agglutinated species (Gooday et al., 2008). Therefore, the benthic foraminifera investigated in the present study are hard-shelled fossilizable species from, in large part, the continental shelf and slope and the immediately adjacent abyssal plain (above the CCD) of the Atlantic Ocean basin. These are the taxa preserved in the fossil record and, consequently, their patterns of species diversity are of importance for understanding both modern ecosystem (Gooday et al., 1992) and paleoecosystem (Thomas & Gooday, 1996) functioning. The samples utilized in this study are derived from many data sets collected over six decades using several sieve sizes (>63 µm, >106 µm, >125 µm, and > 150 µm). Jones & Murray (2017) discussed at length the potential influence of this methodological variation on standing crop. When size-fraction was included in statistical models as a covariate, it was not significant in explaining the standing crop. They concluded that between sample density variation is a result of environmental variation rather than the size-fraction used. Density variation due to the former is orders of magnitude larger than density variation due to the latter.

For depths <200 m (continental shelf) in the northern hemisphere, 158 samples met the criterion of counts >200 individuals (Table 1). Results of a least squares regression analysis are shown in Table 2. In the 200–2000 m depth category (continental slope), 171 samples with N > 200 are distributed over both the northern and southern hemispheres (Table 3). Results of two least squares regression analyses are shown in Tables 4 and 5. For the >2000 m category (abyssal plain) in the southern hemisphere, 82 samples met the criterion of counts >200 (Table 6). Regression analysis results are shown in Table 7. In all depth zones and hemispheres assessed, there was widespread variability but an overall significant trend of decreasing within-habitat diversity of live foraminiferal populations with increasing latitude (Tables 2, 4, 5, and 7; Fig. 2).

The relationships are not as clear in the southern hemisphere slope (Table 5) and the abyss (Table 7) as in other areas. Mean values for H indicate the shelf (Table 1; Fig. 2) and slope (Table 3; Fig 2) have more variability along the latitudinal gradient than the abyss (Table 6; Fig 2). The difference between maximum and minimum values (the range) of H on the shelf is 2.59 (Table 1), on the slope the range is 1.59 (Table 3), and for the abyss is 0.25 (Table 6).

The samples analyzed in this study are arranged into three depth categories - shelf, slope, and abyss (Table 8). Figure 3 indicates a significant and striking increase in mean H with depth (Tables 8, 9) with mean values of H for <200 m, 200–2000 m and >2000 m of 1.94, 2.40, and 3.13, respectively.

Live, total and fossil populations exhibit the same patterns of within-habitat benthic foraminiferal diversity despite seasonality and relative rarity of live specimens and differences in population density or even presence of individual species in live populations owing to a variety of taphonomic circumstances (Murray, 1982; Mackensen et al., 1990). This characteristic of populations encourages the following discussion where we integrate the new data of this paper with published results of studies based on total and fossil populations. In this way, we can address the relevance and significance of oceanic-scale within-habitat diversity to biogeography, species durations, and community structure.

The live foraminiferal data of this paper indicate a significant increase in mean H with depth (Fig. 3). For the same depth categories, shelf (<200 m), slope (200–2000 m), abyss (>2000 m), Jones & Murray (2017) obtained mean density values of 237.4, 199.3, and 64.2 foraminifera per 10 ml, respectively, for live populations. Different sieve sizes were used by researchers and may have introduced bias into the results. However, as Tables 1 and 6 indicate, most of the shelf sieve sizes were 63 µm and all from the abyss were 125 µm. Consequently, any bias would result in underestimating the difference in values of H between depths. The increase in within-habitat diversity with depth in live populations agrees with the pattern of increasing diversity with depth in total populations from 350 samples ranging in depth from 29 m to 5,001 m and extending from the Arctic to the Gulf of Mexico (Buzas & Gibson, 1969; Gibson & Buzas, 1973). In that survey, maximum values of H occurred in samples from abyssal depths. Culver & Buzas (2000) demonstrated a latitudinal diversity gradient for the total population at abyssal depth in both hemispheres of the Atlantic, while Dorst & Schönfeld (2013) noted a similar pattern of diversity increase on the Atlantic shelf and slope off western Europe.

The number of benthic foraminiferal biogeographic entities, provinces, and their component communities, recognized by numerical and statistical analyses, decreases with depth. This pattern has been observed in the western Atlantic Margin of North America (Buzas & Culver, 1980), the Gulf of Mexico (Culver & Buzas, 1981), the Pacific continental margin of North America (Buzas & Culver, 1990) and New Zealand (Hayward et al., 2010). This is because deeper-dwelling species are more widespread and, hence, although the within-habitat diversity may be greater in the abyss, the total number of species is smaller than in shallower areas (Buzas et al., 2014).

The widespread distribution of deep-dwelling versus shallow-dwelling species is also supported by molecular studies. Hayward et al. (2021) showed that the three species of the shallow water genus Ammonia thought to be world-wide in their distribution actually belong to 60 species each with a limited distribution. In marked contrast, molecular studies on a cosmopolitan abyssal-dwelling species (Epistominella exigua) indicate genetic homogeneity across regions of the Arctic, Atlantic, Pacific, and Antarctic Oceans (Lecroq et al., 2009). The widespread abyssal species, Cibicidoides wuellerstorfi, does, however, exhibit some genetic differentiation between different areas (Burkett et al., 2020).

We noted earlier that shelf diversity data are from the northern hemisphere, abyssal data are from the southern hemisphere, and slope data are from both hemispheres. Thus, we can compare hemispheres for the latter data only. The latitudinal diversity gradient is greater in the northern hemisphere (Fig. 2B) than the southern (Fig. 2C). A weaker southern hemisphere latitudinal gradient also characterizes the deep-sea macrofauna, reflecting a higher degree of regional variation in the south (Rex & Etter, 2010).

The latitudinal gradient for within-habitat diversity that we see today at all ocean depths has a long history. For abyssal depths, Thomas & Gooday (1996) suggested the pattern for increasing diversity with decreasing latitude in benthic foraminifera originated at the Eocene-Oligocene boundary ∼34 Ma when the Earth transitioned from “greenhouse” to “ice-house” conditions. Neogene to modern benthic foraminiferal populations from shelf environments of the temperate Atlantic Coastal Plain and the tropical Central American Isthmus indicate that not only has a latitudinal gradient of diversity (measured by Fisher's alpha) been present for at least 10 Ma, but also that it has been increasing over time (Buzas et al., 2002a), by 40% at the temperate region and by 106% at the tropical region.

Species durations of benthic foraminifera (Buzas & Culver, 1984) show the same depth and latitudinal patterns as species diversity. Off the Atlantic coast of North America both partial durations (of living species) and species diversity are greater at lower latitudes and increased water depth: 16 Ma for <200 m (shelf) compared with 26 Ma for >200 m (slope and abyss), and 7 Ma for <200 m Cape Hatteras to Newfoundland compared with 20 Ma for <200 m Florida to Cape Hatteras (Buzas & Culver, 1984). Similar patterns of durations and diversity were documented around New Zealand (Hayward et al., 2010).

Also relatable to depth (shelf, slope, abyss) and, hence, diversity, is community structure, defined quantitatively by Buzas & Hayek (2011) and Hayek et al. (2019) as the mathematical statistical distribution fit to the observed relative abundance vector. Consideration of the decomposition equation for species richness (S), evenness (E), and H, plus their respective regressions on the accumulation of the number of individuals (N), leads to the establishment of three structural types of community. The types can be identified by a measure composed of the slope (β_1H) of the regression of accumulated H vs N within a community. A positive measure denotes community growth, zero denotes the existence of stasis, and a negative measure denotes the existence of a declining community. Global analysis of 72 communities with living and total populations were surveyed (Buzas & Hayek, 2011). For shelf and slope communities, the measure is either mostly positive or 0, while in the abyss either 0 or negative. The average measure is 0.13 for the shelf, 0.14 for the slope, and -0.06 for the abyss (table 21 in Buzas & Hayek, 2011).

In summary, the variables considered above and their relative values (extracted from the new data of the current study and from related earlier studies on benthic foraminiferal distribution and diversity through time) are shown in the contrast between shallow (<200 m) and deep (>200 m) categories presented in Table 10 . The tabulation demonstrates that shallow and deep-dwelling benthic foraminiferal communities are easily discriminated.

Within-habitat diversity is achieved through the interplay of species origination and immigration and species extinction and emigration over time (Buzas & Culver, 1998). To achieve high diversity, a community must maintain a relatively low extinction rate. Species density and the plethora of abiotic and biotic variables that determine its value (Jones & Murray, 2017) are important only as end values. Very low population densities may lead to extinction of species, thereby lowering diversity. Very high densities of organisms may lead to competition among species, and if there is competition among community members for a limited resource, then competitive exclusion demands a reduction in diversity. Relatively low extinction rates, then, suggest low overall ecological extinction from changes in abiotic and biotic variables and low competition among community members to ensure high diversities over time. The time component may require millions of years, fostering longer species durations in high diversity areas (Buzas & Culver, 1984).

Researchers have offered a variety of explanations for observed patterns of the latitudinal and depth diversity gradients (e.g., Pianka, 1966; Rohde, 1992; Rex & Etter, 2010; Jablonski et al., 2017; Gagne et al., 2020), and many of them are not mutually exclusive. Pontarp et al. (2019) have argued that the lack of consensus regarding the underlying causes for a latitudinal diversity gradient is due to the “verbal nature” of hypotheses and the fact that observed patterns can have multiple explanations. They proposed mechanistic linking of eco-evolutionary processes (selection, dispersal, ecological drift, and speciation) to the diversity gradient to better understand the contributions of these processes.

The great variability in values of H with latitude in shallower (<200 m) areas shown in this study suggests a variety of drivers are likely responsible for individual values. Nevertheless, there is a significant trend of decreasing diversity with increasing latitude. Gagne et al. (2020) modeled global diversity for terrestrial and marine species. Their analysis for marine organisms (44,575 species) indicated maximum diversity in the tropics. Depth, water temperature, and sunlight were the principal drivers. Curiously, their data set showed a decrease in diversity with depth, a reflection, perhaps, of the many organisms involved or just a consideration of gamma diversity. We do not consider depth as an environmental variable but, along with latitude and longitude, an attribute that locates a sample in space. It is the change in environmental variables associated with depth that is of primary interest. The variables water temperature and sunlight are reasonable and in accordance with advocates of solar energy or primary production as the principal driver of diversity (e.g., Rohde, 1992). Our data are not extensive enough to address the question of whether there is a decrease in species richness in the marine realm near the equator (Chaudhary et al., 2016; Woolley at al., 2016) or whether this is due to a knowledge gap (Menegotto & Rangel, 2018).

While Jones & Murray (2017) found an overall positive relationship between benthic foraminiferal density and particulate organic carbon (POC), on the shelf it was negative, prompting them to suggest predation (not food) was limiting density on the shelf. If predation on foraminifera (Culver & Lipps, 2003) affects all members of the community equally then, effectively, it limits density so that there is no species competition among community members. The lack of competition as judged by foraminiferal species with a community reacting in concert (pulsating patches) in shallow water was noted by Buzas et al. (2002b). However, for a latitudinal gradient, predation would have to be more severe at high latitudes to obtain the observed pattern. The presence of the gradient on the slope and in the abyss suggests a “trickle-down” ecologic economy where shallower vicissitudes are transferred to the deep ocean. The decreasing difference between maximum and minimum values of H from shelf to slope to abyss support this idea. The suggestion that the larger variability of environmental variables at the higher latitudes (particularly particulate organic matter flux to the sea floor) is responsible for the pattern is attractive (Hessler & Sanders, 1967; Thomas & Gooday, 1996; Rex & Etter, 2010; Cordier et al., 2022). Relative lack of variability explains why diversity is high in shallow tropical settings and nearly uniformly high in the abyss.

Numerous authors of research on benthic foraminifera (e.g., Gooday, 1988; Jorissen et al., 1995; Schmiedl et al., 1997), and on deep-sea communities in general (e.g., Smith et al., 2008), agree that food is an important limiting variable in the deep-sea (Buessler et al., 2007). As might be expected, the density of foraminifera in the deep-sea is much lower than on the shelf and slope (Jones & Murray, 2017). Although we hypothesized that predation reduced densities so that competition was not important in shallower waters, the greater reduction of abyssal densities is evidently still not great enough to cause extinction (but see below the mid-Pleistocene extinction event of elongate benthic species; Hayward et al., 2012). The low abyssal densities have not resulted in competition for food. Perhaps, the inputs from the surface water are so irregular in time and space (Gooday, 1988) that no one community member can have an advantage, allowing many species to cohabit within a community. However, recall that the status of abyssal community structures is at stasis or in decline so that abyssal communities are continually on the brink of extinction. The long species durations of abyssal communities, however, indicate extinctions are rare (background rate of ∼2% myr¹ during the Cenozoic in benthic foraminifera; Hayward et al., 2012) but extinction events occur. For example, the extinction event in the abyssal foraminiferal biota at the Paleocene-Eocene Thermal Maximum (PETM), at 55.5 Ma (Bowen et al., 2015) was accompanied by a negative value signifying a declining community (Hayek et al., 2019). The extinction of many elongate benthic species (25% loss of deep-sea benthic diversity) in the late Pliocene to middle Pleistocene, mostly between 1.2 and 0.55 myrs ago, was likely caused by decrease of specific phytoplankton food flux during global cooling leading up to the mid-Pleistocene Climate Transition (Hayward et al., 2012). Foraminiferal species confined to the abyss have long species durations, and many abyssal species are also distributed on the slope suggesting migration into the abyss from shallower depths (Hayward et al., 2010; Buzas et al., 2014). Such migration occurs with the macrofauna (Rex et al., 2005) and at shallower depths with the foraminifera (Buzas & Culver, 2009).

New data on shelf, slope, and abyssal living benthic foraminifera in the Atlantic Ocean basin demonstrate a latitudinal gradient of within-habitat diversity with increase toward lower latitudes in all depth categories and an increase in diversity with depth regardless of latitude. Similar patterns are seen for dead and total (live plus dead) foraminiferal populations allowing integration of the new data with diversity, community structure, species duration, and biogeographic patterns of Neogene fossil benthic foraminifera. Surprisingly, while density at abyssal depths decreases owing to decreased food supply compared to the shelf and slope, within-habitat diversity is not affected and is high in the abyss.

We acknowledge the many authors who provided data to Professor John W. Murray. Dave Mallinson, Seth Sutton, and Laura de Sousa kindly provided technical support. Reviews by Andrew Gooday and Bruce Hayward were very useful and much appreciated. This work was supported by the NERC National Capability funding to the National Oceanography Centre, as part of the Climate Linked Atlantic Section Science (CLASS) program (Grant Reference NE/R015953/1).