Response of biological productivity to North Atlantic marine front migration during the Holocene

. Marine fronts delineate the boundary between distinct water masses and, through the advection of nutrients, are important facilitators of regional productivity and bio-diversity. As the modern climate continues to change, the migration of frontal zones is evident, but a lack of information about their status prior to instrumental records hinders future projections. Here, we combine data from lipid biomarkers (archaeal isoprenoid glycerol dibiphytanyl glycerol tetraethers and algal highly branched isoprenoids) with planktic and benthic foraminifera assemblages to detail the biological response of the marine Arctic and polar front migrations on the North Iceland Shelf (NIS) over the last 8 kyr. This multi-proxy approach enables us to quantify the thermal structure relating to Arctic and polar front migration and test how this inﬂuences the corresponding changes in local pelagic productivity. Our data show that following an interval of Atlantic water inﬂuence, the Arctic front and its associated high pelagic productivity


Introduction 25
Marine fronts separate different water masses and are a globally ubiquitous feature in the oceans (Belkin et al., 2009). By nature, marine fronts are characterized by strong horizontal gradients of typically correlated water properties, such as temperature, salinity and nutrients. Some fronts influence the position of geostrophic currents that act as conduits for heat, salt and nutrient transport (i.e. density front), whereas others contribute to localized hot spots of primary productivity and biodiversity (i.e. convergent front) (Belkin, 2004;Belkin et al., 2009). Moreover, the combination of downwelling with 30 4 this study's inter-site comparisons, we focus on the biomarker records that have been analysed in both cores (Fig. 2, Bendle and Rosell-Melé, 2007;Cabedo-Sanz et al., 2016;Kristjánsdóttir et al., 2017). Given that sedimentation rates are nearly linear 95 at both locations over the last 8 ka, we do not correct proxy records for sediment flux.
Consequently, the presence of IP25 indicates the occurrence of seasonal sea ice along the sea ice margin whereas the presence of HBI III reflects highly productive Marginal Ice Zones (MIZ) characterized by open water phytoplankton blooms (Belt et 110 al., 2015(Belt et 110 al., , 2019. Further, the proportion of HBI III to HBI IV can provide evidence for the spring phytoplankton bloom .

Archaeal isoprenoid glycerol dibiphytanyl glycerol tetraethers
We present new subsurface temperature (subT) proxy records for each marine sediment core based on the distribution of 115 isoprenoid glycerol dibiphytanyl glycerol tetraether (GDGT). GDGTs are cell membrane lipids biosynthesized by marine ammonia oxidizing Thaumarchaeota (Könneke et al., 2005;Schouten et al., 2013), which on the NIS vary in terms of their cyclization (number of cyclopentane moieties) according to winter and/or annual subT (0-200 m, Harning et al., 2019), and potentially to a lesser extent, ammonia oxidation rates (e.g. Hurley et al., 2016). Total lipids extracts (TLEs) were obtained from freeze-dried marine sediment subsamples (~2-3 g, n=53 for MD99-2269, n=18 for JR51-GC35) by ultrasonication using 120 dichloromethane:methanol (2:1, v/v). TLEs were then suspended in hexane:isopropanol (99:1, v/v), sonicated, vortexed, and filtered using a 0.45 µm PTFE syringe filter. Prior to instrumental analysis, samples were spiked with 10 ng of the C46 GDGT internal standard (Huguet et al., 2006). Isoprenoid GDGTs were identified and quantified via high performance liquid chromatography -mass spectrometry (HPLC-MS) a Thermo Scientific Ultimate 3000 HPLC system interfaced to a Q Exactive Focus Orbitrap-Quadrupole MS after Harning et al. (2019). We adopt the TEX86 L index (Kim et al., 2010) to quantify the 125 distribution of GDGTs  S1). HBI IV is often co-produced by the same open-water diatoms that produce HBI III (e.g. Rhizosolenia spp., Rowland et al., 2001) and the combination of the two has recently been shown to be a useful predictor of spring phytoplankton blooms in the Barents Sea . Analysis of purified lipid fractions containing HBIs was 135 carried out using gas chromatography-mass spectrometry (GC-MS) in total ion current (TIC) and selected ion monitoring (SIM) modes after Belt et al. (2019). HBI IV was identified based on its characteristic GC retention index (RIHP5MS = 2091) and mass spectrum (Belt et al., 2000;Belt, 2018). HBI IV quantification was achieved by comparison of mass spectral responses of its molecular ion (m/z 346) in SIM mode with those of the internal standard (9-OHD, m/z 350) and normalized according to an instrumental response factor. HBI III and HBI IV datasets are then expressed through the following equation 140 where T25 ≥ 1 provides a qualitative proxy measure for the occurrence of spring phytoplankton blooms .
Moreover, Rowland et al. (2001) demonstrated a systematic increase in the amount of HBI III relative to HBI IV (i.e. higher T25) with increasing growth rate of Rhizosolenia setigera cultured at higher temperatures.

Planktic and benthic foraminifera
The modern distributions of planktic and benthic foraminifera from marine surface sediment around Iceland provide a key baseline for interpreting local environmental changes (Johannessen et al., 1994;Rytter et al., 2002;Jennings et al., 2004).
Based on this framework, downcore records from MD99-2269 have focused on certain indicator species at low-(benthic and planktic species, Justwan et al., 2008) and high-resolution (planktic species, Cabedo-Sanz et al., 2016). We expand these 150 records by presenting high-resolution and complete assemblage counts from 273 planktic and 295 benthic foraminifera samples from MD99-2269. The original high-resolution planktic records for T. quinqueloba and N. pachyderma (s) from MD99-2269 have been updated and represent the finalized dataset.
All sub-samples were prepared by wet sieving at 63 µm, 106 µm and 150 µm. Our dataset combines samples that had been air-dried at 35 o C after sieving and stored dry together with ones that were not dried . The 155 former sample set was wetted prior to wet splitting, and all samples were reanalysed for assemblages in a buffered solution.
Planktic (>150 µm) and benthic (>106 µm) foraminiferal assemblages are expressed as percentages of the total planktic and total benthic population. In addition to individual foraminiferal species profiles, we also use our assemblage datasets to estimate 6 bottom water temperature (BWT, benthic foraminiferal dataset) and sea surface temperature (SST, planktic foraminiferal dataset). Temperature estimates were quantified using North Atlantic transfer functions (BWT, Sejrup et al., 2004;SST, 160 Pflaumann et al., 2003) and Weighted Averaging Partial Least Squares Regression (WAPLS) and Maximum Likelihood (ML) techniques. Although the transfer functions assume a relationship with spring/summer temperature likely due to the seasonal flux of phytodetritus that heterotrophic foraminifera feed on, foraminifera productivity likely occurs at other times of the year, which result in annually integrated temperature records.

Statistical analyses
Step 1: To visualize long-term trends in proxy time series, we performed locally weighted smoothing (LOESS) to help minimize the influence from outliers and short-term variability while retaining persistent shifts. The smoothing criterion for each time series was selected automatically after optimization using generalized cross validation. Values along the resulting LOESS curves were extracted using a 60 yr timestep, representing a balance between both the arithmetic mean (~90 yr) and 170 median (~30 yr) of first derivatives for each time series.
Step 2: To detect persistent step shifts within individual time series, we performed Sequential T-test Analysis of Regime Shifts (STARS) on the LOESS-smoothed data (Rodionov, 2004(Rodionov, , 2006. The algorithm was tuned to detect regime shifts on millennial timescales by setting the cut-off length to 60, and we report the timing of shifts identified at the 95% confidence level (p ≤ 0.05). The timing of regime shifts is best interpreted as approximate, as their exact timing and magnitude 175 are affected by the chosen confidence level and cut-off length (e.g. Rodionov, 2006). Comparison of regime shift timing between sediment cores should be approached with caution due to the variable resolution and quality of chronological control.
Step 3: To evaluate the relative similarity between select MD99-2269 productivity and temperature proxy time series, we performed complete-linkage Agglomerative Hierarchical Clustering (AHC) using Dynamic Time Warping (DTW) distance as a curve shape-based dissimilarity. In order to preserve the influence of short-term variability while rectifying missing data, 180 we avoided using LOESS smoothing and linearly interpolated the proxy records across all available horizons (600+ data points).

Trends and regime shifts identified in existing proxy records 185
Alkenone-derived U K' 37 spring SST records from MD99-2269  and JR51-GC35 (Bendle and Rosell-Melé, 2007) reveal different overall patterns of temperature variability during the last 8 ka. The record from MD99-2269 is relatively lower in absolute spring SSTs as well as in the magnitude of shorter-term variability. It features four regime shifts rather than five as found in JR51-GC35 ( Fig. 2a- b). Both U K' 37 records indicate that highest spring SSTs were achieved during the earliest portion of the records with the lowest spring SSTs occurring variably between ~4 and 2 ka BP ( Fig. 2a-b).

T25 210
HBI IV was present above the detection limit in all sediment samples (Fig. S3). The T25 records for both cores show similar trends, with highest values achieved during the early Middle Holocene and progressively lower values towards present . For marine sediment core MD99-2269, the T25 record exhibits four negative regime shifts, whereas for marine sediment core JR51-GC35, the T25 record exhibits five negative regimes shifts ( Fig. 2i-j). The overall decreasing LOESS trends as well as the consistently negative regimes shifts in both records indicate a generally continuous decline in pelagic productivity over 215 the last 8 ka (e.g. . Further, given that HBI III is produced at a relatively higher proportion under conditions of more rapid growth in some diatoms (e.g. R. setigera, Rowland et al., 2001), these records also suggest a long-term decrease in diatom growth rates.

Planktic and benthic foraminifera 220
Seven planktic foraminiferal species were identified in the >150 µm size fraction and used for SST estimation (Fig. S4). Our reconstruction of planktic foraminifera assemblage and SST estimates include low-and high-resolution data before and after 8 ka BP, respectively ( Fig. S4-S6). We focus on two planktic indicator species for our paleoceanographic discussion; T. quinqueloba and N. pachyderma (s). T. quinqueloba is an AF species (Johannessen et al., 1994;Volkmann, 2000;Pados and Spielhagen, 2014) that shows highest abundance during the early Middle Holocene and lowest abundance during the Late 225 Holocene, with four negative regime shifts (Fig. 3c). N. pachyderma (s) is a Polar Water species (Johannessen et al., 1994;8 Jennings et al., 2004) that shows lowest abundances during the early Middle Holocene, and three positive shifts during the Late Holocene when its maximum abundance is achieved (Fig. 3d). In terms of foraminifera-reconstructed SST (S.E.=1.3 o C), there is an overall trend of cooling throughout the last 8 ka from ~10 °C to ~3 °C, although SST from 8 to ~4 ka BP and from ~3 ka BP to present are relatively more stable than between 4 and 3 ka BP (Fig. 4b). Four negative regime shifts are identified 230 throughout the SST record, the most pronounced of which occurs at ~4 ka BP (Fig. 4b).
In the benthic dataset, we focus on the two most abundant species, Cassidulina neoteretis and Cassidulina reniforme, for paleoceanographic interpretations. Around Iceland, C. neoteretis is abundant on the western NIS where it is associated with NIIC Atlantic Water and dramatically decreases eastward where Atlantic Intermediate Water forms the bottom water (Rytter et al., 2002;Jennings et al., 2004). C. reniforme is an Arctic species prevalent on the NIS but has its highest abundance in the 235 low-salinity fjords of NW Iceland (Rytter et al., 2002;Jennings et al., 2004). In our records, C. neoteretis shows maximum abundances during the Middle Holocene and exhibits three regime shifts throughout the last 8 ka (Fig. 3f). C. reniforme peaks during the early Middle Holocene and has five regime shifts throughout the last 8 ka (Fig. 3g). We note that another dominant species is Nonionella iridea, which formed up to 20 % of the benthic fauna (Fig. S5) and included in our statistical analyses to better understand its environmental preferences (see Supporting Information Text S1). Thirty-three of the 65 species in the 240 Sejrup et al. (2004) transfer function are present in MD99-2269, although neither N. iridea nor the agglutinated species were included (Sejrup et al., 2004). The BWT estimates are less variable than the foraminifera-based SST estimates, with temperatures generally falling between ~2.5 and 5 °C and maximum temperatures occurring at ~4 ka BP (S.E.=1.0 o C). Two positive regime shifts are identified in the Middle Holocene followed by two negative regime shifts during the Late Holocene, consistent with the general long-term pattern of warming from 8 to ~4 ka BP and cooling from ~4 ka BP to present (Fig. 4d). 245

Agglomerative Hierarchical Clustering analysis
The AHC dendrogram of DTW distances shows that the 11 selected proxy records from MD99-2269 successfully separate into four clusters (Fig. 5). The first grouping of T25, foraminifera SST, and TEX86 L subT clusters proxies that are influenced by annual near surface to surface temperatures. The spring U K' 37 SST belongs to its own cluster likely due to seasonality 250 differences that result in a relatively smaller range of temperature changes as well as warming since ~2 ka BP (Fig. 4a-c). The third grouping of HBI III, T. quinqueloba and, N. iridea, and to a slightly lesser extent, HBI IV, suggests that this cluster reflects a separate but dominant environmental variable on the NIS, likely pelagic productivity (see Supporting Information Text S1 for discussion on N. iridea). The close relationship of the productivity cluster (blue) to the temperature clusters (purple and green) highlight the fact that SST and productivity are innately connected on the NIS as warmer Atlantic Water with 255 sufficient vertical mixing often carry higher nutrient loads that stimulate primary production (e.g. Thordardóttir, 1984;Zhai et al., 2012). Finally, the fourth grouping of IP25, N. pachyderma (s) and foraminifera BWT likely clusters cold proxy indicators, such as the sea ice biomarker IP25 and planktic foraminifera species associated with Polar Water and bottom waters. The fact that foraminifera BWT clusters in the fourth grouping may relate to the relatively depressed temperatures in the early portion of the record, despite a negative trend from ~4 ka BP to present (Fig. 4d) that would be consistent with trends observed in the 260 9 aforementioned near surface temperature and productivity proxies (Figs. 3 and 4). In any case, we note that the similarity between BWT and Polar Water indicators is weaker than that between IP25 and N. pachyderma (s).

Frontal proxies on the NIS 265
Modern distributions of T. quinqueloba along AF systems suggests that the species feeds in the high nutrient and productive waters within the warmer margins (Volkmann, 2000;Pados and Spielhagen, 2014). In prior work from MD99-2269, Cabedo-Sanz et al. (2016) used Principal Component Analysis to show that HBI III and T. quinqueloba were related, along with biogenic CaCO3, as proxies for surface productivity. Our AHC analysis, which uses the updated T. quinqueloba record, reinforces the relationship previously observed between HBI III and T. quinqueloba time series on the NIS (Fig. 5). Further, 270 given the known environmental preferences of T. quinqueloba, we suggest that the NIS diatom producers of HBI III thrive in similar hydrographic conditions. Although we lack information on diatom sources in our study, the distribution of some HBI III-producing diatoms (i.e. Rhizosolenia hebetata f. semispina) do trace the AF's modern position (Oksman et al., 2019), and are important members of Arctic Water assemblages in the North Atlantic (Koç Karpuz & Schrader, 1990;Oksman et al., 2019). Collectively, this evidence indicates that HBI III and T. quinqueloba track the migration of the AF along the NIS. 275

Environmental controls on benthic foraminifera Nonionella iridea
N. iridea is an often overlooked or missed constituent in benthic faunal assemblages (e.g. Sejrup et al., 2004) due its small size  and because its abundance is underestimated using dry analyses. Experimental work indicates that N.
iridea may feed on seafloor phytodetritus and/or the associated suboxic-hypoxic bacterial populations that can develop, but 280 only dominates the assemblage in response to pulsed phytodetritus delivery (Gooday and Hughes, 2002;Duffield et al., 2015).
The wet picking techniques employed in our study allowed us to quantify downcore variations in N. iridea and better understand its environmental preferences and role in Icelandic benthic communities. N. iridea shows highest taxa abundance during the Middle Holocene, which features one negative regime shift. Lowest abundance occurred during the Late Holocene, which features one negative regime shift followed by one positive (Fig. 6). Overall, the structure of the N. iridea, HBI III and 285 T. quinqueloba records are similar (i.e. decreasing toward present) and the timing of regime shifts between them is consistent ( Fig. 6). This notion of similarity is further supported by group clustering in our AHC analysis (Fig. 5). Although future studies including modern foraminifera distributions that include this species are still needed, this statistical evidence suggests that N.
iridea's food supply may have been influenced by the presence of frontal systems and/or warmer waters. In other words, the enhanced pelagic productivity of the AF likely resulted in increased export of phytodetritus to the seafloor where it could be 290 consumed by N. iridea.
In this interval, planktic foraminifera assemblages and the derived annual SST estimates indicate that NIIC Atlantic Water was 295 entering the MD99-2269 site, whereas Arctic Intermediate Water occupied the lower depths as indicated by the benthic fauna and low BWT estimates (Fig. 4b-c and 7a). Coccolithophore assemblage data indicative of NIIC Atlantic Water from 8-7 ka BP from the same core (Giraudeau et al., 2004) support the warm planktic foraminifera SST reconstruction. Both of our HBI III records (Figs. 2g-h and 3a), as well as % T. quinqueloba (Fig. 3c), document low surface diatom productivity until ~6.1 ka BP, which we interpret to reflect a northward position of the AF relative to the NIS (Fig. 7a). Although overall diatom biomass 300 was low, high T25 values indicate increased diatom growth rates (Fig. 3b), likely driven in part by higher temperatures (Fig. 5).
In addition, total benthic and planktic foraminifera counts were highest in this interval (Fig. S6), which may also indicate a long ice-free season for production of some Arctic species. Although this interval of low phytoplankton productivity was previously interpreted to reflect strong convection in the Iceland Sea and the production of Arctic Intermediate Water (Giraudeau et al., 2004), we note that these interpretations are inconsistent with our planktic foraminifera assemblage dataset. 305 However, given the relatively high orbitally-induced seasonality at this time (Berger and Loutre, 1991), the periodic production of cool and fresh North Icelandic Winter Water (Stefánsson, 1962) may have increased along the northern coastline. The potential contribution of these local waters along with the distal position of the AF would have dampened surface productivity for many frontal marine photosynthesizers (Fig. 7a).
Our datasets from MD99-2269 provide a comprehensive view of corresponding ocean temperatures on the western 310 NIS during this early Middle Holocene interval. U K' 37 and foraminifera SST estimates show respectively high yet decreasing spring and annual surface temperatures ( Fig. 4a-b). In contrast, our TEX86 L -derived subT and foraminifera BWT records show depressed yet increasing temperatures (Fig. 4c-d). Hence, MD99-2269's temperature profile depicts a strong annual vertical temperature gradient during the Early-Middle Holocene (~8 o C, Fig. 4e), resulting from stratification of Arctic Intermediate Water at lower depths and NIIC Atlantic Water at the surface. At the JR51-GC35 core site, we also observe similarly high and 315 stable SST and relatively lower corresponding subT ( Fig. 2b and 2d). However, a lack of BWT estimates in this core prevents us from analysing the vertical temperature gradient on the eastern NIS. Comparison of TEX86 L subT temperature records between the western and eastern NIS show comparable absolute estimates (Fig. 4f), which likely reflects the invasion of NIIC Atlantic Water across the surface to subsurface of the entire NIS (Fig. 7a). 320

6.3 to 3.4 ka BP (local stabilization of the Arctic Front)
At the beginning of this interval, planktic and benthic foraminifera indicate that NIIC Atlantic Water dominated the surface whereas Arctic Intermediate Water occupied the lower depths on the western NIS (Fig. 3). At ~6.1 ka BP, abrupt increases in HBI III and T. quinqueloba abundance from MD99-2269 indicate significantly enhanced pelagic primary production on the western NIS, which although variable, is sustained to ~3.8 ka BP (Fig. 3). Although the HBI III record from JR51-GC35 shows 325 similar trends, the overall increase in inferred pelagic productivity is smaller (Fig. 2h). We interpret this spatial variability to reflect that by 6.1 ka BP, the AF was positioned near MD99-2269 but never or rarely expanded eastward to JR51-GC35, and that NIIC Atlantic Water was less of a contributor to the JR51-GC35 site (Fig. 7b). The increase in certain dinoflagellate and coccolithophore algae species (e.g. Coccolithus pelagicus) have also been used to track the presence of the AF to MD99-2269 from 6.2 to 3.5 ka BP (Giraudeau et al., 2004;Solignac et., 2006). Not only do the dinocyst and C. pelagicus records further 330 support the timing of AF presence around MD99-2269, they also reflect short-term oscillations between NIIC Atlantic Water and EIC Arctic Waters (Giraudeau et al., 2004;Solignac et al., 2006). Hence, although the AF was located around MD99-2269 in this interval, its precise location on the NIS was not static.
In terms of temperature, U K' 37 and planktic foraminifera proxies reflect a decrease in spring and annual SST, respectively, beginning by ~5.3 ka BP followed by a second regime shift at ~4 ka BP (Fig. 4a-b). As nutrient availability was 335 presumably higher, diatom growth rate and cell division inferred from decreasing T25 was likely hampered by these decreasing SSTs (Fig. 3b). In contrast to the SST records, TEX86 L subT on the western NIS shows anomalous warmth from ~6 to 5 ka BP (Fig. 4c). Interestingly, this peak subT warmth closely aligns with the initial increase of HBI III and T. quinqueloba and southeastward migration of the AF over the western NIS at MD99-2269 (Fig. 2). As shown by the foraminifera fauna and temperature proxies, this frontal migration began to de-stratify the water column (Fig. 4e), which may have allowed warmer 340 water to advect to lower depths where Thaumarchaeota live (Fig. 7b). Alternatively, if Thaumarchaeota became stressed over increased competition for NH4 + with frontal primary producers, their ammonia oxidation rates could decrease, which in culture studies has been shown to increase GDGT-based temperatures (Hurley et al., 2016). If this is the case, the lack of increased TEX86 L subT changes at JR51-GC35 may suggest a reduced competition for nutrients and reinforce reduced advection of NIIC Atlantic Water and a distal position of AF on the east NIS. 345 Following the regime shift at ~3.8 ka BP in the HBI III record from MD99-2269, the PF expanded southeast over the NIS (Fig. 7c). At the surface, this is reflected by decreases in spring and annual SST, subT, BWT (Fig. 4a-d), decreases in % T. quinqueloba at the expense of % N. pachyderma (s), lower diatom growth rates (decreasing T25), increased IP25-inferred sea ice (Fig. 3), as well as continued thermal destratification (Fig. 4e). However, beneath the surface, benthic fauna show the progressive submergence of the NIIC Atlantic Water under EIC Arctic Water ( Fig. 3 and 7b). On the eastern NIS at JR51-350 GC35, subT and IP25-derived sea ice records also continue to decrease and increase, respectively, although to a lesser extent than at MD99-2269 (Fig. 2). Combined with records of water mass distribution inferred from radiocarbon reservoir ages  and low-resolution planktic foraminifera assemblage and their δ 18 O (Eiríksson et al., 2000;, these proxies reflect progressive cooling of EIC Arctic Water, rather than the shift in surface water source observed at MD99-2269 (i.e. NIIC to EIC) (Fig. 4f). Additional records of ocean temperature and iceberg rafting along the NIS indicate 355 contemporaneous millennial-scale cooling during this interval (e.g. Moros et al., 2006;Jiang et al., 2015) likely linked to changes in oceanic circulation, such as AMOC slowdown (Thornalley et al., 2009) and/or changes in atmospheric circulation patterns (Orme et al., 2018).
One record that remains challenging to interpret during the Middle Holocene is the U K' 37 SST at JR51-GC35. Bendle and Rosell-Melé (2007) connected the high amplitude variability of U K' 37 SST (Fig. 2b) to the strength of North Atlantic Deep 360 Water (NADW) formation and the Atlantic Meridional Overturning Circulation (AMOC). Although the amplitude of JR51-GC35's U K' 37 SST variability is inconsistent with our additional SST records, the timing of these SST changes is similar to 12 those inferred from paired Mg/Ca-δ 18 O measurements of planktic foraminifera south of Iceland that reflect the relative strength of the AMOC and Subpolar Gyre (Thornalley et al., 2009). However, if changes in the AMOC explain SST variability along the eastern NIS, it would be expected that similar variations would also be present in the MD99-2269 record. One possibility 365 that may explain the contrasting observations is that the warm NIIC positioned over MD99-2269 was fluctuating to and from the JR51-GC35 core site in accordance with variability in AMOC strength (Fig. 7b). Given that local coccolithophore communities have changed during the Holocene and that the species during this interval favor these dynamic surface conditions (Giraudeau et al., 2004), it is also possible that the high-amplitude changes in JR51-GC35's U K' 37 SST may relate to algal community changes associated with the varying EIC vs NIIC source water. 370

3.4 ka BP to present (southward migration of Arctic Front)
Consistently low HBI III concentrations at both core sites throughout the remainder of the records suggest that the PF remained a local feature along the NIS during the Late Holocene ( Fig. 2 and 7c). Further, planktic foraminifera (low % T. quinqueloba and high % N. pachyderma (s)) and rising concentrations of IP25 in MD99-2269 (Fig. 3) suggest that the EIC was now advecting 375 more Polar Water relative to Arctic Water. Accordingly, there are regime shifts toward lower spring and annual SST and subT from ~2 to 1.3 ka BP (Fig. 4), and from toward N. pachyderma (s) dominance at ~1.6 ka BP (Fig. 3). Although the annual vertical temperature gradient on the western NIS had decreased to <1 o C from its early Middle Holocene maximum (Fig. 4e), Late Holocene cooling is observed at all depths (Fig. 4). The low SSTs in combinations with low pelagic nutrient availability is likely responsible for the lowest Holocene diatom biomass and growth rates observed in the HBI III and T25 records, 380 respectively (Fig. 2). These observations fit well within the context of previous proxy research that has shown this period to be the coldest of the Holocene on the NIS (Eiríksson et al., 2000;Andersen et al., 2004;Jiang et al., 2015). Moreover, the periodic advection of sea ice around east Iceland and the increase in % T. quinqueloba in southwest Iceland marine sediment cores show that the AF had migrated to the southern coastline during the peak of this cooling from ~1760 to 1920 CE (i.e. Little Ice Age, Jennings et al., 2001). 385 Several regional proxy records, such as U K' 37 in MD99-2269, have recently noted an upturn in the derived spring/summer SSTs over the last two millennia (e.g. Moossen et al., 2015;Kristjánsdóttir et al., 2017), possibly related to regional changes of atmospheric circulation (Orme et al., 2018). In light of the continued decrease in northern hemisphere summer insolation (Berger and Loutre, 1991) and growth of Icelandic ice caps that require local summer cooling (Larsen et al., 2011;Harning et al., 2016Harning et al., , 2018Harning et al., , 2020Anderson et al., 2018;Geirsdóttir et al., 2019 suggested 390 that this "warming" may have been driven by either the presence of a thin freshwater lid that restricted mixing with deeper colder water and the supply of nutrients and/or a seasonal shift of alkenone production to later summer months following more persistent spring sea ice. Our foraminifera-derived annual SST from MD99-2269 (Fig. 4b) and foraminifera assemblage data ( Fig. S4 and S5) suggest that despite possible changes in seasonality, cold and low salinity Polar Water formed the surface waters even during the spring/summer. These results support both of the mechanisms proposed by Cabedo-Sanz et al. 395 13 (2016) and suggest that drawing conclusions about regional spring SST warming over the last 2 ka should be approached with caution until the nuances of the individual proxies are better understood.

Controls on the Holocene migration of the Arctic and Polar Fronts
The dominant climate forcing during the Holocene has been the first-order, orbitally-driven decrease in northern hemisphere 400 summer and annual insolation (Berger and Loutre, 1991). Our NIS SST records depict a similar first-order decrease ( Fig. 4a- b), consistent with an insolation-driven reduction in northward heat transport (e.g. Andersen et al., 2004). Although variability in the strength of northward heat transport has been modulated by centennial to millennal-scale changes in NADW formation and the stability and strength of the AMOC (e.g. Thornalley et al., 2009), as well as possible changes in atmospheric circulation (e.g. Orme et al., 2018), the progressive migration of the Arctic and Polar Fronts around Iceland suggest that insolation has 405 also remained its primary control. Given that the movement of the North Atlantic frontal systems is manifested by regional temperature and productivity variability in the past, our new records may provide useful analogues for ongoing fossil fueldriven warming.
Following the end of the Little Ice Age in ~1920 CE and a century of first-order warming (http://vedur.is), the AF and PF have returned to positions along Greenland and proximal to the NIS (Fig. 7d), similar to those during the warm early 410 Middle Holocene (Fig. 7b). Accordingly, enhanced primary productivity (increased T. quinqueloba) is noted in recent NIS sediment records, that in addition to the position of the AF and PF straddling the NIS, may also be driven by increased advection of Atlantic Water (Jónsson and Valdimarsson, 2012) and/or freshwater discharge from Greenland (Perner et al., 2019). As global temperature continue to rise and melting of the Greenland Ice Sheet accelerates, freshening of North Atlantic surface waters is expected to continue. Although recent melt may help stimulate this productivity (Perner et al., 2019), there likely 415 exists a threshold where too much freshening will restrict nutrient availability required for photosynthesis, as we observe during the Early Holocene (Fig. 7a). Further, the observed slowdown (Rahmstorf et al., 2015) and conceivable future shutdown of AMOC from Greenland melt (Bakker et al., 2016) may also result in a southward shift of North Atlantic frontal systems.
Although more emprical and modeling research is needed to forecast future climate trajectories, our records highlight the sensitivity of the AF and PF to temperature changes, which hold important implications for the local pelagic productivity and 420 economies that it supports today.

Conclusions
• We present new TEX86 L -derived subT and HBI productivity records from the western and eastern NIS, as well as new high-resolution planktic and benthic foraminifera assemblage (and corresponding temperature) records from the 425 western NIS; all covering the last 8 ka.
• The Arctic Front, a zone of intensified pelagic productivity (as indicated by HBI III and % T. quinqueloba), moderate phytoplankton growth rates and warmer waters, stabilized on the NIS by ~6.1 ka BP, with greater influence on the west than the east. By ~3.8 ka BP, the Arctic Front migrated south of the NIS, allowing for cold, Arctic/Polar Water 14 associated with the Polar Front (as indicated by IP25 and N. pachyderma (s)) to dominate the NIS for the Late 430 Holocene.
• Vertical temperature gradients on the western NIS were largest during the early Middle Holocene and progressively decreased to the lowest temperature gradients during the Late Holocene. Longitudinal temperature gradients suggest that warmer NIIC Atlantic Water was more influential on the west compared to the eastern NIS where cooler EIC Arctic Water dominated. 435 • The Holocene migration of the AF and PF has been primarily controlled by first-order decreases in northern hemisphere insolation, but the productivity it supports is also sensitive to freshening of surface waters. The future balance between these two variables will shape how the local configuration of marine fronts in the North Atlantic will develop under continued regional warming.

Data availability
Upon publication data will be uploaded to the PANGAEA database. During the review of the manuscript, data will be made available upon reasonable request to the authors.

Declaration of Competing Interests 445
The authors declare that they have no conflict of interest.

Acknowledgements
We extend great thanks to the scientific parties and crews of the RRS James Clark Ross and R/V Marion Dufresne II cruises for initial sediment core collection, as well as the numerous authors of previous studies (and in particular Dr. John T. Andrews) 450 who have contributed to both cores' comprehensive datasets and planted the seeds for continued research. We also thank Dr.
Suzanne Maclachlan at the British Ocean Sediment Core Research Facility (BOSCORF), UK, for making sediments available from core JR51-GC35 and Dr. Patricia Cabedo-Sanz for providing us with HBI IV biomarker data. University of Colorado students Allyson Rugg, Ian Courtney, Kelly Curtis, Jessica Scherer and Kylie Smith counted the foraminiferal assemblages in MD99-2269, expanding on the initial assemblage data produced by Nancy Weiner. This study has been supported by the 455 Icelandic Research Council (RANNIS) Grant of Excellence #141573-051 to ÁG and co-authors, as well as internal University of Colorado Boulder funds to JS.

Author Contributions
DJH conceived the study following discussions with AJ and STB. AJ led foraminifera analyses, DK extracted TLEs and 460 performed statistical analyses, and DJH performed GDGT analyses. STB, ÁG and JS supported analytical expenses. DJH wrote the manuscript with contributions and discussion of data interpretations from all co-authors.