Introduction

The Gulf of Maine, a semi-enclosed sea on the east coast of North America (Fig. 1), has recently undergone significant warming, including increases in sea surface temperatures (SSTs) greater than much of the global ocean (0.23 °C yr−1 from 2004–2013)1, extreme marine SST heatwaves, including in 2012 and 20162,3, and subsurface temperature increases, particularly in the late summer, fall and winter, as high as 0.5 °C yr−1 since 20044. This extreme warming has already had a significant effect on Gulf of Maine ecosystems1,2,4. Warming has been linked to the stagnation or collapse in the recovery of several species and fisheries in the Gulf of Maine, including right whales (Eubalaena glacialis)4 and the cod (Gadus morhua) fishery1. Warming has also been implicated in the change in timing of lobster landings, which has significantly impacted the world-renowned lobster industry in this region2.

Fig. 1: Map of study site and broader North Atlantic region.
figure 1

Arrows indicate schematic depictions of major ocean currents while colors indicate average sea surface temperatures. Markers indicate the location of the study site as well as other paleoceanographic records mentioned in the manuscript (see map legend). The site of the Boothbay Harbor sea surface temperature (SST) instrumental record is indistinguishable from the site of this study (green marker) on this map. Bathymetry contours are marked in white at 200, 1200, 2200, and 3200 m depths from Smith and Sandwell88. SSTs are from the World Ocean Atlas 2018 climatological mean (from 1981 to 2010)89. The black box in the main figure shows the Gulf of Maine region highlighted in more detail in the inset. Black arrows in the inset depict surface currents, with gray arrows showing sub-surface currents. LC Labrador Current, DWBC Deep Western Boundary Current, GS Gulf Stream, NEC Northeast Channel, GB Georges Bank.

The extreme warming in the Gulf of Maine and surrounding region has multiple causes, likely including increases in atmospheric temperatures3 and pressures5, increases in Gulf Stream warm core rings in the region6, changes in the Gulf Stream western destabilization point7, and the influence of various internal modes of climate variability1,5,8. Several studies have implicated the position of the Gulf Stream relative to the Gulf of Maine as a driver of some of the extreme warming that this region has undergone recently1,4,5,9,10, with a more northerly Gulf Stream resulting in warmer water throughout the water column in the Gulf of Maine region9. In particular, several studies have noted the relationship between the strength of the Atlantic Meridional Overturning Circulation (AMOC) and the position of the Gulf Stream in the region11,12 and have suggested that a warming western North Atlantic region may be related to a weakening AMOC10,13,14. The deep, returning limb of the AMOC, the Deep Western Boundary Current, results in a northern recirculation gyre south of the Grand Banks11. A stronger AMOC and Deep Western Boundary Current results in a stronger recirculation gyre, which keeps the Gulf Stream farther offshore in that region12. The opposite is true for a weaker AMOC and Deep Western Boundary Current, allowing the Gulf Stream to shift closer to the Gulf of Maine.

The Gulf of Maine warming and the recent changes on the western North Atlantic shelf have been assessed using modeling results9 or very short (<40 years) instrumental records1. The longest, continuous instrumental record available is an SST record from Boothbay Harbor (BBH), Maine, which extends back to 1905. The Gulf of Maine has been generally warming since this instrumental record commenced, undergoing a temperature increase of 2.2 °C century−1 (based on annual averages of the BBH SST record from 1906 to 2018). Very little is known about Gulf of Maine water properties or the drivers of water property variability before the BBH SST record was installed. Wanamaker, et al.15 use oxygen isotopes in shells of the marine bivalve Arctica islandica (ocean quahog) to suggest that the Gulf of Maine has been cooling by 1–2 °C over the last millennium, consistent with other records in the western North Atlantic16. However, the time series presented in the Wanamaker, et al.15 study is not continuous and includes measurements from only four individual specimens. It is therefore not known when the warming documented by the instrumental records began or how unique the warming rate is in the context of the last several centuries. Additionally, uncertainties in the drivers of these documented temperature trends still remain.

Here we present a 300-year reconstruction of Gulf of Maine hydrographic variability using three different proxies of water properties, which are then spatially and temporally assessed in the context of the Last Millennium Ensemble (LME) simulations of the Community Earth System Model (CESM). The reconstructed water properties in this study are derived from chemical proxies from A. islandica shells (Supplementary Fig. 1): oxygen isotopes (δ18O) as a proxy for seawater hydrography (temperature and salinity) and nitrogen isotopes (δ15N) as a proxy for water mass source. Additionally, these water properties are compared to a previously published record of age-corrected radiocarbon (Δ14C)17 measured in the same A. islandica shells as those measured for δ18O and δ15N18. The objectives of this study are to: (1) reconstruct the water mass history in the Gulf of Maine for the last 300 years using multiple proxy records; (2) evaluate this hydrographic variability in the context of the last 1000 years, particularly seeking to investigate possible driving mechanisms using the LME; and (3) determine when recent observed warming in the Gulf of Maine started and what may have triggered this warming.

All three water properties reconstructed in this study differ between the different water masses entering the region: Warm Slope Water (WSW), formed adjacent to the Gulf Stream, and Labrador Slope Water (LSW), transported from the Labrador Sea by the Labrador Current19, are sub-surface water masses that enter the Gulf of Maine through the Northeast Channel (Fig. 1). Scotian Shelf Water (SSW) forms on the Scotian Shelf and is transported into the Gulf of Maine as surface water20,21. These reconstructed water properties can therefore be used to track changes in the relative proportion of these slope waters through time. WSW is warmer19, younger with respect to radiocarbon22,23,24, and has lower nitrogen isotope values of dissolved nitrate25,26 when compared to LSW. SSW, being a surface water mass, has a larger seasonal temperature range, and tends to be younger (with respect to radiocarbon) than LSW but older than WSW, although SSW has overlapping radiocarbon values with both WSW and LSW17. Lower‐Spies, et al.17 detail specific Δ14C ranges expected for the different water masses influencing the Gulf of Maine based on previous studies. In brief, Gulf Stream waters, the component of WSW with the highest Δ14C values, range between −62.4‰ and −22.7‰. SSW Δ14C values range between −71.9‰ to −50.8‰, while LSW has Δ14C values that range between −70.6‰ and −64.4‰. There are sparse data available on the nitrogen isotope composition of dissolved nitrate in SSW but, as a surface water mass, SSW is likely to contain dissolved nitrate that has a higher nitrogen isotope signature than the slope waters.

