Anthropogenic acidi ﬁ cation of surface waters drives decreased biogenic calci ﬁ cation in the Mediterranean Sea

Anthropogenic carbon dioxide emissions directly or indirectly drive ocean acidi ﬁ cation, warming and enhanced strati ﬁ cation. The combined effects of these processes on marine planktic calci ﬁ ers at decadal to centennial timescales are poorly understood. Here, we analyze size normalized planktic foraminiferal shell weight, shell geochemistry, and supporting proxies from 3 sediment cores in the Mediterranean Sea spanning several centuries. Our results allow us to investigate the response of surface-dwelling planktic foraminifera to increases in atmospheric carbon dioxide. We ﬁ nd that increased anthropogenic carbon dioxide levels led to basin wide reductions in size normalized weights by modulating for-aminiferal calci ﬁ cation. Carbon ( δ 13 C) and boron ( δ 11 B) isotopic compositions also indicate the increasing in ﬂ uence of fossil fuel derived carbon dioxide and decreasing pH, respectively. Alkenone concentrations and test accumulation rates indicate that warming and changes in biological productivity are insuf ﬁ cient to offset acidi ﬁ cation effects. We suggest that further increases in atmospheric carbon dioxide will drive ongoing reductions in marine biogenic calci ﬁ cation in the Mediterranean Sea.

S ince the onset of the Industrial Revolution, human-induced climate change has rapidly modified the Earth system with notable impacts on the oceans. Rapidly increasing atmospheric carbon dioxide (CO 2 ) concentrations have resulted in global surface ocean pH decline of approximately 0.1 units compared to pre-Industrial levels 1 , impacting the calcification process of marine calcifiers 2,3 . Since the mid-19th century, enhanced CO 2 emissions have also contributed to global ocean surface warming by 0.76 ± 0. 16 K between 1904 and 2016 4,5 . In the second half of the 20th century, anthropogenically induced warming has also intensified vertical ocean stratification 6 , leading to changes in nutrient cycling and depletion of surface marine productivity 7 . These alterations of pH, temperature, and food availability are each assumed to influence the calcification rates of critical calcifying taxa such as planktic foraminifera 8,9 , but their response to the combination of these factors is still unclear.
The Mediterranean is a primary climate change hot-spot that shows amplified climate responses to global change 10,11 , making it a key region to investigate the response of marine ecosystems to human-induced changes of the chemical, physical, and biological environment. Due to its high alkalinity and the short residence times of its water bodies, the Mediterranean Sea is expected to have an exceptionally high uptake rate of anthropogenic carbon [12][13][14] . Both annual mean surface air temperature and sea surface temperature (SST) have accelerated in the Mediterranean region compared to the global average since 1880 Common Era (CE) 11 . The recent SST warming at a rate of 0.35 K per decade 15 , may result in thermally induced stratification and decreasing biological productivity 16 .
Comprehensive, multi-decadal monitoring of plankton communities started in 1931 17 , while time series measurements of seawater carbonate chemistry 18 and surface ocean productivity 19 began only in the late 20th century. As a result, observational time series are too short to reasonably attribute changes in the ocean's biological productivity to anthropogenically driven climate change 20 . Laboratory experiments give important clues on how calcifying plankton copes with rapid climate change, but it remains difficult to capture the complex natural environment in a laboratory setting. Paleo-reconstructions of plankton calcification and accompanying records of environmental conditions can address these issues, providing records going back to preindustrial times which could capture the impact of changing physical and chemical conditions on the marine ecosystem.
Here, we provide radiocarbon and radionuclide-dated, highresolution marine proxy records from three locations in the western Mediterranean Sea (Fig. 1), covering the transition from the pre-industrial (pre-IE) to industrial era (IE), starting from about 1800 CE. With these we investigate the basin-wide, interspecies (Globigerinoides elongatus and Globigerina bulloides) response of planktic foraminiferal calcification (represented by the measured Size Normalized Weight, SNW) to warming, acidification, and variations in primary productivity, reconstructed using a suite of paleo proxies, over the onset of the Industrial Revolution. We show an unprecedented basin wide SNW decrease of 18% to 34% in the species G. elongatus, and 7% to 35% in G. bulloides over the 20th century. Enhanced uptake of anthropogenic CO 2 by Mediterranean surface waters as shown by decreasing planktic foraminiferal δ 13 C (Suess Effect) and decreasing pH derived from boron isotopes (δ 11 B), brings us to the conclusion that anthropogenic CO 2 impairs foraminiferal calcification.

