Exceptional 20th Century Shifts in Deep-Sea Ecosystems Are Spatially Heterogeneous and Associated With Local Surface Ocean Variability

Traditionally, deep-sea ecosystems have been considered to be insulated from the effects of modern climate change, but with the recognition of the importance of food supply from the surface ocean and deep-sea currents to sustaining these systems, the potential for rapid response of benthic systems to climate change is gaining increasing attention. However, very few ecological time-series exist for the deep ocean covering the twentieth century. Benthic responses to past climate change have been well-documented using marine sediment cores on glacial-interglacial timescales, and ocean sediments have also begun to reveal that planktic species assemblages are already being influenced by global warming. Here, we use benthic foraminifera found in mid-latitude and subpolar North Atlantic sediment cores to show that, in locations beneath areas of major surface water change, benthic ecosystems have also changed significantly over the last ∼150 years. The maximum benthic response occurs in areas which have seen large changes in surface circulation, temperature, and/or productivity. We infer that the observed surface-deep ocean coupling is due to changes in the supply of organic matter exported from the surface ocean and delivered to the seafloor. The local-to-regional scale nature of these changes highlights that accurate projections of changes in deep-sea ecosystems will require (1) increased spatial coverage of deep-sea proxy records, and (2) models capable of adequately resolving these relatively small-scale oceanographic features.