The geochemical records presented in this study reveal that the Gulf of Maine has been warming since the late 1800s, driven at least in part by changes in the proportion of different water masses entering the Gulf of Maine, likely related to changes in western North Atlantic ocean dynamics. Climate model simulations agree with the timing of the warming and suggest that this warming trend reverses at least 900 years of cooling in the region, cooling that appears to be largely driven by volcanic forcings. The simulations further reveal that the Gulf of Maine warming seen in the most recent century is more rapid than during almost all other 100-year periods in the last 1000 years and that such warming is driven, at least in part, by greenhouse gas forcings.

Results

Two separate proxy records of water mass variability (δ18O and δ15N)27 are reconstructed (Fig. 2) from shells collected in the western Gulf of Maine (Fig. 1) and absolutely dated in a previously published chronology using cross-dating techniques (Supplementary Fig. 2)8,28. The δ18O record is an annually resolved record and extends from 1694 to 2013 CE (n = 270), with gaps in data from 1703 to 1746, 1749, 1762 to 1763, and 1767 to 1768. The δ15N record is a decadally resolved record, which extends from 1751 to 2008 CE (n = 78). The amount of time integrated varies from 1 to 21 years (average = 5 years, 2σ = 10 years). Samples dated after 1927 were previously published in Whitney, et al.26. Additionally, these records are compared to a radiocarbon (Δ14C) record29 previously published in Lower‐Spies, et al.17. The Δ14C record is a decadally resolved record and extends from 1685 to 1935 CE. Samples included in the Δ14C record integrate between 3 and 30 years of data (average = 11 years, 2σ = 13 years, n = 34). Both the δ15N and Δ14C data are plotted at the midpoint of the years that each sample integrates.

Fig. 2: Geochemical, climate model, and instrumental time series.
figure 2

a Oxygen isotopes measured in Arctica islandica shells. Markers are the average value for that year. Black lines show the result of the segmented regression analysis. Light shading shows the average 1σ of the data, calculated from those years when more than one shell was sampled, darker shading shows the average 1σ analytical uncertainty. Note that the y-axis is inverted. b The age-corrected radiocarbon record, originally published in Lower‐Spies, et al.17. Horizontal error bars represent the amount of time incorporated in each sample. Shading shows the average 1σ analytical uncertainty. c The nitrogen isotope record. Data after 1926 was previously published in Whitney, et al.26. Markers are the average value for that year. Horizontal bars represent the amount of time incorporated in each sample. Light shading shows the average 1σ of the data, calculated for those years when more than one shell was sampled, darker shading shows the average 1σ analytical uncertainty. Note that the y-axis is inverted. d The CESM-LME Gulf of Maine annual temperature anomaly time series, calculated at 35 m depth. Anomalies are calculated using the average of the entire record shown here. The bold line is the ensemble member mean, the thin lines show each of the 13 ensemble members with full transient forcing. The linear regression lines show the result of the segmented regression analysis on this record. e Temperature from instrumental records (BBH) and gridded data products (ERSST31, HadISST32, OISST34, EN436).

To assess the utility of using δ18O as a sea surface temperature proxy in the Gulf of Maine, we compared this record to various Gulf of Maine instrumental temperature records and gridded data products. Table 1 shows the maximum Pearson correlation coefficients, after accounting for autocorrelation, between the δ18O record presented in this study and an instrumental dataset (the BBH SST record30) as well as various gridded data products (Extended Reconstructed Sea Surface Temperature – ERSST - v531, Hadley Center Sea Ice and Sea Surface Temperature – HadISST- v1.132,33, Optimum Interpolation Sea Surface Temperature – OISST-v2.134,35, EN4.2.1 data product36,37 for surface and 35 m) for the time period of overlap between all datasets (1983–2013). Supplementary Table 1 shows the Pearson correlation coefficients for the full range of overlap between the oxygen isotope record and each individual data product. See the Methods section for a detailed discussion of these data products, including the spatial domain over which the temperature average has been calculated for the purpose of these analyses. Because it is not known when precisely during the year A. islandica in this region precipitate shell material, we correlated the oxygen isotope record with the instrumental data products averaged over 12-month periods with various leads and lags. Table 1 reports the maximum correlation coefficient resulting from these analyses. Although the season of maximum correlation varied somewhat, most data products had the highest correlation with the oxygen isotope record for a 12-month period starting sometime in the fall months. Of note, for correlations calculated over the 1983–2013 period, the highest correlations occurred when the gridded data products lagged the oxygen isotope record by a year. This suggests that temperatures at the collection site lead average Gulf of Maine temperatures. The fact that almost all maximum correlation coefficients are significant at the 0.05 significance level suggests that the oxygen isotope record presented in this study is a useful proxy of annual average Gulf of Maine SSTs. Variability in the correlations between the instrumental records and the δ18O record could arise from the differences in the region that each record represents (see Methods), the depth at which those instrumental measurements are taken from (all are SST except for the EN4 at 35 m data product) and the sample quality and depth of the instrumental data themselves through time (each dataset uses different data sources, which vary in coverage spatially and temporally. See Methods for more details).

Table 1 Correlations between the oxygen isotope record and available observational records from 1983 to 2013.

Given the statistically significant correlations between the δ18O record with instrumental temperature records, we then analyze the oxygen isotope record along with the water mass source proxies, δ15N and Δ14C, to understand hydrographic changes in the Gulf of Maine through time. All three geochemical records show centennial scale variability since the late 17th century (Fig. 2).

Before the mid-1800s, the δ18O record is in broad agreement with the previously published Δ14C record, with decreasing δ18O (indicative of warming or freshening waters) and increasing Δ14C (potentially indicative of the increasing presence of WSW) from the beginning of the record to the early 1700s, followed by increasing δ18O (colder or saltier waters) and decreasing Δ14C (increased presence of LSW) until the mid-1800s. These shifts in water composition are accompanied by decreasing δ15N values at the beginning of that record (1750) until the early-mid 1800’s, likely indicating shifts in the nitrogen cycle in the Gulf of Maine.