Results and discussion
Human induced sea surface warming and productivity changes in the western Mediterranean. The western Mediterranean Sea is a hotspot of rising SST, at a rate above the global average 21 , and anthropogenically-driven productivity changes during the IE 16 . Warming of western Mediterranean surface waters over the 20th century is shown by alkenone derived SST data and (HadISST) satellite time series observations (Figs. 2a, 3c). A mid-20th century cooling period by > 1 K in the Mediterranean from the end of the 1940s to the end of the 1980s is confirmed, as seen in other western Mediterranean paleotemperature records 16,22 . This cooling episode is followed by rapid warming over the last 20 years. Any offsets between instrumental and paleorecords might be due to low data coverage for SST calibrations in the central Mediterranean (Supplementary Note 1). Indeed, while the Strait of Sicily alkenone record indicates warming just above average over the last two decades, other recently published paleorecords in the Alboran and Balearic Sea 16 consistently show accelerated 20th century sea surface warming (Fig. 2a), with unprecedented warming rates at present in the Mediterranean Sea 23 . High sea surface warming rates at our sites confirm a rapid, anthropogenically induced sea surface warming in the western Mediterranean.
Accelerated anthropogenic climate change has altered marine productivity in the Mediterranean Sea, as reflected by plankton community changes 16 . The sedimentation rate of planktic foraminiferal tests is positively related to primary productivity 24 , while alkenone concentrations are closely related to phytoplankton biomass 25 . Both foraminiferal tests and alkenones thus serve as proxies of changes in past ocean productivity and reflect surface productivity levels at all three study sites (Supplementary Note 1). In the Strait of Sicily, test accumulation rates decrease between the 1850s and 1940s, then increase rapidly, which, alongside increasing C 37 alkenone concentrations, indicates enhanced productivity during the second half of the 20th century (Fig. 2b). Biogeochemical model results 26 and historical records 27 , confirm productivity increase in Mediterranean surface waters during the 20th century, and suggest enhanced anthropogenic nitrate and phosphate loads as the main driver of surface water eutrophication. In contrast, a decrease in productivity is observed in the Alboran and Balearic records during the second half of the 20th century (Fig. 2b). At these sites, enhanced stratification through 20th century SST increase likely explain surface productivity losses, which are also recorded through changes in planktic foraminiferal species composition 16 . Our results and previous studies thus show regionally distinctive patterns of 20th century productivity changes for the western Mediterranean.
Anthropogenic CO 2 intrusion into western Mediterranean waters recorded by planktic foraminiferal test chemistry. The δ 13 C and δ 11 B signatures of planktic foraminiferal tests provide sensitive tracers of anthropogenic carbon uptake and acidification of the western Mediterranean Sea. While fluctuations in carbonate chemistry, productivity, and hydrography may modulate the carbon isotope signal, these influences are relatively small compared to the change in the isotopic composition of ambient seawater due to the input of low δ 13 C CO 2 from fossil fuel burning (the Suess Effect) [28][29][30] , making δ 13 C a good tracer for natural versus anthropogenic environmental changes.
We use the two morphotypes of planktic foraminiferal species Globigerinoides elongatus (syn. Globigerinoides ruber white sensu lato) and G. ruber albus (syn. G. ruber white sensu stricto) as well as Globigerina bulloides to evaluate environmental alterations in the upper water column, as they are surface mixed layer dwellers 31 . As G. bulloides thrives in late winter/early spring and G. ruber in late spring/summer 32,33 , we also cover most of the seasonal change in sea surface temperature and primary productivity.
At all three sites analyzed here, G. elongatus and G. bulloides show a notable δ 13 C drop during the IE, against a backdrop of pre-industrial natural variability (Figs. 2c and 3d), which we refer to here as negative anthropogenic carbon isotope excursions Fig. 2 Changes of sea surface conditions and planktic foraminiferal size normalized weight (SNW) over the Common Era at three western Mediterranean core sites (Alboran and Balearic Seas, and Strait of Sicily). a Sea surface temperature (SST) changes derived from alkenones (dark blue line) and the HadISST data set (light blue line). b Marine surface productivity represented by total test accumulation rate (purple line), alkenone concentration (grey symbols), and species composition (G. bulloides versus G. inflata) and species abundance (G. inflata + G. truncatulinoides) changes (brown line) 16,122 . c Suess effect displayed through δ 13 C (‰ VPDB) signature measured in planktic foraminiferal tests of G. elongatus (red dashed line) and G. bulloides (blue dashed line) with error bars indicating 0.1 ‰ uncertainty (1 SD level) and displayed as 3 point running average. d Surface ocean pH estimated with CO2SYS (black line-see methods) and boron isotope measurements (Orange symbols) in tests of G. elongatus (Alboran Sea) and G. ruber albus (Strait of Sicily), with vertical error bars indicating one standard deviation of uncertainty and horizontal error bars displaying age error range. e Planktic foraminiferal calcification expressed as SNW for G. elongatus (red) and G. bulloides (blue) with error bars indicating one standard deviation of uncertainty, displayed as 3 point running average. Grey shaded area in Alboran and Balearic time series represents the pre Industrial Era, prior to 1800 CE with a compressed time axis. (nCIEs). nCIEs show most marked decreases during the latter half of the 20th century with a mean amplitude of −0.60‰ and range between −1.24‰ and −0.14‰ ( Supplementary Fig. 1). We attribute these nCIEs to the increase of the surface ocean concentration of anthropogenic CO 2 34 . The measured decline of δ 13 C in our foraminiferal tests during the second half of the 20th century describes a more pronounced Suess Effect signal when compared with results from the Northwest Atlantic (range = −0.174‰ to −0.757‰; weighted average = −0.45‰) 30 .
Ocean uptake of anthropogenic CO 2 increases hydrogen ion (H + ) concentrations and hence lowers pH 35 . Boron isotope ratios of tests of G. ruber albus and G. elongatus reveal a notable pH decline in the Alboran Sea of −0.14 units between 1697 and 1995 CE, which is in accordance with previous estimates over the industrial period ranging from −0.06 to −0.16 units 36 , and also notable ocean acidification in the Strait of Sicily, with reconstructed pH decreasing by 0.22 units from 1830 to 2007 CE. These results are largely consistent with the acidification expected from anthropogenic CO 2 invasion ( Supplementary  Fig. 2), which we calculate by assuming equilibrium between surface waters and the atmosphere alongside modern alkalinity and observed SST change at each site (see Methods); we note that absolute pH values from δ 11 B are lower than recent pH levels ( Whilst the nCIEs primarily reflect uptake of isotopically light CO 2 in the surface ocean, falling [CO 3 2− ] and pH levels tend to increase δ 13 C in foraminiferal tests 44 , and could thus dampen the Suess Effect driven nCIE signal we observe. Based on laboratory and field data [44][45][46] , western Mediterranean pH and [CO 3 2− ] changes since 1900 are estimated to account for approximately 0.24‰ to 0.33‰ (G. ruber) and 0.48‰ (G. bulloides) of δ 13 C test changes (Supplementary Table 1). Ambient water temperature changes have a minimal, negative effect on δ 13 C in foraminiferal tests 45,47 , accounting for −0.12 ‰ in G. bulloides tests, when basing estimations upon SST warming rates in the western Mediterranean (Supplementary Table 1). Variability in upwelling strength and biological productivity can also influence δ 13 C in calcite tests at the species scale 45 , as can algal or bacterial symbionts, which change the microenvironment of planktic foraminifera near the test surface impacting their δ 13 C signal [48][49][50] . Strong pH and SST changes, differential upwelling, and productivity patterns at our three study sites are thought on the whole to slightly counter the Suess Effect and attenuate the full magnitude of the signal of decreasing δ 13 C in our planktic foraminiferal tests 29 .
We conclude that the input of isotopically light anthropogenic CO 2 in the surface ocean has caused nCIEs during the IE in both G. bulloides and G. elongatus in the western Mediterranean Sea.  At the same time, enhanced CO 2 uptake by the ocean drives acidification, as indicated by boron isotope measurements and model estimates, attenuating the Suess Effect signal. Sea surface warming, upwelling, and changes in marine surface production during the 20th century have minor attenuating influences on the δ 13 C signal.
Anthropogenic impacts on planktic foraminiferal calcification. SNW is a useful tool for understanding the differential response of planktic foraminifera to climate change stressors like ocean warming, acidification, and productivity changes. SNW provides a normalized measure of chamber wall thickness and density, reflecting carbonate production and preservation of planktic foraminiferal tests 47 , and is used as a proxy for past [CO 3 2− ] 9,51 , though has also been suggested to provide a potential proxy for growth rate in planktic and benthic foraminifera 52,53 . Although secondary influences such as physico-chemical conditions of ambient seawater, and autecological properties can influence SNW, these do not seem to prevent the application of the proxy 8,54 . We found that SNW of G. elongatus and G. bulloides decreased at all three stations throughout the last 200 yr., unprecedented against the backdrop of our~1500 yr. long records in the Alboran and Balearic Seas (Figs. 2e; 3f). Basin-wide relative SNW loss during the 20th century in G. elongatus ranges from 18% to 34% and in G. bulloides from 7% to 35% ( Supplementary  Fig. 3).
Decreasing ocean pH and [CO 3 2− ] due to enhanced CO 2 uptake from the atmosphere can impede the biological precipitation of carbonate tests 55 , depending on the speciesspecific response. Commonly, G. elongatus and G. bulloides reduce their test calcite mass with lower [CO 3 2− ], calcite saturation (Ω), and pH, shown by numerous studies based on sediment and water samples, sediment trap time series, as well as culture experiments 9,51,56,57 . Despite this, in the natural environment, it is important to consider the combined influence of the various environmental drivers on planktic foraminiferal calcification, with studies using cultures, sediment traps, plankton nets and marine sediments attributing SNW changes to surface ocean carbonate chemistry, alongside temperature and productivity 8,9,51,54,[57][58][59][60] .
Sea surface warming decreases the solubility of atmospheric CO 2 , elevating [CO 3 2− ] in surface water, while higher temperature can also foster enzymatic activity, resulting in faster growth rates and enhanced calcification rates 61 . Both effects cause an increase in planktic foraminiferal SNW, a relationship observed in our target species in a number of earlier studies 8,9,51,56,58,62 . Sea surface warming in the western Mediterranean, as shown in our records (Figs. 2a; 3c), should have enhanced SNW in both species over the IE if temperature were the leading influence. The negative correlation between SST and SNW of G. elongatus and G. bulloides at all three study sites during the IE, is indicative of an overriding effect of declining [CO 3 2− ] due to anthropogenic carbon input during the last 150 yr., masking any small opposing effect of temperature changes. Although temperature rise may have somewhat dampened the decrease in SNW, our results substantiate the dominant role of [CO 3 2− ] decline due to ocean acidification as the main driver of test calcite production.
Marine primary productivity in surface waters is assumed to influence test calcite production and therefore SNW 8,63 . Test calcification is energy demanding 64 , as greater food availability could allow foraminifera to allocate more energy to it. Chlorophyll a concentrations may serve as a proxy for marine primary productivity and food availability, as phytoplankton is considered to be the basis of planktic foraminiferal alimentation 65,66 . However, the relationship between marine primary productivity and SNW for our targeted species varies according to interspecific and intraspecific diversity and thus different calcification responses to changing trophic conditions. Previous studies detect weak positive/negative correlations between chlorophyll a and SNW for symbiont-barren G. bulloides/symbiont bearing G. elongatus in the Atlantic Ocean 8,58 , while a reverse correlation is observed for those species in the Mediterranean Sea 63 . Evidence for a minor effect of varying productivity on the calcification response of planktic foraminifera is given by the fact that, despite opposing productivity signals during the second half of the 20th century (Fig. 2b), we see consistent SNW decreases in both species at all three stations (Figs. 2e; 3f).
Changes of the chemical test microenvironment due to photosynthetic activity also influences calcification of planktic foraminifera 48 . Higher photosynthetic activity increases pH and reduces CO 2 concentrations near the shell, resulting in enhanced calcification rates 49 . If enhanced stratification led to enhanced light availability during the 20th century, this would likely have increased SNW of planktic foraminiferal species associated with symbionts capable of performing photosynthesis, though this effect is evidently masked by the dominant role of [CO 3 2− ] on calcification. Diagenetic processes can bias the significance of SNW as a proxy for planktic foraminiferal calcification, as dissolution/precipitation of carbonates on tests might be interpreted as reduced/enhanced mass of planktic foraminiferal shells. Across all three western Mediterranean records (Fig. 1), scanning electron microscope images ( Supplementary Figs. 5-7) of both species G. bulloides and G. elongatus show approximately the same degree of dissolution from top to bottom core depth with minor surface corrosion and breakup of the surface layer structure of the test wall 67,68 . Therefore, calcite saturation level of subsurface to bottom waters at the three core sites are assumed to not have significantly changed, and weight loss of settling tests in the deep water column as discussed by Schiebel et al. 69 may have been the same over the entire time interval discussed here. As western Mediterranean deep and bottom waters are supersaturated with regard to calcite 13,37,70 , dissolution at the sediment pore water interface is unlikely. The post-sapropel (<6 kyrs) western Mediterranean is a mesotrophic sea with limited changes and overall decreasing biological productivity over the past 150 years 16,71 . Accordingly, supra-lysoclinical dissolution as well as corrosive porewater conditions, which may infer post depositional test weight reduction in surface seafloor sediments through decomposing organic matter can be assumed absent 72 . This is in line with previous sediment trap and core top studies from the western Mediterranean Sea, reporting no signs of dissolution in fossil planktic foraminifera assemblages 32,33,73 . In contrast, early diagenetic calcite coating is reported for both species from core top samples at certain sites in the western and central Mediterranean Sea 74,75 . However, no signs of diagenetic overgrowth and secondary precipitated inorganic carbonates (e.g., inorganic calcite crusts, inorganic calcite crystals, overgrown pores) are present in the tests of G. elongatus and G. bulloides from the three core sites analyzed here ( Supplementary Figs. 5-7). Accordingly, due to the constant mild surface etching dissolution of similar magnitude of all tests observed across all three records, the relative foraminiferal SNW changes at the three core sites are expected to result from changes in biogenic calcification at the sea surface rather than diagenetic processes during sedimentation.
Driven by industrial carbon dioxide emissions, enhanced anthropogenic CO 2 penetration into the western Mediterranean Sea and resulting reductions in pH and [CO 3 2− ] are shown to be the key drivers of planktic foraminiferal SNW loss during the 20th century (Fig. 3a-f, Supplementary Fig. 2), corroborated by multiple regression analysis, discounting of multicollinearity between the independent parameters controlling SNW changes (see methods). Reduced SNW of similar magnitude during the 20th century have been observed in the tropical Atlantic, Southern Ocean, the California Current System, and the western Mediterranean Sea, showing similar amplitudes of relative calcification decrease 57,73,76,77 . Sea surface warming, productivity and solar activity changes are assumed to play only a negligible role, with rising SST and solar irradiance maxima potentially slightly attenuating the impact of anthropogenic CO 2 on planktic foraminiferal SNW. Accordingly, different oceanographic settings at the three selected core sites are responsible for the variability in 20th SNW trends (Supplementary Note 3), while the effects of diagenetic processes biasing the planktic foraminiferal SNW signal can be ruled out.

Conclusion
Accelerated anthropogenically-induced surface water pH decline in the western Mediterranean Sea is confirmed by boron isotopes and drives decreasing SNW of planktic foraminiferal tests. Negative carbon isotope excursions during the IE in tests of the planktic foraminifera species G. elongatus and G. bulloides display the signal of the Suess Effect-the anthropogenic CO 2 -driven shift to lower δ 13 C-and highlight the influence of anthropogenic emissions as the most important driver of foraminiferal calcite mass decline. Alkenone-derived SST and instrumental temperature records show unprecedented anthropogenic sea surface warming, that potentially mitigates this effect slightly. This basin-wide, interspecies multiproxy approach, covering multidecadal to centennial time scales, thus yields a much-needed high-confidence, low-variability assessment of the impacts of anthropogenic climate change on key marine species, in a critically sensitive region.

Methods
Material. The study is based on high resolution multicore records collected in the central and western Mediterranean Sea with a MC400-Multicorer system during the MedSeA cruise (Mediterranean Sea Acidification in a changing climate) on 2nd May to 2nd June 2013 onboard the R/V Àngeles Alvariño. Core sampling of MedSeA-S3-c1 and MedSeA-S23-c1 was previously described in ref. 16 . MedSeA-S7-c2 was retrieved in the Strait of Sicily (Lat. 37.7080°N, Long. 12.40553°E) at a water depth of 263 m, with a core length of 46.5 cm, sliced every centimeter.
Age Model development. Age models of cores MedSeA-S3-c1 and MedSeA-S23-c1 were estimated through a combination of radiocarbon and radionuclide dating, previously described in ref. 16 . The age model of MedSeA-S7-c2 (Supplementary  Table 2; Supplementary Note 4) was estimated through radionuclide analysis, determining total 210 Pb activity by measuring its alpha-emitter daughter nuclide 210 Po, following the methodology described in ref. 78 . 209 Po was added as an internal tracer, before sample aliquots of 200-300 mg were totally digested in acid media by using an analytical microwave oven and Po isotopes plated on silver discs in HCl 1 N at 70°C while stirring for 8 h. Po emissions were subsequently counted through α-spectrometers equipped with low background silicon surface barrier (SSB) detectors (EG&G Ortec) for 4 × 10 5 seconds. The difference between total 210 Pb and the constant 210 Pb at depth describes the concentration of excess 210 Pb, which was used to estimate maximum mean sediment accumulation rates and the age-depth model by applying the Constant Flux: Constant Sedimentation (CF:CS) model 79,80 . The activities of 137 Cs (661 keV) were determined by γ spectrometry in a coaxial high-purity Ge detector (EG&G Ortec) calibrated with the SRM-4276 solution standard supplied by the National Institute of Standards and Technology. The quality of the results determined by gamma and alpha spectrometry was evaluated by participation in IAEA proficiency tests and continuous analysis of certified and replicate materials.
Planktic foraminifera. Samples were dried at 60°C for approximately 24 hours to obtain dry bulk sediment mass. After washing over a 63 μm screen using distilled Elix water, the size fraction larger than 63 μm was oven dried at 60°C for 1 − 2 hours. Samples were dry sieved and divided by a riffle splitter (comparable to an ASC Sample Microsplitter, ASC Scientific) at the chosen size fractions: 150-250 μm, 250-315 μm, and > 315 μm. Total test accumulation rate was obtained by counting of more than 300 specimens in the >150 μm size fraction in all samples.
Size Normalized Weight. In order to obtain the calcite mass of planktic foraminiferal tests, a size related measure of test wall thickness and density-the Size Normalized Weight (SNW)-was obtained, which can be calculated by means of Eq. (1): Where C is the species average mass (μg) for each sample, A represents the average minimum Feret diameter (μm) for each sample and B is the average minimum Feret diameter (μm) of all specimens of one species of all samples. The equation (Eq. 1) includes the term to account for the relative SNW changes of the same species. Our methodology for measuring SNW is based on the sieve-based weight (SBW) technique, however, our approach also integrates measurements of individual test sizes in SNW estimations, a parameter used for the measurement-based weight (MBW) technique 82 . All specimens of a single sample are weighed together, while the average test size per sample is calculated by the minimum Feret diameter of each individual specimen. The SNW average standard deviation for each station (Alboran Sea: 7.0%; Balearic Sea: 8.3%; Strait of Sicily: 8.6%) is below the methodological 11% error reported by Beer, Schiebel 82 . Narrowing down the size fraction of analyzed tests to 250-315 µm ensured sufficient numbers of adult specimens which are large enough for SNW measurements. By applying the test weight sizenormalization technique described above, we were able to factor out size related mass changes, mainly influenced by temperature variations. SNW was determined for G. elongatus and G. bulloides. Bulk population measurements were performed on as many tests as available but not more than 10 and at least two specimens per sample per species, representing a total of 1,990 analyzed tests out of 224 samples included in this dataset. Selected tests were free of visible clay particles or organic matter and abnormally formed specimens were excluded from analysis. We weighed selected specimens together by triplicate with a Sartorius CP2P microbalance at Universitat Autònoma de Barcelona (UAB) and a Mettler Toledo XP6U Ultra Micro Balance at Max Planck Institute for Chemistry (MPIC) in Mainz, Germany, in environmentally controlled weighing rooms. Tests were positioned umbilical side up, photographed with a Canon EOS 650 D camera device attached to a Leica Z16 AP0 at UAB and with an Olympus UC90 attached to an Olympus SZY16 stereoscope at MPIC to measure the minimum Feret diameter of each individual test, using the software ImageJ 1.53k 83 and Olympus Stream Essentials 2.1. After measuring mean weight and mean minimum Feret diameter we calculated SNW for each sample using Eq. (1).
Stable isotope measurements. For chemical analysis, adult specimens were selected from a narrow size range (250-315 µm) to minimize uncertainty in paleoceanographic interpretation due to vital effects, following suggestions by Spero and Lea 89 . δ 13 C and δ 18 O. Bulk population measurements, including 5-10 specimens per species sample (not necessarily the same specimens as used for SNW measurements), were performed on G. elongatus and G. bulloides. Sample weight was determined through a Mettler Toledo XP6U Ultra Micro Balance at MPIC. Selected tests were crushed and soaked in ultrapure water (Milli-Q), before ultrasonication for up to 20 seconds, to remove clay contamination.
Samples were analyzed on a Thermo Delta V mass spectrometer equipped with a GASBENCH preparation device.~8-50 µg of CaCO 3 sample, placed in a Hefilled 4.5 ml exetainer vial was digested in water-free H 3 PO 4 at a temperature of 70°C. Subsequently the CO 2 -He gas mixture was transported to the GASBENCH in Helium carrier gas. In the GASBENCH, water vapor and various gaseous compounds were separated from the He-CO 2 mixture prior to sending it to the mass spectrometer in nine separate peaks. Isotope values of these individual peaks were averaged and reported as δ 13 C and δ 18 O relative to V-PDB. A total of 20 replicates of two in-house CaCO 3 standards were analyzed in each run of 55 samples. CaCO 3 standard weights were chosen so that they span the entire range of sample weights of the samples. After correction of isotope effects related to sample size the reproducibility of these standards is typically better than 0.1 ‰ (1 SD) for δ 18 O and 0.1 ‰ (1 SD) for δ 13 C.
Samples of 8 to 3 µg of CaCO 3 were analyzed on the same equipment using a cold trap technique that analyzes the entire sample in one single peak 90,91 . With this setting, 30 CaCO 3 standards were analyzed in each run of 30 samples. CaCO 3 standard weights were chosen so that they span the entire range of sample weights of the samples. Reproducibility of these <8 µg standards with the cold trap technique, after sample size correction, is typically better than 0.1 ‰ (1 SD) for δ 18 O and 0.1 ‰ (1 SD) for δ 13 C. δ 11 B. For boron isotope analysis, samples of G. ruber (250-315 µm) were prepared in a class 100 clean lab at the University of St Andrews STAiG laboratories. The two morphotypes of the species were analyzed separately. Samples were first crushed between two glass slides before transferring to acid cleaned centrifuge tubes in preparation for clay removal by ultrasonication. This was followed by oxidative cleaning in 1% hydrogen peroxide buffered with 0.1 M NH 4 OH and a weak acid leach in 0.0005 M distilled HNO 3 according to established protocols 92,93 . After addition of 100 µl boron-free Milli-Q, samples were dissolved by incremental addition of 0.5 M distilled HNO 3 (total volume of acid per sample 30-50 µl). A 5% cut of the dissolved sample was analyzed by ICP-MS for a suite of trace elements (Supplementary Table 3) including Mg/Ca on an Agilent 8900 following Foster 94 and Rae, Foster 93 with the addition of trace hydrofluoric acid in the wash to improve washout 95 . Samples were run alongside matrix-matched bracketing and consistency standards to monitor long-term instrument reproducibility which is reported at 95% confidence (2 SD). Low Al/Ca in all samples confirms the absence of residual clays. Elevated Mn/Ca was recorded in one sample but it does not display elevated Mg/Ca or any other anomalous trace element or isotopic signal.
After trace element analysis, a few samples with low boron concentrations, as determined from B/Ca, meant it was necessary to combine adjacent dissolved samples in MedSeA-S7-c2 core in order to reduce δ 11 B uncertainty (Supplementary Table 4). Boron was isolated from the sample matrix using boron-specific Amberlite IRA 743 resin 96,97 using a batch method. For this technique, 2 ml centrifuge tubes containing 50 µl of resin were pre-cleaned by addition and removal of 1 ml 0.5 M HNO 3 (repeated 3 times) and 1 ml boron-free Milli-Q (repeated twice) and then pre-conditioned with 100 µl 0.1 M ammonium acetate buffer at pH 8.5. Before loading, the sample was first buffered to pH 5 with 0.1 M ammonium acetate. After pipetting onto the resin, the tubes were agitated on a VWR mini vortex and the solution pipetted off after the resin settled. Following rinses with boron free Milli-Q water, the concentrated sample was eluted in 500 μl 0.5 M distilled HNO3 and the boron isotopic composition analyzed by MC-ICPMS with sample-standard bracketing against NIST 951 according to methods described by 93,94 . Addition of Hydrofluoric acid at 0.3 M helped to improve washout and reduce boron evaporation 95 . Matrix matched carbonate consistency standards were run alongside samples to monitor reproducibility which is ± 0.23 ‰ (2 SD).
pH Calculations CO2SYS. To estimate unknown variables in the seawater marine carbonate system, we used the CO2SYS program (version 2.0) 98 . Through a set of given input conditions (e.g., temperature, salinity, atm. CO 2 concentration), this program can estimate the remaining unknown variables of a two "master" carbonate system (TA, DIC, pH, pCO 2 ). We estimated surface ocean pH and carbonate ion concentrations for the three study sites, using seawater pCO 2 and total alkalinity (TA). TA was assumed constant over the last century, accounting for 2400 µmol/kgSW in the Alboran Sea 99 , 2478 µmol/kgSW in the Strait of Sicily 13 , and 2550 µmol/kgSW in the Balearic Sea 100 . Surface seawater pCO 2 was derived from temperature, salinity and pressure dependent solubility coefficient (K 0 ) following the formulation of Weiss 101 and atmospheric CO 2 measurements from an interpolated stack including the Mauna Loa time series 102 and reconstructed CO 2 concentrations 103,104 . Time-specific K 0 values were determined as input conditions for down core calculations using surface pressure, SST and salinity for each study location. SST was based on alkenone (U k' 37 ) derived paleo estimates at each study location 16,105 and salinity was defined as the constant mean (1950-2019 CE) of Hadley EN4 (version EN.4.2.2) subsurface salinity objective analysis 106 . Constants for our CO2SYS carbonate system estimations were the dissociation constant (K 1 and K 2 ) defined by 107 , refitted to seawater scale by 108 , the dissociation constants K HSO4 and K B according to 109 , and the chlorinity relationship and boron concentration based on 110 .
Boron isotope transformation. To calculate pH, the δ 11 B of G. ruber albus and G. elongatus was first converted to seawater δ 11 B borate using the calibration of 111 for the 250 − 300 µm size fraction: δ 11 B borate ¼ ðδ 11 B foram À9:52 ± 1:51Þ=ð0:60 ± 0:08Þ The boric acid dissociation constant is determined by both salinity and temperature and thus calculated pH can be susceptible to variations in these assigned values. Salinity was set at modern levels (37 PSU) for the region and temperature was calculated from Mg/Ca using the species specific calibration for G. ruber from 112 .
Mg=Ca ¼ 0:97 ± 0:50 x expð0:063 ± 0:019 * TÞ ð 3Þ The pH derived from Mg/Ca temperatures were compared to those from temperatures determined from Hadley, U k' 37 and δ 18 O SST in (Supplementary  Table 4). A sample depth of 25 m was given to approximate the upper habitation depth of the species below the sea surface 31 .

Statistical analysis
General additive model. Statistical analysis has been performed using the statistical software R 113 . General additive models (GAM) have been used in paleoenvironmental time series as a superior alternative tool to conventional statistical trend estimations 114 . We used a GAM to better show the basin wide, interspecies specific change in SNW and δ 13 C, and western Mediterranean SST changes. The complex, non-linear GAM fits as a smoothed function, which allows us to identify the period of significant SNW, δ 13 C and SST change based on proper accounting for model uncertainty. The general and applied mathematical form of the general additive model can be described as: SNW ¼ 4:081e À15 þ 6:997 þ ε i where ε i $ Nð0; σ 2 Þ ð 5Þ δ 13 C ¼ 0:02111 þ 6:238 þ ε i where ε i $ Nð0; σ 2 Þ ð6Þ SST ¼ À0:02993 þ 8:033 þ ε i where ε i $ Nð0; σ 2 Þ ð 7Þ SNW and δ 13 C values were standardized and merged, combining G. bulloides and G. elongatus at all three stations. Standardized SST is based on alkenone measurements and the calculated mean from HadISST dataset across the three stations. Uncertainty in the estimated model fits were addressed through 95% confidence intervals, after confirming gaussian error distribution via histogram plot, Q-Q plot and Kolmogorov-Smirnov test for standardized datasets of SNW (pvalue = 0.5501), δ 13 C (p-value = 0.1738) and SST (p-value = 0.7126). Additionally, a GAM fit was applied 115 , by means of the the mgcv package (version 1.8.42). To avoid under smoothing of the fit, the smoothing term was calculated using restricted maximum likelihood 116,117 .
Regression analysis. To quantify changes of SNW and δ 13 C during the 20th century, linear regression analysis was applied including all available data points post 1900 CE. Estimations of SNW and δ 13 C (negative anthropogenic carbon isotope excursions; nCIEs) change are defined by the difference between 1900 CE and 2013 CE, based on the regression model accounting for statistical significance (Supplementary Figs. 1, 2). To determine statistically significant relationships between the explanatory variable (SNW) and the independent variables (productivity, SST, pH and δ 13 C) we performed multiple linear regression analysis for time windows prior (pre-Industrial Era) and post 1800 CE (Industrial Era). Abundance changes of indicator species, along with changes in alkenone concentrations and planktic foraminiferal test accumulation rates (Fig. 2b) are included in the productivity variable. SST changes are based on the alkenone and HadISST datasets (Fig. 2a). The pH variability is calculated from pCO 2 and total alkalinity on CO2SYS (Fig. 2d), and δ 13 C changes are based on carbon isotope measurements (Fig. 2c). To generate evenly spaced time series, we standardized data for each station and parameter before combining station data for each parameter by interpolating nonlinear GAM fits (method see above) using generalized cross-validation (GCV) method for smoothing parameter estimation to avoid over smoothing of the fit. Goodness of the fit for each parameter is shown through deviance explained (Supplementary Table 5). Due to high collinearity, we have removed pH from the model including data post 1800 CE. Both parameters are highly correlated particularly during the Industrial Era as the δ 13 C signals mirrors the input of anthropogenic CO 2 , resulting in pH decrease (see discussion). The conditions of application required for the multiple linear regression analysis have been verified by using performance package (version 0.10.3). It can be assumed that there is no problematic amount of Multicollinearity since all parameters included in the model show a Variance Inflation Factor below ten 118 . Information on multiple linear regression models and their significance values can be found in the supplementary material (Supplementary Tables 6, 7), confirming a significant correlation between independent parameters and the explanatory variable, and substantiating the effect of ocean acidification on SNW.