Seven fossil mollusk specimens were collected by one of the authors (S.G.) for analysis from temporary exposures of the Oorderen Member of the Lillo Formation in the Antwerp area (Belgium; see Fig. 1). These were specimens of A. benedeni benedeni (specimen ID: SG126 and SG127) from Deurganckdoksluis (51°16′49″N, 4°14′56″E; collected in 2013), P. complanatus (SG115) and Glycymeris radiolyrata (SG116) from Deurganckdok (51°17′24″N, 4°15′37″E; collected on 2 February 2001), and Pygocardia rustica (SG105), Ostrea edulis (SG113), and Arctica islandica (SG107) from Verrebroekdok (51°16′16″N, 4°12′53″E; collected in 1999–2000). At these localities, the Oorderen Member, from bottom to top, is divided into basal shell bed, Atrina level, Cultellus level (SG115 and SG116), and benedeni level (SG126 and SG127) (37). Specimens SG105, SG113, and SG107 were collected ex situ within the Oorderen Member. The estimated age range for the Oorderen Member is 2.72 to 3.3 Ma (see Fig. 1) (35, 67).

Field observations and a detailed assessment of the composition of the invertebrate fauna reveal that the fossil-bearing sediments of the Oorderen Member were deposited between 30 and 50 m water depth during warm, highstand intervals (35, 36). Previous studies have argued that the SNS undergoes summer stratification during the mPWP, which would cause molluscan shell reconstructions to underestimate summer temperatures (24, 47). However, the lack of dysoxic faunas in the Oorderen Member suggests the absence of persistent stratification (35). Furthermore, the presence of sedimentary structures indicative of tidal currents and strong vertical mixing during deposition of the Oorderen Member argues against a large surface to seafloor temperature gradient (39). Therefore, it is likely that the mollusks recorded temperatures close to the SSTs year round. This assessment is supported by the evidence of warm spring and cool autumn temperatures, suggesting a direct response of the water temperature experienced by the mollusks to the seasonal cycle in SATs (see Results; Fig. 3). However, we cannot fully exclude that spatially and temporally restricted temperature stratification did occur during the lifetime of the studied specimens (24, 47). Note that, even if summer stratification occurred, the summer temperatures recorded by the mollusks would underestimate the true mPWP summer SST and the occurrence of enhanced summer warming during the mPWP would still be supported by the data.

Shells were partially embedded in epoxy resin before polished thick sections were prepared to expose a cross section through the axis of maximum shell growth [following (68); see fig. S3]. Preservation of the original shell calcite and aragonite was verified using a combination of scanning electron microscopy (SEM), cathodoluminescence microscopy, electron backscatter diffraction microscopy (EBSD), micro–x-ray fluorescence (μXRF), and x-ray diffraction [XRD; see (67) and figs. S4 to S12]. Only shells with excellent preservation as demonstrated by these methods were considered for clumped isotope analysis.

Carbonate was sampled along transects in the direction of growth in transects on cross sections through the outer shell layers of all shells except P. complanatus (SG115) and A. benedeni benedeni (SG126 and SG127), whose thin outer shell layers necessitated sampling on the outside of the shell. A combination of handheld drilling (outside) and micromilling (in transects) was used to sample along narrow growth increments to minimize time averaging (Supplementary Materials). The number of replicates per specimen ranges from 61 to 157 (Fig. 2): SG105, n = 98; SG107, n = 63; SG113, n = 90; SG115, n = 100; SG116, n = 157; SG126, n = 79; SG127, n = 61. Some samples were replicated more than once, causing multiple replicates to originate from the same location in the shell, and some of these replicate measurements were discarded as outliers during data processing (see below). Measurements from A. benedeni benedeni (SG126 and SG127) include data from (67) as well as additional analyses on the same specimens carried out for this study. Small (70 to 160 μg) aliquots of each sample were digested in phosphoric acid for 600 s at 70°C in a Kiel IV carbonate preparation device, after which the resultant CO₂ gas was purified using two cold fingers and a Porapak trap (69) before the clumped isotope composition (Δ₄₇) was analyzed using two Thermo MAT253 mass spectrometers. A long-integration dual-inlet (LIDI) workflow (69, 70) was used for all except 53 aliquots from P. rustica (SG105); the latter aliquots were measured using sample-standard measurement cycles [“click-clack” mode; e.g., (71)]. Clumped isotopic ratios were corrected for intensity-based Δ₄₇ offsets based on intensity-matched ETH-3 standards (see section S3.1) before being standardized using the ETH standards (72), which were run in approximately 1:1 ratio with the samples (73) and corrected for ¹⁷O concentration following (74). Long-term analytical precision was monitored using IAEA-C2 and Merck reference materials (typical standard deviation of 0.04‰; see table S2) after outlier removal based on metadata on instrument performance (see criteria in section S3.2).

For each mollusk specimen, samples were internally dated relative to the seasonal cycle using ShellChron (75), after which measurements within a specimen were grouped in four 3-monthly bins (“seasons”). Summer was defined as the consecutive 3-month period with the lowest mean Δ₄₇ value, calculated separately for each specimen (see section S3.2 and figs. S20 to S32). The fact that these 3-monthly bins do not always line up with the minima and maxima in the stable oxygen isotope profiles through the shells shows that the common assumption that these extreme stable oxygen isotope values represent summer and winter seasons does not always hold true [see (76)]. Mean seasonal Δ₄₇ values were obtained through weighted averages of seasonal means and uncertainties considering uncertainty within and variability between specimens. Seasonal temperatures during the mPWP were reconstructed from Δ₄₇ for each seasonal group using the clumped isotope calibration by (77), propagating uncertainty on Δ₄₇ values following the procedure described in the supplement of (78).

Seasonal SST from the SNS was extracted from the corresponding local ocean grid cells (51°N to 55°N, 2°E to 4°E) in the PlioMIP2 model ensemble (see Fig. 2). SNS SST records were obtained from the same area within the ERSSTv5 product (41) and supplemented with local SST observations from the NIOZ monitoring station (42). Future projections matching the SSP2-4.5 scenario (48) were produced for CMIP6, and the regional SST outcomes were exported from the Intergovernmental Panel on Climate Change (IPCC) WG1 Interactive Atlas (79). Seasons were defined as weighted means of the 3 months containing the largest number of days in the season according to the astronomical definition: winter: January to March; spring: April to June; summer: July to September; autumn: October to December. Data processing was carried out using the open-source computational software R (80), and scripts are provided in (81).