During the mid to late 1800s, the region experienced a rapid shift in hydrographic properties resulting in relatively high δ18O (indicative of colder or more saline waters), high δ15N (potentially indicative of the increased presence of LSW), and low Δ14C (again potentially indicative of the increased presence of LSW)17 waters. This hydrographic composition peaks in the 1870s, as shown by segmented regression analysis on the δ18O record (which produces a breakpoint of 1875). From the 1870s until present day, waters in the region have gradually shifted towards a composition with larger contributions from waters with low δ18O, low δ15N, and higher Δ14C (this latter record only contains data up until 1935 due to atomic bombing testing in the 1950s), suggesting increasing WSW in the region. Notable exceptions from these general trends include decreased Δ14C values around 1900, suggesting increased LSW in the region.

To put these geochemical results into a broader temporal and spatial context, we analyzed fully-coupled climate model simulations from the Community Earth System Model-Last Millennium Ensemble (CESM-LME)38, which includes 13 different fully-forced ensemble member simulations in addition to several smaller single-forcing ensembles where only one external forcing (greenhouse gas, volcanic, solar, ozone/aerosol, land use change, and orbital) is allowed to vary through time39. The CESM-LME is forced by reconstructions of external forcings over the last millennium. Further details of these climate model simulations can be found in the Methods section. The ensemble member mean for the Gulf of Maine region, calculated as annual (January–December) temperature at 35 m depth, shows a gradual cooling in the region over the last millennium followed by a relatively rapid warming beginning at the turn of the 20th century (Fig. 3). Therefore, the decadal-to-centennial scale shift in water mass properties recorded by the shell-based isotopic records coincides with the end and reversal of a long-term cooling trend in the Gulf of Maine suggested by the climate model simulations. Specifically, the segmented regression analysis of the ensemble member mean, performed using one breakpoint (i.e., the year at which the largest change in trend throughout the record occurs), suggests this trend reversal occurs in 1904, timing which is corroborated by the majority of individual ensemble members and is slightly delayed but in general agreement with the changes in temperature suggested by the oxygen isotope record (where segmented regression analysis suggests a trend reversal in 1875; Fig. 4). Multi-decadal trends are even more similar between Gulf of Maine oxygen isotope records (both this study and Wanamaker, et al.15) and δ18O values calculated using temperature and salinity variables from the climate model simulations (Fig. 5 and Supplementary Fig. 3; see the Methods section for an explanation of this analysis), with each time series suggesting cooling Gulf of Maine conditions until around or just before the turn of the twentieth century, followed by rapid warming. Values are normalized in Fig. 5 for better comparison between records as the oxygen isotopes calculated using the climate model simulations demonstrate significantly lower variability than the oxygen isotope records (a common caveat of climate model simulations40). Differences in the beginning of the oxygen isotope record from this study and the modeled oxygen isotope record (i.e., for the period between ~1750 and 1800; see Fig. 5b) can likely be attributed to internal variability: As shown by the differences in the individual ensemble members in Fig. 5b, internal variability significantly influences the oxygen isotope time series.

Fig. 3: Simulated Gulf of Maine temperature and AMOC over the last millennium.
figure 3

Bold lines are ensemble member means while thin lines show each individual ensemble member for each simulation from the Community Earth System – Last Millennium Ensemble (CESM-LME). Trend lines show the results of the segmented regression analysis for the fully-forced Gulf of Maine temperature time series. See the Methods section of main text for full explanation of how each time series was calculated. Anomalies are calculated using the average over the entire record of each dataset. AMOC Atlantic Meridional Overturning Circulation. GHG greenhouse gas, VOLC volcano.

Fig. 4: Distribution of breakpoints of Gulf of Maine temperature and oxygen isotope time series.
figure 4

a Breakpoint years calculated from segmented regression analysis of single and fully forced simulations as well as the oxygen isotope record over the period of overlap between the records (1694–2005). The histogram plots the breakpoint year of ensemble members of the fully-forced simulations, with a bin-size of ~60 years. The oxygen isotope record breakpoint (green vertical line) is 1875 and calculated over the full record (1694–2013). The calculated breakpoint for each single forced ensemble member are represented by colored points on the x-axis.Years of ensemble member mean breakpoints (for both full forcing and single forcing simulations) are represented as vertical lines. b As in a but segmented regression analysis is performed over the entire millennium period (1000–2005). Bin size is still ~60 years. Oxygen isotope breakpoint included for reference. SSI spectral solar irradiance, VOLC volcanic, GHG greenhouse gas. Breakpoint years calculated from the climate model simulations are the result of segmented regression analysis of annual temperatures at 35 m depth in the Gulf of Maine for individual ensemble members (colored points; both full forcing and single forcing), the ensemble member mean (vertical lines; both full forcing and single forcing).

Fig. 5: Oxygen isotope records in the Gulf of Maine over the last 1000 years.
figure 5

a Time series of oxygen isotopes measurements from this study and Wanamaker, et al.15 compared to those calculated from temperature and salinity output from the fully forced CESM-LME simulations, averaged over the Gulf of Maine region at 35 m depth (see the Methods section in the main text for details). The dashed rectangle represents the time period of the graph shown in (b) for ease of comparison. Note both y axes are inverted. b As in (a) but for the shorter period of 1685–2013. The legend shown in (a) applies to both (a) and (b). All oxygen data are normalized (z-scores) so that they can be plotted on the same scale.

This most recent century (1901–2000) of warming in the Gulf of Maine is faster (0.44 °C century−1) than almost any other 100-year period in the past 1000 years (Fig. 6), as indicated by Monte Carlo simulations of the fully-forced Gulf of Maine climate model temperature simulations. The Monte Carlo simulations suggest the warming trend is greater than 99% of the trends of the ensemble member mean, randomly sampled 10,000 times. The slowest warming trend over the past 100 years for the individual ensemble members (0.15 °C century−1) is greater than 78% of the trends of the individual ensemble members. The fastest warming trend of the individual ensemble members (0.95 °C century−1) is greater than 99% of the trends.

Fig. 6: Linear regression analysis of temperature trends in the Gulf of Maine.
figure 6

Histograms show the distribution of slopes of linear regression analysis of 100-year periods in the 1000 year Gulf of Maine temperature simulation of the fully forced ensemble member mean (a) and individual ensemble members (b). The ensemble member mean and each of the individual ensembles was sampled 10,000 times using Monte Carlo simulations. Black vertical lines show the slope of the linear regression of the most recent 100 years of the Gulf of Maine temperature simulations for the ensemble member mean (a) and the individual ensemble members (b).