INTRODUCTION
The deep-sea benthic habitat is one of the most extensive on Earth. It spans continental slopes, seamounts, deep-ocean trenches, and the abyssal plains. Habitat types range from soft sediment environments to rocky slopes, boulder fields, and hydrothermal vents. Habitat forming species such as cold-water corals and sponges add additional three-dimensional structure and provide nursery and shelter for many other species (Roberts et al., 2006;Howell et al., 2016). New species are frequently discovered during scientific expeditions.
When compared with the surface ocean environment, the deep ocean is relatively isolated from the habitat changes most typically associated with anthropogenic greenhouse gas emissions, i.e., temperature change, ocean acidification and changes in mixed layer depth. However, direct and rapid links between the surface and deep ocean are increasingly being recognized. The main source of food to organisms living in the deep ocean is organic matter sinking as a result of surface productivity, and the amount and type of this export productivity is a driver of the abundance and diversity of benthic species (Sun et al., 2006;Corliss et al., 2009;Gooday et al., 2012;Diz and Barker, 2016). The volume of this export productivity is controlled by numerous factors (temperature, positions of fronts, nutrient supply, mixed layer depth, size, and type of matter), which are all affected by climate change. In addition, the presence or absence, and strength, of deep-sea currents plays an active role in determining the bottom habitat, and in stimulating food supply for suspension feeders. The strength of such currents can respond quickly to wind and buoyancy forcing of the surface ocean. The deep ocean is therefore not immune to the effects of anthropogenic climate change.
However, deep-ocean ecosystems are difficult to monitor at spatial scales beyond the light penetration limit of submersibles, or over long time periods. With the exception of a few dedicated monitoring stations (Gooday et al., 2010) and rare high resolution sedimentary records (Dijkstra et al., 2015(Dijkstra et al., , 2017, observations of deep-sea ecosystems are made from research expeditions that rarely return to exactly the same locations. It has therefore been difficult to demonstrate that modern climate change is affecting deep-sea ecosystems beyond the range of natural variability. Over much longer timescales, deep-sea sedimentary records containing benthic foraminifera, as well as collections of ancient coral skeletons, have revealed that both types of organism have responded rapidly to climatic change during the last 100,000 years. For example, cold-water coral mounds expanded northward in the North Atlantic and species changed their ranges in the Drake Passage during the last deglaciation, keeping pace with deglacial climate change (Henry et al., 2014;Margolin et al., 2014). Benthic foraminifera inhabit a number of niches in deep-sea sediments, including epifaunal species sensitive to the presence of a bottom layer of detritus derived from surface productivity, and infaunal species sensitive to pore water oxygen concentration, itself a function of the concentration of organic matter within the sediment. Previous work has shown that benthic foraminifera responded to past rapid climate changes, such as abrupt events during the last glacial period (Thomas et al., 1995;Ohkushi et al., 2013;Moffitt et al., 2014Moffitt et al., , 2015, and to more modest climatic change during the Holocene (Ohkushi et al., 2013;Moffitt et al., 2014Moffitt et al., , 2015Perner et al., 2015;Smart et al., 2019;Palmer et al., 2020). Now, with a growing collection of high-resolution sediment archives that capture recent climate variability, we are able to explore changes in benthic ecosystems during the twentieth century and place them in the context of thousands of years of prior ocean history. Here, we present data on benthic foraminifera populations over the Holocene period from the North Atlantic, a region of the ocean that has been responding to industrial-era (∼1850 CE onward) climate change in a variety of ways, mainly via changes in the strength and positions of surface and deep ocean currents and associated changes in temperature and productivity (Thornalley et al., 2018;Osman et al., 2019;Spooner et al., 2020). The primary aim of this research is to undertake an initial investigation into whether the anomalous changes in the oceanography of the 20th century North Atlantic have altered benthic ecosystems, and if so, where. Building on the work of Spooner et al. (2020), which investigated 20th century changes in ecology and hydrology in the northeast Atlantic surface ocean, we examine benthic and planktic records from across the North Atlantic basin. We supplement the (mostly planktic) data of Spooner et al. (2020) with new planktic and/or benthic data from the northeast Atlantic (shallow Iceland Basin, the Reykjanes Ridge, and the Eirik Drift) and present new data from the northwest Atlantic (the Scotian Slope and the Laurentian Slope). In this basin-wide overview, we focus on the common features and relationships observed across our sites for the Holocene and the industrial era (past ∼10 kyr). We show that, in areas most affected by surface ocean changes (i.e., the northeast and western margin of the Atlantic close to major water mass transitions), benthic foraminifera underwent large 20th century changes. In some cases, and in line with certain surface ocean data, these changes appear exceptional in the last 10,000 years.

MATERIALS AND METHODS
Benthic foraminifera assemblages were analyzed from six locations spanning the North Atlantic between 35 • N and 60 • N. With the exception of Reykjanes Ridge, two cores were analyzed from each location: a gravity or piston core spanning thousands of years, and a companion multi-or box-core containing the sediment-water interface (Table 1 and Figure 1).
RAPiD-17-5P/EN539-MC16-A/B are three cores from within 100 m of one another from the Iceland Basin, within the path of the Iceland Scotland overflow water (ISOW; Moffá-Sanchez et al., 2015;Spooner et al., 2020). We shall refer to these cores as "IB deep." EW9302-29GGC and EN539-MC25-A are also from the Iceland Basin, but lie in shallower water closer to Iceland and are bathed by Eastern Mode Water (Hátún et al., 2016). We shall refer to these cores as "IB shallow." The core RAPiD-21-3K, "Reyk R, " is from the eastern Reykjanes Ridge and lies in downstream ISOW, where flow speed is modulated by Labrador Sea Water (Boessenkool et al., 2007). The cores RAPiD-35-14P/RAPiD-35-25B (Moffa-Sánchez et al., 2014;Moffá-Sanchez et al., 2015), "Eirik D, " are from the Eirik Drift, located in the central subpolar gyre within the path of Denmark Straits overflow water (DSOW). These cores have been spliced together in previous studies (RAPiD-35-COM;Moffá-Sanchez et al., 2015), although in our present study a separate sub-core of RAPiD-35-25B has been used: a slab taken at the time of the box-core retrieval, obtained by pressing a plastic tray used for kasten core sampling into the side of the box-core. The "slab" sub-core has a more highly resolved top few centimeters compared to the Moffa-Sánchez et al. (2014) and Moffá-Sanchez et al. (2015) sub-core whose uppermost portion underwent significant compaction during sub-coring and its storage in an upright vertical position for several years prior to sampling. The remaining cores are from the Scotian Slope and Laurentian Slope in the western North Atlantic, a region which has undergone significant change in upper ocean temperature and inferred circulation (Thornalley et al., 2018). KNR197-10-MC45-C and KNR197-10-44GGC, "Scot S, " are within the path of upper Labrador Sea Water (LSW) and KNR158-4-9GGC and KNR158-4-10MC, "Laurent S, " are within the path of LSW. Age models for IB deep, IB shallow, Eirik D, and Laurent S were taken from the literature cited in Table 1 ( Spooner et al., 2020). Given the uncompacted core top of our subcore for RAPiD-35-25B, and that the box-core sediment-water interface was collected intact, we assume a modern age for the sediment-water interface, and then apply the age constraints of Moffa-Sánchez et al. (2014) by transferring the age-depth tie points onto our subcore depths via the alignment of the common features in their percent coarse fraction records (Supplementary Figure 1). A new age model for Reyk R was constructed using 210 Pb and 14 C data from Boessenkool et al. (2007) and Sicre et al. (2011). 210 Pb data was taken from the companion core RAPiD-21-12B (Boessenkool et al., 2007). The percent coarse fraction records of both cores suggest that these data can be applied to Reyk R (Supplementary Figure 2; Miettinen et al., 2012;Moffa-Sánchez and Hall, 2017). 14 C data were calibrated using the MARINE13 calibration curve (Reimer et al., 2013), individually using Calib7.1 (Stuiver et al., 2021), and by using a Bayesian framework in OxCal with a Poisson distribution model (Ramsey, 2008). One date at 206-207 cm was excluded from the age model due to poor agreement. A new age model for KNR197-10-MC45-C was constructed using 210 Pb data measured on bulk sediment at University College London (Supplementary Figure 3). The age model for KNR197-10-44GGC was constructed using δ 18 O measurements measured on Neogloboquadrina pachyderma at the University of Cambridge and 210 Pb data measured on bulk sediment at University College London. N. pachyderma δ 18 O measurements were tuned to published N. pachyderma δ 18 O data from the previously radiocarbon-dated cores OCE326 14GGC and OCE326 26GGC (Supplementary Figure 4; Keigwin et al., 2005). The core top of KNR197-10-44GGC was shown to overlap with the base of KNR197-10-MC45-C using a combination of planktic foraminifera species counts and grain size data (Supplementary Figure 5).
Benthic foraminifera were typically counted in the >150 µm fraction (rather than >63 µm) because the cores sites are located at drift sites where the relatively fast flowing bottom currents have the potential to transport smaller benthic foraminifera. Counts on the >63 µm fraction were conducted at Reyk R, where sediment grain size analysis indicates the presence of only relatively fine silts (Supplementary Figure 6; Boessenkool et al., 2007). In EN539-MC25-A, higher resolution counts were also made in the >212 µm size fraction (Supplementary Figure 7). Choice of size fraction can impact the relative abundance of species found at different sites. For example, the epifaunal foraminifera Epistominella exigua is generally smaller than >150 µm, but is abundant in the 63-150 µm size fraction. Therefore, we compare general trends across core sites rather than comparing species assemblages, and our main arguments rely on the total accumulation of benthic foraminifera.  Table 1). (A) Colors indicate the mean annual climatological sea-surface temperature (SST; • C) from the World Ocean Atlas (Locarnini et al., 2013). Purple circles indicate A, OCE-326-30GGC/MC29 (Emerald Basin; Keigwin et al., 2003), and B, OCE-326-14GGC/MC13 (Laurentian Fan; Keigwin and Pickart, 1999). Furthermore, at site IB shallow, and more limited data at site Reyk R, we note that similar trends are observed for the total benthic abundance obtained from >63, >150, and >212 µm (Supplementary Figures 6-8), suggesting the robustness of these trends regardless of specific size fraction. Benthic foraminifera counts included calcifying species, and excluded agglutinating forms due to their typical poor preservation downcore, with the exception of the well-preserved Miliolida Sigmoilopsis schlumbergeri. We include planktic foraminifera counts from the literature or newly counted on the >150 µm size fraction to discern changes in the surface ocean. We provide new, unpublished counts in EW9302-29GGC, RAPiD-35-14P and 25B, Laurent S and Scot S.
Accumulation rates (AR) of benthic and planktic foraminifera were calculated according to the formula: where #/g is the number per wet gram of sediment, and ρ is the density of the wet sediment, calculated from the dry mass and mass of water assuming an average grain density of 2.65 g/cm 3 , and water density of 1.027 g/cm 3 . Water content data were not collected for EW9302-29GGC or KNR158-4-10MC, and we therefore estimated the water content in these cores from the base of EN539-MC25-A (50%) and the whole of KNR158-4-9GGC (50%), respectively. Water content for the samples taken from RAPID-35-25B was estimated using the samples in Moffa-Sánchez et al. (2014). The sedimentation rate at the top of EW9302-29GGC is complicated by a ∼5 cm thick ash layer between 6 and 11 cm. Therefore, despite the radiocarbon dates suggesting a very high sedimentation rate of 84 cm/kyr, we assign a sedimentation rate to the top sample (0-1 cm) of 27 cm/kyr, the mean sedimentation rate of the Holocene section of the core (±4.3 cm/kyr 2 SD).
Accumulation rates and % abundance data were plotted on linear time-series to examine temporal trends and were accompanied by density plots to more clearly evaluate the distribution of Holocene (10 kyr b2k to 1750 CE) versus industrial era (1900 CE onward) values in each dataset.

RESULTS AND DISCUSSION
Our results demonstrate that benthic foraminifera have been responding to ongoing climate change during the industrial era, and that the state reached by the benthic fluxes can be unusual for the Holocene. In addition, the regions that show the largest magnitude changes in the benthic accumulation rates are also subject to significant surface ocean changes involving the shifting properties/locations of major ocean currents (Figures 1, 2). Thus, the benthic response to climate change is not an en masse basinwide response, but is dependent on local surface ocean conditions and can therefore differ depending on location.

Northeast Atlantic
In the northeast Atlantic close to Iceland, we find a range of 20th century benthic responses. At IB deep, the abundance of benthic foraminifera declined during the 20th century and reached a Holocene low (0.2 cm −2 yr −1 ) during the last few decades (this study; Spooner et al., 2020). This result is compared to a Holocene maximum of up to 4 cm −2 yr −1 (Figures 3A,B). Epifaunal and infaunal genera show the same trends in constant relative abundance (Supplementary Figure 9), so the decline in accumulation rate cannot be explained by a core top expression of habitat zoning. The decline in abundance is not accompanied (E) Greenland ice core methanesulfonic acid (MSA) record of subpolar surface ocean productivity, orange (Osman et al., 2019); (F) AMOC index of Rahmstorf et al. (2015); (G) 5-year averaged subarctic Atlantic freshwater storage anomalies (±2σ; anomalies relative to 1955) from Curry and Mauritzen (2005); (H) Surface ocean temperature changes in the northwest Atlantic slope region based on % abundance of the polar species N. pachyderma in cores OCE-326-14GGC/MC13 [Laurentian Fan; (Keigwin and Pickart, 1999); age model updated to include a tie at 4.5 cm for the 1929 Grand Banks turbidite which deposited a red clay layer at this site], OCE-326-30GGC/MC29 [Emerald Basin; (Keigwin et al., 2003); using the published 14 C and 210 Pb data for the age model], and Laurent S (Laurentian Slope; Thornalley et al., 2018; and this study), for visualization of the relative temperature changes OCE-326-14GGC/MC13 has been offset by +10%; (I) Benthic foraminifera flux in Laurent S (this study); (J) Benthic foraminifera flux in Scot S (this study). Vertical gray-shaded bands indicate post-1850, 1950, and 1970 periods. *Note, for this figure an additional age tie was included for IB deep (20.75 cm = 1848 CE, decreasing the core age by 28 years, well within the age model uncertainty) based on alignment of the % T. quinqueloba record with the remarkably similar MSA subpolar productivity record, the latter benefiting from a more precise age model. by significant changes in the percentages of any major group of benthic foraminifera [e.g., high food group (Smart et al., 2019), phytodetritus group, epifaunal, infaunal, Miliolida] at either core site (Supplementary Figure 9), suggesting that selective taphonomic bias also does not affect our results (Stefanoudis et al., 2017). The abundance trend is similar to the collapse in planktic foraminifera abundance in the northern Iceland Basin over the 20th century (Figures 3F,G), as well as a change in the dominant planktic species (Figures 3K,L), suggesting declining export productivity at the eastern edge of the subpolar gyre (SPG). It has previously been argued that this change in productivity was caused by a shift in local frontal activity associated with the retreat of subpolar water and its replacement by transitional water with its origins in the subtropical gyre, as evidenced by the steep decline in the abundance of the subpolar and high productivity-associated planktic foraminifer Turborotalita quinqueloba (Figures 3K,L), and the increase in transitional-subtropical Orbulina universa (Figures 3P,Q) (Spooner et al., 2020). Slightly further north at site IB shallow, there is some evidence for a higher recent accumulation rate of benthic foraminifera, but this is not beyond other events seen in the Holocene (Figures 3C,D). In addition, the number of benthic foraminifera per gram of sediment does not show such a prominent late 20th increase as seen in IB deep (Supplementary Figure 16). Moreover, the observed slight increase in benthic foraminiferal accumulation rate at IB shallow does not exceed our uncertainty estimates. Overall, the evidence for a significant 20th century shift in benthic foraminifera at this site is not very robust. This finding reflects our earlier results that showed little change in the planktic foraminifera assemblage and abundance (most notably T. quinqueloba) over this time interval (Figures 3H,I,M,N,R,S; Spooner et al., 2020).
Approximately 500 km to the southwest, at Reyk R, the 20th century is characterized by an increase and then a decrease in Turborotalita quinqueloba % abundance; Orbulina universa % abundance (P-T). Turborotalita quinqueloba % abundance in IB deep and Reyk R, and Orbulina universa % abundance data in IB deep, are from Spooner et al. (2020). Benthic abundance data (benthics g -1 ) for IB deep and planktic abundance data (planktics g -1 ) for IB deep and Reyk R were previously published (Spooner et al., 2020) and have been converted to benthic and planktic fluxes, respectively (this study). Vertical red-shaded band indicates industrial era. Gray-shading in the benthic flux time series plots represents 2σ error.
both the benthic and planktic accumulation rates (Figures 3E,J). The benthic decline after around 1980 leads to the lowest benthic accumulation rate in the record (spanning 1500 years). The most abundant benthic foraminifera at Reyk R (>63 µm) was Epistominella exigua, an opportunistic species that dwells in the phytodetritus layer; however, the recent decline in benthic foraminifera is seen in infaunal and epifaunal species, and photodetritus and non-photodetritus associated species (Supplementary Figure 12). Although more muted than the surface ocean changes recorded at our sites closer to Iceland, after ∼1980, Reyk R also records a decline in the high-productivity, subpolar planktic species T. quinqueloba ( Figure 3O) (Spooner et al., 2020) and an increase in the warm, more-oligotrophic associated O. universa (Figure 3T), suggesting a moderate decrease in the local surface ocean productivity. Additionally, Reyk R lies downstream of IB deep in the Iceland-Scotland overflow water (ISOW), so it is possible that the effects of changes in surface oceanography and export in the northern Iceland Basin were propagated by deep-sea currents to the deep ecosystems of the Reykjanes Ridge, compounding any local decrease in the surface-to-deep food supply.

The Labrador Sea
In the central SPG (Eirik D, Figures 1,4), the range of variability in benthic foraminifera abundance over the last 500 years is not significantly different to that of the whole Holocene. There is also no difference in the mean benthic flux for the 20th century compared with the Holocene (Figure 4B). The cores from this location are not as exceptionally high sedimentation rate as those from the eastern north Atlantic, so it is possible FIGURE 4 | Ten thousand-year Labrador Sea sediment records and density distribution plots comparing the distribution of Holocene (black) versus industrial era (red) data from Eirik D (A,C,E) and HU-90-013-013 (Bilodeau et al., 1994). (A-C) Benthic foraminifera flux; (D,E) planktic foraminifera flux; (F,G) Turborotalita quinqueloba % abundance (black, panel F only) and Neogloboquadrina pachyderma % abundance (gray). Vertical red-shaded band indicates industrial era. Gray-shading in the benthic flux time series plots represents 2σ error. that some shorter term climate signals have been smoothed by bioturbation; however, the presence and exponential decline of excess 210 Pb in the top 5 cm and bomb radiocarbon in the core-top sample in the lower resolution subcore used by Moffa-Sánchez et al. (2014) suggest bioturbation was not too extensive; we estimate that the top 4 cm of our subcore captures the 20th century (Supplementary Figure 1). Holocene benthic foraminifera abundances at Eirik D are similar in magnitude and variability to benthic fluxes reported for nearby core HU-90-013-013 from the Greenland Rise (Bilodeau et al., 1994). The two youngest samples (∼1550 and 1750 C.E.) from HU-90-013-013 show unexpectedly high values compared with Eirik D, but this could be due to disturbance of the core top. Overall, the null hypothesis that benthic foraminifera have not significantly changed during the industrial era is consistent with the existing data.
Eirik D planktic foraminifera abundance (Figures 4C,D) and the relative abundances of key indicator species such as N. pachyderma (Figures 4E,F) show a slight decrease and increase, respectively, in the industrial era relative to the Holocene. The increase in % N. pachyderma began well before the onset of the industrial era, around 3000 years b2k ( Figure 4E). Overall, however, planktic foraminifera do not show significant industrial era shifts in these cores, suggesting only small changes in surface ocean conditions. In the northeast central Labrador Sea, the core site is not in a position where significant changes in water mass might be expected, being influenced by subpolar waters and the topographically constrained extension of the East Greenland Current. The site is located close to the socalled "warming hole" -the subpolar region where there has been an absence of recent warming (a relative cooling compared to the global mean SST) and which has been related to an industrial-era weakening of AMOC (Caesar et al., 2018); the lack of any recent warming is consistent with our observation of no major industrial-era changes in the planktic foraminifera assemblages. Furthermore, recent observational results indicate that the Labrador Sea contributes minimally to basin-wide overturning (Lozier et al., 2019;Petit et al., 2020), and convection in the Labrador Sea is not closely related to the basin meridional overturning circulation due to density compensation (Zou et al., 2020). Therefore changes in AMOC can be achieved without major shifts in the circulation and overturning of the Labrador Sea. Overall, we conclude that Eirik D has likely experienced relatively similar surface productivity conditions throughout the Holocene and industrial era.

Northwest Atlantic Slope
In the northwest Atlantic, at the Laurentian Slope site Laurent S, benthic foraminifera accumulation was quite stable for the last 8.5 ka at around 0.5 cm −2 yr −1 , but an increase in abundance up to typically 1-2 cm −2 yr −1 emerged during the late 19th century and 20th century (Figures 5A,B). This increase in overall abundance of benthic foraminifera is not associated with any clear change in the relative or absolute abundance of different species groups (Supplementary Figure 14), suggesting that the change is a broad response to increased food supply but remaining within the same type of benthic environment. Notably, a rapid increase in the relative abundance of the planktic foraminifera N. pachyderma from 25 to 50% occurred over the same time interval (Figures 5I,J).
Just slightly further south, at site Scot S on the Scotian Slope, benthic fluxes were low in the Holocene, ∼0.01-0.1 cm −2 yr −1 and then increased from the late 19th century to typical values of 1-2 cm −2 yr −1 for the 20th century (Figures 5C,D). Unlike Laurent S, planktic accumulation rates were higher in the 20th century, 2 cm −2 yr −1 , compared with the Holocene, <0.1 cm −2 yr −1 (Figures 5G,H). Also, % N. pachyderma abundances were variable in the Holocene and decreased in the 20th century (Figures 5K,L).
The elevated 20th century benthic foraminifera accumulation rates at both the Laurentian and Scotian Slope sites (Figures 5A-D), which experienced contrasting surface ocean changes (Figures 5E-L), may be related to inferred changes in ocean frontal gradients/activity. In general, the northwest Atlantic slope region is characterized by northeastward flowing warm slope water influenced by the Gulf Stream, and southwestward flowing cold Labrador slope water. The confluence of these two water masses results in sharp surface ocean thermal gradients and strong frontal activity -in the region of our sites termed the Slope Water jet/front ) -resulting in highly productive waters and enhanced food supply to the benthos (Belkin et al., 2009). Evidence presented here, and in previous work (Thornalley et al., 2018, Extended Data Figure 5), reveals that after the late 1800s, there was a greater proportion of cold Labrador slope water on the Laurentian Slope at site Laurent S (Figures 5I,J, 2H), while sites located further south and west [Laurentian Fan sites, site Scot S (Figures 5K,L), and Emerald Basin (Figure 2H)], where there is greater dominance of warm Gulf Stream derived waters, show warmer watersthe opposite signal. This pattern implies that in comparison to earlier centuries, the 20th century is associated with increased temperature gradients across the region, which will have resulted in enhanced surface ocean frontal activity. Here, we suggest that this inferred frontal activity likely promoted increased export productivity from the surface ocean to the deep sea and thus the observed greater benthic productivity at our slope sites.
Although the northwest Atlantic shelf and slope region is an area with complex ocean circulation, we suggest that the broad regional oceanographic changes described above may be related to previously reported weakening of both the AMOC (Thornalley et al., 2018) and the SPG (Spooner et al., 2020). AMOC weakening since ca. 1850 ( Figure 2F) is thought to cause a northward shift of the Gulf Stream, causing warming at the southwest sites and cooling at the Laurentian Slope site (Thornalley et al., 2018). Furthermore, weakening of the SPG has also been associated with greater penetration of cold subpolar water into the Newfoundland basin, which may also help explain the cooling observed at our open-ocean Laurentian Slope site (Bersch, 2002;Hátún et al., 2005). Earlier work by Pickart et al. (1999) identified a similar mode of variability operating on shorter timescales (interannual-decadal) at the Scotian Slope, characterized by a northward shift of the Slope Water front at the same time as a strengthening of the off-shelf, southwestward flowing, cold Labrador Current and the offshore, northeastward flowing, warm Slope Water jet.
It is important to note that the benthic ecosystem changes we observe at our offshore, open-ocean, relatively deep slope sites (∼1 and 1.8 km water depth) should be considered distinct from reported changes occurring on the shallow continental shelf, including the Gulf of Maine, where there are different ocean currents and dynamics [notably the inshore shelf branch(es) of the Labrador Current]. Moreover, on the shelf, benthic ecosystems can be controlled directly by changes in the upper ocean circulation that bathes the seafloor at the shallow water depths, with a reported recent warming of the shelf causing substantial changes in shallow marine ecosystems (Pershing et al., 2015;Kavanaugh et al., 2017).

Timing and Cause of the Benthic Responses to Industrial Era Climate Change
Our results demonstrate that North Atlantic benthic foraminifera exhibited regional responses to changes in climate and ocean circulation during the industrial era. However, for locations where benthic fluxes show significant shifts in the industrial era compared to the Holocene, the timing of these changes are not always synchronous between sites (Figure 2). Our ability to link the timing of benthic foraminifera changes to specific changes in North Atlantic ocean circulation is hindered by limitations in our knowledge of when, what and how circulation changed during the industrial era. Notwithstanding these limitations, and by greatly simplifying the evidence, existing proxy data suggest that there have been two major industrial era shifts, the first (Phase 1) occurring at ∼1850-1900, and then a series of shifts occurring through the mid-late 20th century (Phase 2) with evidence for more marked changes at ∼1970 onward in some datasets, which may be linked to increased freshwater in the subarctic Atlantic (Figure 2; Rahmstorf et al., 2015;Moore et al., 2017;Caesar et al., 2018Caesar et al., , 2021Thornalley et al., 2018;Osman et al., 2019;Lower-Spies et al., 2020;Spooner et al., 2020). There has also been large amplitude recent decadal variability through the 1990s and early 21st century (Robson et al., 2016;Holliday et al., 2020). The first shift (Phase 1) has been associated with weakening of LSW flow in the Deep Western Boundary Current [DWBC LSW (Thornalley et al., 2018); please note that recent work has suggested that density changes in the Labrador Sea may be governed by "upstream" processes occurring in the subpolar NE Atlantic (Lozier et al., 2019;Zou et al., 2020)], a weakening of the AMOC, a westward contraction and weakening of the SPG, and exceptional shifts in the surface circulation of the northwest Atlantic slope and shelf. The late 20th century changes (Phase 2) appear most pronounced in records linked to the upper ocean subpolar gyre, while little-to-no change is recorded in the DWBC LSW data, the off-shelf surface temperature records of the northwest Atlantic or the Zhang subsurface AMOC temperature fingerprint index (Thornalley et al., 2018; Figure 2). Similarly, previous work has also noted that the spatial pattern of temperature change in the subpolar North Atlantic over recent decades differs from the longer-term, century trend (Spooner et al., 2020). Therefore, there are likely different underlying mechanisms and circulation patterns associated with the two major industrial era shifts as well as recent decadal variability, that are not yet fully resolved or understood. Figure 2A shows that at IB deep in the northeast Atlantic, there was only a modest decrease in benthic foraminifera flux during the early stage of the industrial era, but a steeper decline occurred after 1950 (Phase 2), which has been linked to changes in SPG circulation (Spooner et al., 2020). At Reyk R, also in the northeast Atlantic, benthic fluxes increased slightly during the second half of the 19th century (Phase 1), remained stable for most of the 20th century, before declining sharply during Phase 2 from the 1980s onward ( Figure 2B). The recent decline coincides with the more extreme changes observed at IB deep and the enhanced incursion of warm oligotrophic waters into the Iceland Basin ( Figure 2C). This may be interpreted as a threshold response at Reyk R, such that only the more extreme recent conditions led to a significant change at this site. The contrasting benthic shifts at Reyk R between Phase 1 (increase) and late Phase 2 (decrease) may also be related to the underlying circulation differences between these two phases, noting specifically that the 1980s through to the early 2000s was associated with strong basin-wide warming, whereas Phase 1 had eastern subpolar surface warming but central subpolar cooling (Spooner et al., 2020). Labrador Sea (Eirik D) and the Iceland Slope (IB shallow) north of the SPG show no clear signal during the recent most decades. Taken together, this suggests that the major changes in the central Iceland Basin (IB deep) and along the eastern flank of the Reykjanes Ridge (Reyk R) are not due to broad changes in the SPG water masses, nor in the subpolar mode waters around the eastern periphery of this gyre, but instead reflect shifts in a major frontal system.
Both northwest Atlantic slope sites (Scot S and Laurent S) record a sharp rise in benthic foraminifera abundance during Phase 1 (∼1850-1900; Figures 2I,J), which we have linked to the concomitant exceptional regional surface ocean changes that resulted in enhanced thermal gradients and inferred frontal activity. It is hypothesized that the underlying driver for this regional circulation shift was weakening of the AMOC and/or SPG, although continued work is required to elucidate the details of circulation changes in this region and their link to basin-wide circulation features. At Scotian Slope site Scot S there is a decline in benthic foraminifera abundance from ∼1930 onward (Phase 2), the cause of which remains uncertain, although it may relate to the documented increasing dominance of warm slope water in the region and its impact on the regional ecosystem (Sherwood et al., 2011). At Laurentian Slope site Laurent S there is no major trend or shift in benthic foraminifera abundance during Phase 2 (Figure 2I), mirroring the surface ocean records (Figure 2H), although the core does not span the most recent decades.

Assumptions in the Interpretation of Benthic Foraminifera Abundance Data
Our main findings are based on the assumption that benthic foraminifera accumulation rates are primarily controlled by food supply from the surface ocean, and covary with surface export productivity (Sun et al., 2006;Jorissen et al., 2007;Corliss et al., 2009;Gooday et al., 2012;Diz and Barker, 2016). However, the amount of organic material reaching the seafloor also depends upon the efficiency of the surface-to-deep organic matter flux, which can vary according to the size of organic matter and its packaging into larger aggregates. Benthic foraminifera may also respond more strongly to seasonal blooms of export of food to the seafloor, as opposed to the integrated annual flux, and certain species feed heavily on food supplied by lateral advection in deep sea currents. Yet, the prevailing picture at all our sites is that there are only minor changes in the species assemblage and, instead, the changes in total foraminiferal accumulation rates are reflected across all the major groups of benthic foraminifera (notwithstanding a small change in Reyk R; see Supplementary Figure 12). We suggest this may indicate that there has been no overall change in the style and type of food delivery (and other environmental conditions) at our sites, rather, there has just been a change (increase or decrease, depending on site) in the overall amount of food. This finding may partly be a limitation of typically only counting from the >150 µm fraction, and greater % relative abundance changes, for example in smaller phytodetritus species, may have been obtained with examination of the >63 µm fraction (avoided due to concerns that smaller foraminifera could be transported by strong bottom currents). Yet at site Reyk R, which lies in the path of downstream ISOW controlled by LSW (Boessenkool et al., 2007) and where we counted the >63 µm fraction, we also saw no major shifts in the relative abundance of species. Further investigation of new cores throughout the North Atlantic can test if there are deep ocean sites where major shifts in species assemblage have occurred. Available published datasets and newly generated grain size records from this study (Supplementary Figures 17-20) also suggest that changes in the strength of the deep sea currents associated with our different sites, which could have altered the supply of organic matter to our sites by variability in lateral advection, were not a strong control on benthic foraminifera accumulation rates. This is further supported by our observation of no major shifts in the species assemblage, which likely would be expected if there was a significant shift in the proportion of food delivered by lateral advection rather vertical sinking.

CONCLUSION
This study has analyzed benthic foraminifera abundances in high accumulation rate sediment cores from the North Atlantic to investigate if the previously documented anomalous ocean circulation changes of the industrial era have impacted benthic ecosystems. Our results reveal that there are spatially variable changes in benthic foraminiferal abundances during the industrial era, with abundance levels in the 20th century at many sites emerging as exceptional compared to anything seen throughout the prior several thousand years of the Holocene. The most substantial changes have occurred in regions where there have been concomitant anomalous shifts in surface ocean conditions such as temperature and surface productivity. We infer that shifts in surface ocean frontal activity altered the food supply to the benthos, such that areas with decreased (increased) surface frontal activity saw a decrease (increase) in benthic abundance and sites where little surface ocean variability occurred also recorded little-to-no change in benthic foraminiferal abundance.

FUTURE WORK AND IMPLICATIONS
This study represents a preliminary investigation into industrial era changes in North Atlantic deep-sea benthic ecosystems, specifically benthic foraminifera. The obvious next step is to better map out the recent and ongoing shifts in benthic foraminifera reported here, and test the hypothesis that anomalous benthic shifts occur at sites where there have been large magnitude surface ocean changes, particularly ocean fronts. Prior work has been limited by the availability of sufficiently high-resolution sediment records, and future work will rely upon obtaining and examining cores from many more high sedimentation rate sites, with obvious benefits from gathering new cores that can capture some of the large magnitude ocean changes that have occurred in the most recent decades. However, this information will inherently be biased since it depends on selecting sites with high accumulation rates, likely placing restrictions on the types of benthic community that can be examined. Therefore, as far as possible, a wide range of benthic environments should be investigated for marine sediment core work, and this information should be combined with evidence of deep sea environment and ecosystem change from other archives, such as deep sea corals. As greater time overlap builds between our proxy archives and time series of observed ecological changes, robust testing of the relationships between these two types of data will be possible. Finally, in addition to ecological investigation, additional work is also needed to better constrain and understand the physical oceanographic changes that occurred during the industrial era, for example examining the regional expressions of changes in the AMOC and SPG circulation, and how they may vary in the future.
Our results have broad implications for marine ecosystem management and planning. The impact of changes in large scale circulation features like the AMOC have been modeled previously, for example, a weakening AMOC has been linked to a broad decline in marine productivity (Schmittner, 2005;Whitt and Jansen, 2020). Simple predictive theory suggests that under global warming, reduction in northward nutrient transport is the primary suppressor of Subarctic Atlantic Ocean productivity, while local reductions in vertical mixing have an important, but secondary, effect (Whitt and Jansen, 2020). The spatial heterogeneity that we observe in 20th century shifts in deep sea benthos, linked to surface ocean variability, cautions that accurate projection of future deep sea benthic ecosystems will require models that can capture the localto-regional scale features of surface ocean circulation -wellillustrated by the complexity of the northwest Atlantic shelf and slope region (Saba et al., 2016;Claret et al., 2018). Marine ecosystem projections based on CMIP5 climate models, which do not adequately resolve many of these features (e.g., Saba et al., 2016), are probably missing likely hotspots of future benthic ecological change. Based on the results of our study, improved modeling capability will be of particular value where deep sea benthic sites of special importance for ecological reasons or socio-economic reasons coincide with regions expected to undergo large shifts in surface ocean properties and frontal migration.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
DT conceived and designed the project. PS, JW, EP, ND, DF, RG, TL, FS, FP, and DT collected the data. CO'B, PS, and JW conducted the data analysis. PS, CO'B, and DT produced the figures and wrote the manuscript. All authors contributed to initial data interpretation.

FUNDING
Funding was provided to DT by NERC Project ReconAMOC (NE/S009736/1); the Leverhulme Trust, and ATLAS and iAtlantic projects. This project has received funding from the European Union's Horizon 2020 Research and Innovation Programme under grant agreement Nos 678760 (ATLAS) and 818123 (iAtlantic).