Single forcing model simulations of the Gulf of Maine region suggest that most of the long-term, pre-industrial cooling and multi-decadal variability seen in the fully forced simulation of average Gulf of Maine temperatures appears to be driven by volcanic forcings (Fig. 3; r = 0.56, p < 0.00001, n = 1000), with lesser influence from other forcings (Supplementary Fig. 4). This variability includes significant warming periods after cool periods triggered by volcanic eruptions, especially after 1260, where the 100-year trend is 0.44 °C century−1. This warming reflects the largest 100-yr trends in the Monte Carlo simulations, slightly larger than the warming trend seen in the Gulf of Maine in the last century (when considering more than two significant digits). An analysis of AMOC strength over the last millennium in the fully-forced climate model simulations also shows similar multidecadal and centennial trends to the Gulf of Maine temperature simulations (Fig. 3).

The only portion of the simulated Gulf of Maine temperature record that does not seem to be significantly influenced by volcanic forcings is the most recent period, since approximately 1850 (Fig. 3). The considerable trends of both greenhouse gases and aerosols during this period explain most of the trends in the fully forced output for the Gulf of Maine region (Supplementary Fig. 4). However, none of the single forcing simulations show the cold period in the late 1800s, present in a majority of the ensemble members from the fully-forced simulations (Supplementary Figs. 5 and 6) as well as the oxygen isotope record. The absence of an obvious external driver for this period is shown clearly by an analysis of extreme years of single-forcing simulations (Supplementary Fig. 6). We hypothesize that the single-forcing simulations do not have enough ensemble members to enable discernment of a forced signal separate from internal variability during this time period by looking at the single-forcing ensemble member mean41,42.

The resolution of the climate model we use in this study precludes an in-depth analysis of the circulation within the Gulf of Maine. Therefore possible changes in the transport of different water masses into the Gulf of Maine through time as interpreted in the isotopic records cannot be directly compared to the model. However, changes in temperature within the model can be compared to the isotopic records, as done in the preceding paragraph. Additionally, analyzing broader-scale North Atlantic conditions during extreme temperature years in the Gulf of Maine (defined as the top and bottom 10% of annual average temperatures; see the Methods for a description of how these analyses were performed) can inform our understanding of the relationship between Gulf of Maine conditions and changes in North Atlantic Ocean circulation. Such an analysis of the fully forced individual ensemble members (Supplementary Fig. 5) suggests that cold periods in the records, of which the cold period in the late 1800s dominates, correspond to anomalous upper-ocean (top 200 m) temperature, salinity, and current velocity patterns in the Gulf Stream region near the Gulf of Maine as well as significantly anomalous North Atlantic Oscillation (NAO)-like sea level pressure patterns in the North Atlantic (Fig. 7). Similar patterns are seen if only looking at ocean variables at the surface (not shown). A similar analysis using the ensemble member mean (i.e., the variability associated with external forcings) shows similar, although reduced, upper-ocean temperature, salinity, and ocean current anomaly patterns but does not show NAO-like sea level pressure patterns (Supplementary Fig. 7). Such anomaly patterns suggest a northward (southward) shift in the Gulf Stream during extreme warm (cold) years in the Gulf of Maine as can also be seen in the meridional temperature profiles (Supplementary Fig. 8). Similarly, the three geochemical hydrographic proxy records derived from A. islandica shells discussed in this study suggest a northward shift of the Gulf Stream over the 20th century as they suggest Gulf of Maine warming concurrent with increases in the proportion of WSW in the region.

Fig. 7: North Atlantic conditions during extreme temperature years in the Gulf of Maine.
figure 7

Simulated anomalous conditions in the North Atlantic during periods with extremely warm (ad) and cold (eh) temperature in the Gulf of Maine in the individual ensemble members of the climate model (averaged together for this figure) over the period of overlap with the geochemical records in this study (1685–2005): a, e: Temperature (upper 200 m average), b, f: Salinity (upper 200 m average), c, g: Sea level pressure (SLP), d, h: ocean current velocity (upper 200 m average). Stippling in ac and eg identify anomalies significant at the 95% significance level. Only vector anomalies significant at the 95% significance level were plotted in d and h.

Discussion

A comparison between the δ18O record and the climate model simulations reveals that the oxygen isotope record presented in this study captures the ending and then reversal of a long-term, multi-century cooling trend in the region. This initial cooling trend seems to be largely driven by volcanic forcings as suggested by the model simulations. This cooling over the last millennium followed by rapid warming is consistent with CESM-LME global SST (Supplementary Fig. 9) and CESM-LME Northern Hemisphere surface temperature simulations, which are also attributed to volcanic forcings prior to 1850 and greenhouse gas and aerosol forcings after 185039. Otto-Bliesner, et al.39 note that the CESM-LME simulations tend to exhibit a more significant amount of cooling after a volcanic eruption than temperature proxy reconstructions suggest occurred, pointing to uncertainties in climate responses to volcanic eruptions43. Supplementary Fig. 9 shows that while Gulf of Maine temperatures have similar trends to the global average SSTs in the CESM-LME model, the particularly cold period in the late 1800s in the Gulf of Maine is not present in the global average SSTs. This difference in relative temperature a century ago means that the Gulf of Maine has warmed much more significantly in the last 100 years than the average global surface ocean as the two records demonstrate similar temperature anomalies today in the model simulations.

The influence of volcanic activity on ocean surface temperatures globally and on those in the North Atlantic in particular over the last 1000 years has been previously documented. McGregor, et al.44 used 57 SST proxy records from marine archives to demonstrate a long-term cooling trend in global oceans over the last 2000 years which they attribute to volcanic forcings by comparing their synthesis data to climate model simulations from 801 to 1800. In particular, the authors propose that the increased frequency of volcanic eruptions in the first half of the last millennium might be driving the majority of this long-term cooling through a reduction in mixed-layer water temperatures. Mann, et al.45 also find significant influence of volcanic eruptions on SSTs, in particular for the North Atlantic. Through an analysis of CMIP5 simulations over the last millennium, the authors concluded that volcanic forcings were the primary driver of the previously identified Atlantic Multidecadal Oscillation (AMO), which describes a 40–60 year periodicity of SSTs in the North Atlantic.

The general cooling trend in the model simulations, for both the Gulf of Maine and, to a lesser extent, the average global surface ocean, from a warmer period in the beginning of the last millennium reflects the Medieval Climate Anomaly to Little Ice Age transition resolved by multiple climate model simulations (including the CESM-LME global surface temperature simulations39) and paleoproxy climate reconstructions and generally attributed to multiple significant volcanic eruptions in the first half of the last millennium. The cooling found in the current study may also in part be related to longer-term cooling found throughout the Holocene in the western North Atlantic surface waters that in part feed the Gulf of Maine46. This cooling of 4–10 °C over the Holocene found in three different locations through alkenone paleothermometry was hypothesized to be related to some combination of decreased solar insolation, increased Labrador Sea convection and/or a southward shift of the Gulf Stream.

As previously noted, Wanamaker, et al.15 also found cooling in the Gulf of Maine over the last 1000 years (Fig. 5). This cooling was determined by using oxygen isotopes measured in 4 different clam shells that spanned different periods of the last 1000 years and were collected near the study site for this study but at a shallower depth (30 m compared to 38 m depth for this study). While the Wanamaker et al. oxygen isotope record is not significantly correlated with the oxygen isotope record presented here, the two reflect similar trends of cooling in the Gulf of Maine up until the late 1800s/early 1900s, followed by warming. The lack of correlation is possibly related to a multitude of factors, including stratification influences and possible dating errors in the Wanamaker, et al.15 study, which was not crossdated.

In addition to long-term cooling trends, it is clear that the mid- to late-1800s were a time of dramatic change in the North Atlantic, as documented both in the Gulf of Maine geochemical records presented and discussed in this study, as well as other records of temperature and ocean circulation changes throughout the North Atlantic (Figs. 1 and 8). These include temperature records from the Scotian Shelf47, the Gulf of St. Lawrence12, and the western Labrador Sea48, as well as proxy records more directly attributed to ocean circulation changes, including off of Cape Hatteras49. Some of this variability has been associated with changes in atmospheric circulation patterns related to the beginning of atmospheric warming47,48, while others have been attributed to changes in the strength of the AMOC49, specifically a weakening of the AMOC as a result of anthropogenic climate change and/or the end of the Little Ice Age.

Fig. 8: Reconstructions of North Atlantic variability over the last 300 years.
figure 8

a Oxygen isotope record from this study (left y-axis, note inverted y-axis) and radiocarbon record (right y-axis) originally published in Lower‐Spies, et al.17. b Nitrogen isotope record from this study and Whitney, et al.26 (left y-axis) and a nitrogen isotope record from deep-sea corals collected outside of the Gulf of Maine and published in Sherwood, et al.25 (right y-axis). c The Labrador Sea algal time series (LATS) from Moore, et al.48, calculated by averaging the normalized average growth rates and Mg/Ca ratios of two different coralline algae samples collected from the western Labrador Sea. d Reconstructed SSTs from alkenones in sediment cores from Emerald Basin on the Scotian Shelf from Keigwin, et al.47. e Oxygen isotopes measured in benthic foraminifera from sediment cores collected from the Gulf of St. Lawrence from Thibodeau, et al.12. f Sortable silt records from sediment cores collected off of Cape Hatteras from Thornalley, et al.49. The locations of previously published records are marked in Fig. 1.

A tilefish die off in the Mid-Atlantic Bight, just south of the Gulf of Maine, in 1882, around the time that the geochemical records and model simulations in this study suggest the presence of the coldest waters in the region during the last 1000 years, has been attributed to significantly cold water related to increased LSW in the region as a result of a particularly negative NAO50. Certainly, atmospheric variability associated with the NAO has long been suggested as playing a dominant role in regulating water properties in this region50,51,52,53,54,55,56,57,58. The model simulations presented here do suggest some role in atmospheric circulation patterns, in this case specifically the NAO, during extreme temperature years in the Gulf of Maine.

These significant oceanographic changes at the end of the 1800s can be further assessed using the multiproxy isotopic analysis presented in this study. This analysis provides a more holistic approach to interpreting paleoceanographic conditions compared to using individual proxy records. While the δ18O record provides reconstructions of temperature and salinity, combining the δ18O record with the δ15N and Δ14C records has the added benefit of tracking changes in source water into the Gulf of Maine. For the 30 years prior to the 1870s, the three isotope proxy records taken together suggest that the proportion of LSW in the Gulf of Maine was increasing17. Comparison between the Δ14C record and other records in the western North Atlantic region47 suggests decreased SSW into the region during this time period, particularly during an abrupt change in the Δ14C record in the 1860s47. Decreased SSW in the Gulf of Maine would allow for an increase in slope water to enter at depth. The δ18O record suggests there was a steady cooling in the Gulf of Maine during this time period, with water temperatures in the early 1800s warmer than during the turn of the 20th century, suggesting that increased slope water in the region was likely mainly composed of LSW.

This multiproxy approach also benefits the interpretation of other portions of the geochemical records. Earlier in the records, between the early 1700s and the mid-1800s, the oxygen isotope and radiocarbon records suggest a cooling and possible increase in LSW through time, while the nitrogen isotope record reflects decreasing nitrogen isotope signatures of dissolved nitrate. Without the oxygen and radiocarbon records, the nitrogen isotope record during this time period might be thought to indicate increasing WSW in the region. However, the other isotopic records do not suggest that this is possible, nor do the climate model simulations. We therefore put forth an alternative hypothesis and suggest that decreasing nitrogen isotopic composition of dissolved nitrate during this time period might be indicative of local mixing processes. Specifically, a decrease in mixed-layer temperatures as suggested by the oxygen isotope record and the model simulations might lead to a better-mixed and consequently more oxygenated water column, which in turn would result in dissolved nitrate at the surface becoming more depleted in 15N through time as (1) dissolved nitrate with decreased nitrogen isotope signatures from depth is mixed to the surface and (2) nitrate at depth undergoes less denitrification (which would increase the nitrogen isotopic signature of dissolved nitrate) due to increased oxygenation of the water column.

Towards the more recent portion of the reconstructions, the similarity in trend between the δ18O, δ15N, and Δ14C records suggests that the warming seen over the last century was driven in part by potential changes in the proportion of water masses entering the Gulf of Maine, with the proportion of WSW in the Gulf of Maine increasing since the 1870s. This increase in WSW suggests changes in North Atlantic ocean dynamics over the last ~150 years, possibly as a result of increased greenhouse gases as suggested by the climate model simulations.

This analysis is limited to some extent by the fact that the radiocarbon record does not extend to present day due to the bomb pulse in the 1950s. While the radiocarbon record in the beginning of the 20th century is consistent with an increase of WSW over the last century, any interpretation of ocean dynamic changes during the later half of the 1900s based on this record alone is not possible. We note here that it is also possible that the decrease in nitrogen isotopes found in both the Gulf of Maine record as well as the Sherwood, et al.25 record of corals from outside of the Gulf of Maine (Fig. 8b) could be related to anthropogenic atmospheric nitrogen deposition26,59,60. Further data are needed from source waters of the Gulf of Maine to understand potential drivers in the nitrogen isotope signatures of these waters through time.

The warming in the Gulf of Maine since the dramatic cooling at the end of the 19th century is likely caused by several potential factors, including changes in North Atlantic Ocean dynamics suggested by the geochemical reconstructions and specifically more recent changes in the Gulf Stream position as discussed in the Introduction. Climate model simulations of the region suggest that greenhouse gas forcings have largely driven the warming since the breakpoint in the record at 1904 (Supplementary Fig. 4). However, the greenhouse gas-forced warming has been dampened by the ozone and aerosol forcing over this same period as those forcings contribute significant cooling in the region, as seen in the respective single-forcing experiments. The strength of the greenhouse gas and ozone/aerosol forcings means that the volcanic forcing contributes much less to Gulf of Maine temperature variability after roughly 1800.

Several studies have suggested that some of the dramatic hydrographic changes seen in the western North Atlantic near the Gulf of Maine, particularly the warming in the Gulf of St. Lawrence12 and the decrease in nitrogen isotopes outside of the Gulf of Maine25, have been caused by a significant weakening of the AMOC over the last 50–150 years25,61,62. However, a weakened AMOC in the 21st century is still debated63,64, particularly as there are many other proxy records that do not suggest a recent weakening of the AMOC16.

The climate model’s AMOC evolution in this study suggests that the AMOC strengthened in the mid-1800s before weakening again, in agreement with the isotope records presented here. Figure 3 shows that the modeled AMOC over much of the last 1000 years shows similar patterns to the Gulf of Maine temperature time series from both the fully-forced simulations as well as the volcanically forced simulations. This relationship is likely the result of the influence that volcanic forcings have on both Gulf of Maine temperatures as well as AMOC strength65,66 and is not necessarily an indication that AMOC is driving the Gulf of Maine temperature changes in the climate model.

However, it seems probable that changes in AMOC strength amplified volcanically forced changes in Gulf of Maine temperatures due to the relationship between AMOC variability and the position of the Gulf Stream57,67,68,69. The position of the Gulf Stream in part dictates the proportion of slope water entering the Gulf of Maine67,70. When the Gulf Stream is in a more northerly position (weaker AMOC), and therefore closer to the Gulf of Maine, more WSW is brought to the Northeast Channel, either because there is more WSW in the Slope Sea region (directly outside of the Northeast Channel)70 or because more WSW is transported to the Tail of the Grand Banks (south of Newfoundland) and transported to the Gulf of Maine via southward flowing currents10,71. The analysis of extreme events in the climate model simulations clearly shows that the position of the Gulf Stream plays a significant role in dictating temperatures within the Gulf of Maine (Fig. 7 and Supplementary Figs. 7 and 8). Increased WSW in the Gulf of Maine suggested by the geochemical records presented here could be a result of a northward shift in the Gulf Stream over the last 100 years, although other changes in ocean dynamics cannot be ruled out.

Of note, the simulations of AMOC strength analyzed here do not suggest that AMOC has weakened significantly in recent years as a result of greenhouse gas forcings. These findings, therefore, suggest that the warming shown in the simulations for the Gulf of Maine is largely driven by increased greenhouse gas forcings influencing other components of regional ocean conditions. The similar trends between the model simulations and the isotope and instrumental records point to recent warming trends suggested by these records also being driven primarily by increased greenhouse gases. However, other primary drivers cannot be excluded.

Regardless of whether the AMOC has begun to weaken, most models suggest that the AMOC will weaken in the future72,73,74. Saba, et al.9 demonstrate, using Geophysical Fluid Dynamics Laboratory (GFDL) global climate models, that the robust relationship between AMOC and Gulf of Maine temperatures will lead to extreme warming in the Gulf of Maine under increased CO2 conditions as AMOC weakens (the upper ocean warmed three times faster than the average global ocean under a CO2 doubling scenario). Thus, in combination with ongoing warming in the Gulf of Maine, likely resulting from increased greenhouse gas forcing, and resulting changes in ocean dynamics as reported in this study, along with the other more recent changes noted in Gulf Stream behavior6,7,10, Gulf of Maine temperatures will also be increasingly influenced by a weakening AMOC. The simulated temperatures presented here suggest that the Gulf of Maine has warmed faster in the last 100 years than any other period in the last 1000 years, further supported by the long-term trend in the BBH SST record. Given future projections of atmospheric greenhouse gas concentrations and AMOC strength, this warming trend in the Gulf of Maine is likely to continue, leading to continued and potentially worsening ecologically and economically devastating temperature increases in the region in the future.

Methods

Proxy archive and geochemical reconstructions

Arctica islandica shells (Supplementary Fig. 1) are a valuable marine proxy archive as they are long-lived75,76,77, precipitate their shell in annual increments78, and have been shown to faithfully record environmental conditions8,22,26,75,79. A. islandica precipitate their shell in oxygen isotopic equilibrium with seawater, thereby recording seawater temperature78,80, and reliably record both nitrogen isotopic composition of dissolved nitrate in both the carbonate and periostracum layer of the shell26 and radiocarbon isotopic composition of dissolved inorganic carbon in seawater in the carbonate of the shell22,81,82. A. islandica shells used in this study were collected from near Seguin Island (western Gulf of Maine; 43.7°N, 69.8°W) at 38 m depth in 2010, 2011, 2012, and 2014. The site was chosen in part because of its proximity to the longest continuous instrumental record in the Gulf of Maine (the Boothbay Harbor SST record referred to throughout this text) as well as a site used in a previously published shell-based oxygen isotope reconstructions15 and as part of a larger network of sampling locations. 34 shells were crossdated into a chronology extending from 1761 to 2013, with several shells individually dated before 1761 (Supplementary Fig. 2). Details on processing and dating of A. islandica shell material, which includes the development, via crossdating techniques83,84, of a master chronology, can be found in Griffin28 and Wanamaker, et al.8. Briefly, shells were cleaned with deionized water, blocked in epoxy and sectioned along their axis of maximum growth. Shell blocks were then polished using a Buelher variable speed grinding-polishing wheel. One block was then etched using 1% HCl and an acetate peel made and imaged for dating purposes, and the other block was reserved for isotope analysis sampling. All processing and analysis took place at the Stable Isotope Lab at Iowa State University unless otherwise noted. While the chronology only extends back to 1761 as the express population signal (EPS) drops below 0.85 before this date8, individual shells have been dated with a high degree of confidence farther back in time. Some of the proxy records presented in this study therefore extend back before 1761.

For oxygen isotope analyses, individual annual increments of shells were milled on the shell block (Supplementary Fig. 1) using a Mechantek micromill with a Leica GZ6 microscope and Brasseler USA carbide drill bits. The resulting carbonate powder samples weighed between 100 and 300 μg. Samples were placed in glass vials and put in a sample oven at 50 °C to burn off any water. Vials were then capped and flush filled with helium gas before samples were acidified with phosphoric acid (H3PO4). Samples were analyzed for oxygen isotopes on a ThermoFinnigan Delta Plus XL mass spectrometer interfaced with a GasBench II and CombiPal autosampler. NBS18, NBS19, and LSVEC were used as reference standards to correct samples using a regression method, resulting in oxygen isotope data standardized to VPDB (Vienna Pee Dee Belemnite). Results are reported in delta notation. The long term, average 2σ analytical uncertainty for δ18O is ±0.18‰. δ18O values measured in increments dated to the same year from different shells were averaged together for the final time series presented in this study. Sample depth in the time series ranged from one to four shells each year. The annually-resolved record is continuous after 1769.

Nitrogen isotope analyses were performed using periostracum material, which is the quinone-tanned, protein rich layer on the outside of the shell (Supplementary Fig. 1). Details on the sampling of this material as well as the validation of periostracum as a useful proxy for nitrogen isotopes of dissolved nitrate in water, can be found in Whitney, et al.26. Additionally, Whitney, et al.26 demonstrate through compound specific isotope analysis of amino acids that variability in this nitrogen isotope record is likely not related to changes in trophic position. The amount of time integrated by each sample varied. Samples were analyzed for bulk nitrogen isotopes at SIL on a ThermoFinnigan Delta Plus XL mass spectrometer interfaced with an Costech Elemental Analyzer. Caffeine (IAEA-600), IAEA-N2, Cellulose (IAEA-CH-3), and Acetanilide (laboratory standard) were used as reference standards to correct samples using a regression method. The resulting nitrogen isotope values were standardized relative to Air (atmospheric N2). Results are reported in delta notation. The long-term, average 2σ analytical uncertainty for δ15N was ±0.36‰. δ15N values measured in increments dated to the same mid-point year from different shells were averaged together for the final record. Sample depth of the time series ranged from one to three shells.

Details on sample collection for radiocarbon analysis can be found in the original manuscript in which it was published17. Briefly, carbonate material was collected from either the outer face of the shell or from the shell block (Supplementary Fig. 1). Samples were sent to the National Ocean Sciences Accelerator Mass Spectrometry (NOSAMS) facility at Woods Hole Oceanographic Institution in Woods Hole, Massachusetts and sampled using Accelerator Mass Spectrometry. Radiocarbon data are presented here as age corrected Δ14C, calculated using the absolute age of the sampled shells as determined through chronology building described above. The equation for Δ14C, as defined by Stuiver and Polach85, is:

$$\triangle {\!\,}^{14}C=(Fm*{{{{{\rm{e}}}}}}^{\lambda (y-x)}-1)\times 1000$$
(1)

Fm is fraction modern, which is equal to the percent difference of a sample compared to 95% of the radiocarbon concentration of the standard NBS Oxalic Acid I, NBS Oxalic Acid I having been normalized to δ13C = −19‰86. λ is defined as 1*(true mean half-life)−1, which is 1.2096 × 10−4 for radiocarbon. y is CE 1950 and x is the year that the shell material formed. The average 2σ analytical uncertainty in Δ14C was ±5.8‰. Sample depth of the time series was no more than 1 shell.

The δ18O record was compared to the Boothbay Harbor sea surface temperature instrumental record, which is located ~16 km away from the shell collection site (43.8°N, 69.6°W). The record begins in March 1905, so correlations were calculated starting in 190630. Gridded SST datasets, where observations that have been interpolated using varying techniques onto rectilinear grids, were also compared with the δ18O record. These datasets included the Extended Reconstructed Sea Surface Temperature (ERSST) v5 beginning in 185431, Hadley Center Sea Ice and Sea Surface Temperature dataset (HadISST) v1.1, starting in 187032,33, Optimum Interpolation Sea Surface Temperature (OISST) v2.1 starting in 198234,35, as well as surface and 35 m depth temperatures from the EN4.2.1 data product36,37. Although the EN4 records extends back to 1900, data uncertainty and weights associated with this dataset suggest the data are sparse in this region prior to 1960. Therefore, we only compare the post 1960 EN4 dataset to our δ18O record. Correlations between the various temperature datasets and the δ18O record were performed by averaging the entire Gulf of Maine area for each gridded dataset. Because each dataset is on a slightly different rectilinear grid, the area over which the average was calculated varied slightly: HadISST (42.5°N-46.5°N, 70.5°W-66.5°W, 1° horizontal resolution), ERSST(42.0°N-46.0°N, 70°W-66°W, 2° resolution), OISST(41.875°N-46.625°N, 70.875°W-65.875°W, 0.25° resolution), EN4 (42.0°N-46.0°N, −71.0°W-66.0°W, 1° resolution).

To help place the geochemical proxy timeseries into spatial and temporal context, we used simulations from the Community Earth System Model-Last Millennium Ensemble (CESM-LME)38,39. The CESM is a fully-coupled, global climate model developed by the National Center for Atmospheric Research (NCAR). The LME is a set of experiments performed by the CESM1 Paleoclimate Working Group at NCAR. The LME uses CESM1 version 1 at a ~1° resolution in the ocean using the sea-ice version of the Parallel Ocean Program version 2 (POP2) and a 1.9° latitude by 2.5° longitude resolution for the atmosphere based on the Community Atmosphere Model version 5 (CAM5). 13 ensemble members have full forcings by orbital, solar, volcanic, greenhouse gas, aerosol/ozone, and land use/land cover changes over the last millennium. Additional simulations are forced by a single varied forcing with the remaining forcings held constant at their 850 level, except for the ozone and aerosol forcing, which is kept at its 1850 level and the volcanic forcing, which is not included in single-forcing runs except the volcano-only simulations. The greenhouse gas single-forcing runs have three ensemble members, the volcanic single-forcing runs have five ensemble members, the solar single-forcing runs have four ensemble members, and the ozone/aerosol single forcing runs have five ensemble members. Land use/land cover and orbital single forcing runs are not discussed in this study because they had little to no influence on these time scales on Gulf of Maine conditions. Additional information on the CESM-LME simulations can be found in Otto-Bliesner, et al.39.

Analysis done for this study employed both the 13 fully-forced ensemble members as well as the single forcing simulations of volcanic, greenhouse gas, solar, and ozone/aerosol forcings. All CESM-LME output used in this study was first regridded from a curvilinear grid to a 1° × 1° rectilinear grid. Seawater temperatures and salinities at 35 m depth (the closest depth to where the clams were collected at 38 m depth) in the Gulf of Maine region (average of data in the region bounded 41.5°N, 44.5°N, 73.5°W, and 66.5°W) were annually-averaged. To assess the changes in the Gulf of Maine over the last 1000 years, the period between 996 and 2005 C.E. was selected for both the fully forced and single forcing runs and an 11-year running average of the data was performed. The initial period of 996–2005 C.E. was chosen so that the 11-year running average represented the full millennium: 1001–2000 C.E. For analyses where climate model simulations were directly compared to isotope records (e.g., for the breakpoint analysis and calculated oxygen isotopes), a running average was not calculated.

To calculate what would be expected for oxygen isotopes values in shells given the conditions in the model simulations, the CESM-LME temperature and salinity time series from the Gulf of Maine, each calculated as described above, were used, along with the linear relationship between salinity and oxygen isotopes of seawater (δ18Oseawater) in the Gulf of Maine as determined by Whitney, et al.18:

$${{{{{{\rm{\delta }}}}}}}^{18}{{{{{{\rm{O}}}}}}}_{{{{{{\rm{seawater}}}}}}}=0.5* {{{{{\rm{Salinity}}}}}}-17.3$$
(2)

Calculated δ18Oseawater along with the temperature time series was input into the modified equation from Grossman and Ku87 to produce a calculated oxygen isotope time series:

$${{{{{\rm{Temperature}}}}}}({\deg}\!{{{{{\rm{C}}}}}})=20.60-4.34{{{{{\rm{x}}}}}}[{{{\updelta }}}^{18}{{{{{{\rm{O}}}}}}}_{{{{{{\rm{shell}}}}}}}-({{{\updelta }}}^{18}{{{{{{\rm{O}}}}}}}_{{{{{{\rm{seawater}}}}}}}-0.27)]$$
(3)

The AMOC from the CESM-LME was calculated by choosing the eulerian mean of the Atlantic Ocean region, taking an annual average (January–December), and then selecting the maximum overturning stream function of each year for the region between 20°N and 90°N.

To assess what Gulf of Maine conditions might indicate about the broader North Atlantic region, the top and bottom 10% of annual temperatures from the fully forced Gulf of Maine time series (processed from the climate model simulations as described in the preceding paragraphs) were calculated for the period of overlap with the geochemical records (1685–2005 CE before an 11-year running average was applied). Conditions outside of the Gulf of Maine in the climate model simulations were then investigated for those identified extreme years. Sea level pressure and temperature, salinity, and meridional and zonal current velocity averaged over the upper 200 m were all investigated. Like for the Gulf of Maine time series, all model simulation output for the North Atlantic region was annually averaged and a 11-year running average was applied. The above analyses were done on both the ensemble member mean, to investigate the forced component of the signal, as well as the individual ensemble members (only averaged at the end of the analysis for visualization purposes as in Fig. 7), to investigate internal variability.

Statistics

The correlations between δ18O (the only annually-continuous geochemical record presented here) and the various temperature datasets based on instrumental data (BBH, ERSST, HadISST, OISST, and EN4) were assessed using simple linear regression optimized with least squares, Pearson’s r correlation coefficients, and associated p values (significant at p < 0.05), after accounting for autocorrelation using AR(1) modeling (using the Pyleoclim 0.6.2 package in Python). Linear interpolation was used to fill in infrequent missing monthly BBH data (17 months of missing data in the 107 year record). For the purposes of identifying periods of maximum correlations, including investigating leads and lags, correlations were performed over the period of overlap, discounting the first year so that the temperature leading the oxygen isotope record could be assessed. Pearson correlation coefficients were also used to compare the fully-forced simulated Gulf of Maine temperature record to single forcing Gulf of Maine temperature records.

Segmented regression analysis was used to determine breakpoints for changes in trend for both the δ18O data as well as the CESM-LME fully forced ensemble member mean (averaged from October-September at 35 m depth for the Gulf of Maine region). This analysis employs simple linear regression optimized with least squares to find the two regression lines that incorporate all of the data and produce the least errors. The place in the time series where the change in regression line occurs is termed the breakpoint and can be interpreted as the point where climatic trends change. One breakpoint with two linear regression lines was specified for both the δ18O data and climate model simulated time-series’.

Monte Carlo simulations were performed on the Gulf of Maine temperature time series from the climate model simulations to evaluate how unusual the recent warming trend in the Gulf of Maine is in the context of the last 1000 years. For both the fully-forced ensemble member mean and each of the individual, fully-forced ensemble members, 100-year periods were selected at random from the 1000 year record, and a simple linear regression analysis was performed. The Monte Carlo simulations were performed on the time series up to the 1900–1999 C.E. period but not including 1901–2000 C.E. so that this latter period could be evaluated in the context of the Monte Carlo simulations.

To determine the significance of the climate properties in the North Atlantic during extreme years in the Gulf of Maine, bootstrapping with replacement was used to compare the conditions during extreme Gulf of Maine years to all years in the record over the time period analyzed (1685–2005). When assessing North Atlantic condition anomalies during extreme years of each individual ensemble, significance of those North Atlantic anomalies when averaged together was determined by first stacking the ensemble members to make a long time series composed of all 13 ensemble members and then performing bootstrapping with a replacement on those ensemble members. Significance was determined at the 95% significance level.