Freshwater pearl mussels from northern Sweden serve as long-term, high-resolution stream water isotope recorders

The stable isotope composition of lacustrine sediments is routinely used to infer Late Holocene changes in 15 precipitation over Scandinavia and ultimately, atmospheric circulation dynamics in the North Atlantic realm. However, such archives provide only low temporal resolution (ca. 15 years) precluding the ability to identify changes on inter-annual and quasi-decadal time-scales. Here we present a new, high-resolution reconstruction using shells of freshwater pearl mussels, Margaritifera margaritifera, from three rivers in north Sweden. We present seasonally to annually resolved, calendar-aligned stable oxygen and carbon isotope data from ten specimens covering the time interval from 1819 to 1998. The studied bivalves 20 formed their shells near equilibrium with the oxygen isotope signature of ambient water and thus reflected hydrological processes in the catchment as well as changes, albeit damped, in the isotope signature of local atmospheric precipitation. The shell oxygen isotopes were correlated significantly with the North Atlantic Oscillation index (up to 56 % explained variability), suggesting that moisture from which winter precipitation formed originated predominantly in the North Atlantic during NAO+ years, but in the Arctic during NAO– years. The isotope signature of winter precipitation was attenuated in the stream water 25 and this damping effect was eventually recorded by the shells. Shell stable carbon isotope values did not show consistent ontogenetic trends, but rather oscillated around an average that ranged from ca. -12.00 to -13.00 ‰ among the studied rivers. Results of this study contribute to an improved understanding of climate dynamics in Scandinavia and the North Atlantic sector and can help to constrain eco-hydrological changes in riverine ecosystems. Moreover, long isotope records of precipitation and streamflow are pivotal for improving our understanding and modeling of hydrological, ecological, biogeochemical and 30 atmospheric processes. Our new approach offers a much higher temporal resolution and superior dating control than data from existing archives.


Introduction
Multi-decadal records of δ 18 O signals in precipitation and stream water are important for documenting climate change impacts on river systems (Rank et al., 2017), improving mechanistic understanding of water flow and quality controlling processes 35 (Darling and Bowes, 2016) and testing earth hydrological and land surface models (Reckerth et al., 2017;Risi et al., 2016; see also Tetzlaff et al., 2014). However, the common sedimentary archives used for such purposes typically do not provide the required, i.e., at least annual temporal resolution (e.g., Rosqvist et al., 2007). In such studies, the temporal changes in the oxygen isotope signatures of meteoric water are encoded in biogenic tissues and abiogenic minerals formed in rivers and lakes (e.g., Teranes et al. 2001;Leng & Marshall 2004). In fact, many studies have determined the oxygen isotope composition of 40 diatoms, ostracods, authigenic carbonate and aquatic cellulose preserved in lacustrine sediments to reconstruct Late Holocene changes in precipitation over Scandinavia and ultimately, atmospheric circulation dynamics in the North Atlantic realm (e.g., Hammarlund et al., 2002;Andersson et al., 2010;Rosqvist et al., 2004Rosqvist et al., , 2013. With their short residence times of a few months (Rosqvist et al., 2013), the hydrologically connected, through-flow lakes of northern Scandinavia are an ideal region for this type of study. Their isotope signatures -while damped in comparison to isotope signals in precipitation -directly respond to 45 changes in precipitation isotope composition (Leng & Marshall 2004). However, even the highest available temporal resolution of such records (15 years per sample in sediments from Lake Tibetanus, Swedish Lappland; Rosqvist et al., 2007) is still insufficient to resolve inter-annual to decadal-scale variability, i.e., the time-scales at which the North Atlantic Oscillation (NAO) operates (Hurrell, 1995;Trouet et al., 2009). The NAO steers weather and climate dynamics in northern Scandinavia and determines the origin of air masses from which meteoric waters form. While sediment records are still of vital importance 50 for century and millennial-scale variations, new approaches are needed for finer scale resolution on the 1-100-year time-scale.
One underutilized approach for hydroclimate reconstruction is the use of freshwater mussels as natural stream water stable isotope recorders (Dettman et al., 1999;Kelemen et al., 2017;Pfister et al., 2018Pfister et al., , 2019. Their shells can provide seasonally to annually resolved, chronologically precisely constrained records of environmental changes in the form of variable increment 55 widths (= distance between subsequent growth lines) and geochemical properties (e.g., Nyström et al. 1996;Schöne et al. 2005a;Geist et al., 2005;Black et al., 2010;Schöne and Krause, 2016;Geeza et al., in press a, b;Kelemen et al., in press). In particular, similar to marine (Epstein et al., 1953;Mook & Vogel, 1968;Killingley & Berger, 1979) and other freshwater bivalves (Dettman et al., 1999;Kaandorp et al., 2003;Versteegh et al., 2009;Kelemen et al., 2017;Pfister et al., 2019), M. margaritifera forms its shell near equilibrium with the oxygen isotope composition of the ambient water (δ 18 O w ) (Pfister et al., 60 2018;Schöne et al., 2005a). If the fractionation of oxygen isotopes between the water and shell carbonate is only temperaturedependent and the temperature during shell formation is known or can be otherwise estimated (e.g., from shell growth rate), reconstruction can be made of the oxygen isotope signature of the water from that of the shell CaCO 3 (δ 18 O s ). Freshwater pearl mussels, Margaritifera margaritifera, are particularly useful in this respect, because they can reach a lifespan of over 200 years (Ziuganov et al., 2000;Mutvei & Westermark, 2001) offering an insight into long-term changes of freshwater ecosystems 65 at unprecedented temporal resolution.
Here we present the first, absolutely dated, annually resolved stable oxygen and carbon isotope record of freshwater pearl mussels from three different rivers in northern Sweden covering nearly two centuries . We test the ability of these freshwater pearl mussels to reconstruct the NAO index and associated changes in precipitation provenance using shell oxygen 70 isotope data (δ 18 O s ). We evaluate how shell oxygen isotope data compare to δ 18 O values in stream water (δ 18 O w ) and precipitation (δ 18 O p ) as well as limited existing environmental data (such as stream water temperature) and how these variables relate to each other. We hypothesize that δ 18 O w and δ 18 O p values are positively correlated to each other as well as to the North Atlantic Oscillation index. In other terms, during positive NAO years, oxygen isotope values in stream water and shells are higher than during negative NAO years, when shell and stream water oxygen isotope values tend to be more depleted in 18 O. 75 In addition, we explore the physical controls on shell stable carbon isotope signatures. Presumably, these data reflect changes in the stable carbon isotope value of dissolved inorganic carbon, which in turn, is sensitive to changes in primary production.
We leverage past work in Sweden that has shown that the main growing season of M. margaritifera occurs from mid-May to mid-October, with fastest growth rates occurring between June and August (Dunca & Mutvei, 2001;Dunca et al. 2005;Schöne et al., 2004aSchöne et al., , b, 2005a. We use changes in annual increment width of M. margaritifera to infer water temperature, because 80 growth rates are faster during warm summers and result in broader increment widths (Schöne et al., 2004a(Schöne et al., , b, 2005a. Since specimens of a given population react similarly to changes in temperature, their average shell growth patterns can be used to estimate climate and hydrological changes. Consequently, increment series of specimens with overlapping lifespans can be crossdated and combined to form longer chronologies covering several centuries (Schöne et al., 2004a(Schöne et al., , b, 2005a.

2 Material and methods
We collected ten specimens of the freshwater pearl mussel, M. margaritifera from three rivers in Norrland (Norrbottens län), northern Sweden ( Fig. 1+2; Table 1). Bivalves were collected between 1993 and 1999 and included nine living specimens and one found dead and articulated (bi-valved; ED-GJ-D6R). Four individuals were taken from the river Nuortejaurbäcken (NJB), two from Grundträsktjärnbäcken (GTB) and four from Görjeån (GJ) (Fig. 1). Since M. margaritifera is an endangered species 90 (Moorkens et al., 2018), we refrained from collecting additional specimens that could have covered the time interval since collection and manuscript preparation, but relied on bivalves that we obtained -with permission -for a postdoctoral project of one co-author (AEM, formerly known as Elena Dunca) and the doctoral thesis of another co-author (SMB).
The bedrock in the studied catchments is dominated by orthogneiss and granodiorite. The vegetation at GTB (ca. 90 m a.s.l. 95 [above sea-level]) and GJ (ca. 200 m a.s.l.) consisted of a mixed birch forest, while conifers, shrubs and bushes dominated at NJB (ca. 400 m a.s.l.). The studied rivers were thus rich in humin acids. The studied streams were fed by small upstream open (flow-through) lakes.

Sample preparation
The soft tissues were removed immediately after collection and shells were then air-dried. One valve of each specimen was 100 wrapped in a protective layer of WIKO metal epoxy resin #5 and mounted to a Plexiglas cube using Gluetec Multipower plastic welder #3. Shells were then cut perpendicular to the growth lines using a low-speed saw (Buehler Isomet) equipped with a diamond-coated (low-diamond concentration) wafering thin blade (400 µm thickness]. One specimen (ED-NJB-A3R) was cut along the longest axis, all others along the height axis from the umbo to the ventral margin (Fig. 2a). From each specimen, two ca. 3 mm-thick shell slabs were obtained and mounted onto glass slides with the mirroring sides (= the portions that were 105 located to the left and right of the saw blade during the cutting process) facing upward. This method facilitated the temporal alignment of isotope data measured in one slab to growth patterns determined in the other shell slab. The shell slabs were ground on glass plates using suspensions of 800 and 1200 grit SiC powder and subsequently polished with Al 2 O 3 powder (grain size: 1 µm) on a Buehler G-cloth. Between each grinding step and after polishing, the shell slabs were ultrasonically cleaned with water. 110

Shell growth pattern analysis
For growth pattern analysis, one polished shell slab was immersed in Mutvei's solution for 20 min at 37-40° C under constant stirring (Schöne et al., 2005b). After careful rinsing in demineralized water, the stained sections were air-dried under a fume hood. Dyed thick-sections were then viewed under a binocular microscope (Olympus SZX16) equipped with sectoral dark field illumination (Schott VisiLED MC1000) and photographed with a Canon EOS 600D camera (Fig. 2b). The widths of the 115 annual increments were determined to the nearest ca. 1 µm with the image processing software (Panopea, © Peinl and Schöne).
Measurements were completed in the outer portion of the outer shell layer (oOSL, prismatic ultrastructure) from the boundary between the oOSL and the inner portion of the outer shell layer (iOSL, nacreous ultrastructure) perpendicularly to the previous annual growth line (Fig. 2c). Annual increment width chronologies were detrended with stiff cubic spline functions and standardized to produce dimensionless measures of growth (standardized growth indices, i.e., SGI values [σ]) following 120 standard sclerochronological methods (Helama et al., 2006;Butler et al., 3013;Schöne, 2013). Briefly, for detrending, measured annual increment widths were divided by the data predicted by the cubic spline fit. From each resulting growth index, we subtracted the mean of all growth indices and divided the result by the standard deviation of all growth indices of the respective bivalve specimen. This transformation resulted in SGI chronologies. Due to low heteroscedasticity, no variance correction was needed (Frank et al., 2007). Uncertainties in annual increment measurements resulted in an SGI error of ± 0.06 125 σ.

Stable isotope analysis
The other polished shell slab of each specimen was used for stable isotope analysis. To avoid contamination of the shell aragonite powders (Schöne et al., 2016), the cured epoxy resin and the periostracum were completely removed prior to sampling. A total of 1,551 powder samples (32-128 µg) were obtained from the oOSL by means of micromilling (Fig. 2c) 130 under a stereomicroscope at 160 × magnification. An equidistant sampling strategy was applied, i.e., within each annual increment the milling step size was held constant (Schöne et al., 2005c). We used a cylindrical, diamond-coated drill bit (1 mm diameter; Komet/Gebr. Brasseler GmbH and Co. KG, model no. 835 104 010) mounted on a Rexim Minimo drill. While the drilling device was affixed to the microscope, the sample was handheld during sampling. In early ontogenetic years, up to 16 samples were obtained between successive annual growth lines. In the latest ontogenetic portions of specimens ED-NJB-135 A2R (last year of life) and ED-GJ-D6R (last nine years of life), each isotope sample represented two to three years.
Stable carbon and oxygen isotopes were measured at the Institute of Geosciences at the University of Frankfurt/Main (Germany). Carbonate powder samples were digested in He-flushed borosilicate exetainers at 72 °C using a water-free phosphoric acid. The released CO 2 gas was then measured in continuous flow mode with a ThermoFisher MAT 253 gas source 140 isotope ratio mass spectrometer coupled to a GasBench II. Stable isotope ratios were corrected against an NBS-19 calibrated Carrara marble (δ 13 C = +2.02 ‰; δ 18 O = -1.76 ‰). Results are expressed as parts per thousand (‰) relative to the Vienna Pee Dee Belemnite (V-PDB). The long-term accuracy based on blindly measured reference materials with known isotopic composition is better than 0.05 ‰ for both isotope systems. Note that no correction was applied for differences in fractionation factors of the reference material (calcite) and shells (aragonite; verified by Raman spectroscopy), because the 145 paleothermometry equation used below (Eq. 2) did also not consider these differences (Füllenbach et al., 2015). However, the correction of -0.38 ‰ would be required if δ 18 O values of shells and other carbonates were compared with each other.

Instrumental data sets
Shell growth and isotope data were compared to a set of environmental variables including the station-based winter (DJFM) NAO index (obtained from https://climatedataguide.ucar.edu, last access: 9 April 2019) as well as oxygen isotope values of 150 river water (δ 18 O w ) and weighted (corrected for precipitation amounts) oxygen isotope values of precipitation (δ 18 O p ). Data on monthly river water and precipitation came from the Global Network of Isotopes in Precipitation (GNIP) and Rivers (GNIR), available at the International Atomic Energy Agency from their website at https://nucleus.iaea.org/wiser/index.aspx (last access: 1 April 2019). Furthermore, monthly air temperature (T a ) data came from the station Stensele, available at the Swedish Meteorological and Hydrological Institute (https://www.smhi.se). From these data, monthly river water temperature (T w ) was 155 computed by using the summer air-river water temperature conversion by Schöne et al. (2004a)  ( 1 )

Weighted annual shell isotope data 160
Since shell growth rate varied during the growing season with fastest biomineralization rates occurring during June and July (Dunca et al., 2005), annual growth increments are biased toward summer, and powder samples taken from the shells at equidistant intervals represent different amounts of time. To compute growing season averages (henceforth referred to as 'annual averages') from such intra-annual shell isotope data (δ 18 O s , δ 13 C s ), weighted (henceforth denoted with an asterisk) annual means are thus needed, i.e., δ 18 O s * and δ 13 C s * values (Schöne et al., 2004a). The relative proportion of time of the 165 growing season represented by each isotope sample was computed from a previously published intra-annual growth curve of juvenile M. margaritifera from Sweden (Dunca et al., 2005). For example, if four isotope samples were taken between two winter lines at equidistant intervals, the first sample would represent 22.38 % of the time of the main growing season duration, the second, third and forth 20.28, 24.47 and 32.87 %, respectively (Table 2). Accordingly, the weighted annual mean isotope values (δ 18 O s *, δ 13 C s *) were calculated by multiplying these numbers (= weights) with the respective δ 18 O s and δ 13 C s values 170 and dividing the sum of the products by 100. The four isotope samples of the example above comprise the time intervals of 23 May-22 June, 23 June-21 July, 22 July-25 August, and 26 August-12 October, respectively. Missing isotope data due to lost powder, machine error, air in the exetainer etc. were filled in by linear interpolation in 20 instances. We assumed that the timing and rate of seasonal growth remained nearly unchanged throughout lifetime of the specimens and in the study region (see also Sect. 4). 175

Reconstruction of oxygen isotope signatures of river water on annual and intra-annual time-scales
To assess how well the shells recorded δ 18 O w values on inter-annual time-scales, the stable oxygen isotope signature of river water (δ 18 O wr *) during the main growing season ('annual' δ 18 O wr *) was reconstructed from δ 18 O s * data and the arithmetic average of (monthly) river water temperatures, T w , during the same time interval, i.e., 23 May-12 Oct. Through this approach, the effect of temperature-dependent oxygen isotope fractionation was removed from the δ 18 O s * data. For this purpose, the 180 paleothermometry equation of Grossman and Ku (1986; corrected for the VPDB-VSMOW scale difference following Gonfiantini et al., 1995)  .

185
Since air temperature data were only available from 1860 onward, T w prior to that time were inferred from age-detrended and standardized annual growth increment data (SGI values) using a linear regression model similar to that introduced by Schöne et al. (2004a). In the revised model, SGI data of 25 shells from northern Sweden (fifteen published chronologies, provided in the article cited above, and ten new chronologies from the specimens studied in the present paper) were arithmetically averaged for each year and then regressed against weighted annual water temperature, henceforth referred to as annual T w *. The annual 190 T w * data consider variations in seasonal shell growth rate. 6.29, 25.49, 24.52, 21.92, 16.88 and 4.90 % of the annual growth increment were formed in each month between May and October, respectively. The values were multiplied with T w of the corresponding month and the sum of the products divided by 100 to obtain annual T w * data. The revised (shell growth vs. For coherency purposes, we also applied this model to post-1859 SGI values and computed river water temperatures that were subsequently used to estimate δ 18 O wr(SGI) * values.

200
To assess how well the shells recorded δ 18 O w values at intra-annual time-scales, we focused on two shells from NJB (ED-NJB-A4R and ED-NJB-A6R), which provided the highest isotope resolution of one to two weeks per sample during the few years of overlap with GNIP and GNIR data. Note, only for this bivalve sampling locality, monthly instrumental oxygen isotope data were available from the GNIP and GNIR data sets (data by Burgman et al., 1981 temperatures were computed as weighted averages, T w *, from monthly T w considering seasonal changes of the shell growth rate. For example, if four powder samples were taken from the shell at equidistant intervals within one annual increment, 6.29 % of the first sample were formed in May and 18.63 % in June (sum ca. 25 %). The average temperature during that time interval is computed with these numbers as follows: (T w of May × 0.0629 + T w of June × 0.1863) / 25. 6.86 % of the second sample from that annual increment formed in June and 17.97 % in July. Accordingly, the average temperature was: (T w of 215 June × 0.0686 + T w of July × 17.47) / 25 etc. Note that annual δ 18 O wr * values can also be computed from intra-annual δ 18 O wr * data, but this approach is much more time-consuming and complex than the method described further above. Yet, both methods produce nearly identical results.

Stable carbon isotopes of the shells
Besides the winter and summer NAO index, weighted annual stable carbon isotope data of the shells, δ 13 C s * values, were 220 compared to shell growth data (SGI chronologies). Since the δ 13 C s * values could potentially be influenced by ontogenetic effects, the chronologies were detrended and standardized (δ 13 C s(d) *) following methods typically used to remove ontogenetic age trends from annual increment width chronologies (see e.g., Schöne, 2016). Detrending was done with cubic spline functions capable of removing any directed trend toward higher or lower values through lifetime.

Results 225
The lengths of the annual increment chronologies of M. margaritifera from the three studied rivers (Nuortejaurbäcken, Grundträsktjärnbäcken and Görjeån) ranged from 21 to 181 years and covered the time interval of 1819 to 1999 AD (Table 1).
Since the umbonal shell portions were deeply corroded and the outer shell layer was missing -a typical feature of long-lived freshwater bivalves (Schöne et al., 2004a ; Fig 2a) -the actual ontogenetic ages of the specimens could not be determined and may have been up to 10 years higher than listed in Table 1. 230

Shell growth and temperature
The ten new SGI series from NJB, GTB and GJ were combined with fifteen published annual increment series of M. margaritifera from the rivers Pärlälven, Pärlskalsbäcken and Bölsmanån (Schöne et al., 2004a(Schöne et al., , b, 2005 to a revised Norrland master chronology. During the 50-year calibration interval of 1926-1975 (the same time interval was used in the previous study by Schöne et al., 2004Schöne et al., a, b, 2005, the chronology was significantly (p < 0.05; note, all p-values of linear regression 235 analyses in this paper are Bonferroni-adjusted) and positively correlated (R = 0.74; R 2 = 0.55) to the weighted annual river water temperature (T w *) during the main growing season (Fig. 3). These values were similar to the previously published coefficient of determination for a stacked record using M. margaritifera specimens from rivers across Sweden (R 2 = 0.60 (Schöne et al., 2005); note that this number is for SGI vs. an arithmetic annual T w ; a regression of SGI against weighted annual T w returns an R 2 of 0.64). 240

Shell stable oxygen isotope data
The shell oxygen isotope curves showed distinct seasonal and inter-annual variations (Fig. 4+5). The former were particularly well developed in specimens from GTB and NJB (Fig. 4), which were sampled with a very high spatial resolution of ca. 30 µm (ED-GTB-A1R, ED-GTB-A2R, ED-NJB-A4R, ED-NJB-A6R). In these shells, up to 16 samples were obtained from single annual increments translating into a temporal resolution of one to two weeks per sample. Typically, the highest δ 18 O s values 245 of each cycle occurred at the winter lines and lowest values about half way between consecutive winter lines (Fig. 4). The largest seasonal δ 18 O s amplitudes of ca. 2.20 ‰ were measured in specimens from GTB (-8.68 to -10.91 ‰) and ca. 1.70 ‰ in shells from NJB (-8.63 to -10.31 ‰).
Weighted annual shell oxygen isotope (δ 18 O s *) values fluctuated on decadal time-scales (common period: ca. 8 years) with 250 amplitudes larger than those occurring on seasonal scales, i.e., ca. 2.50 and 3.00 ‰ in shells from NJB (-8.63 to -11.10 ‰) and GTB (-7.84 to -10.85 ‰), respectively (Fig. 5a+b). The chronologies from GJ also revealed a century-scale variation with minima in the 1820s and 1960s and maxima in the 1880s and 1990s (Fig. 5c). The δ 18 O s * curves of specimens from the same locality showed notable agreement in terms of absolute values and visual agreement (running similarity), specifically specimens from NJB and GTB (Fig. 5a+b). However, the longest chronology from GJ showed only little agreement with the 255 remaining three series from that site (Fig 5c). The similarity among the series also changed through time (Fig 5a-c). In some years, the difference between the series was less than 0.20 ‰ at NJB (N = 4) and GTB (N = 2) (1983) and 0.10 ‰ at GJ (1953;N = 4), whereas in other years, differences varied by up to 0.82 ‰ at NJB and 1.00 ‰ at GTB and GJ, respectively. Average shell oxygen isotope chronologies of the three studied rivers exhibited a strong running similarity (passed the 'Gleichläufigkeitstest' by Baillie and Pilcher, 1983, for p < 0.001) and were significantly positively correlated to each other 260 (R 2 values of NJB vs GTB = 0.34, NJB vs GJ = 0.40 and GTB vs GJ = 0.36; all at p < 0.0001).

Shell stable oxygen isotope data and instrumental records
At NJB -the only bivalve sampling site for which measured stream water isotope data were available from nearby localities 1977 followed by a slight increase of ca. 0.50 ‰ until 1980. In contrast to the damped stream water signal (average seasonal range during four years for which both stream water and precipitation data were available [1975,1976,1978,1979] Despite the limited number of instrumental data, seasonally averaged δ 18 O wr * data showed some -though not always statistically significant -agreement with δ 18 O w and weighted δ 18 O p data (corrected for precipitation amounts), respectively, both in terms of correlation coefficients and absolute values (Table 3). These findings were corroborated by the regression analyses of instrumental δ 18 O p values against δ 18 O w values (Table 3). For example, the oxygen isotope values of summer 280 (June-Sep) precipitation were significantly (Bonferroni-adjusted p < 0.05) and positively correlated to those of shell carbonate precipitated during the same time interval (98 % explained variability in both specimens, but only at p < 0.05 in ED-NJB-A6R). Likewise, δ 18 O w and δ 18 O p values during summer were positively correlated to each other (R = 0.91), though less significant (p = 0.546). Strong relationships were also found for δ 18 O wr * and δ 18 O w values during the main growing season as well as annual δ 18 O wr * and Dec-Sep δ 18 O p values. The underlying assumption for the latter was that the δ 18 O wr * average value 285 reflects the combined δ 18 O p of snow precipitated during the last winter (received as meltwater during spring) and rain precipitated during summer. Instrumental data supported this hypothesis, because stream water δ 18 O values during the main growing season were highly significantly and positively correlated to Dec-Sep δ 18 O p data (Table 3). On the other hand, changes of the isotope signal of winter (Dec-Feb) snow were only weakly and not significantly mirrored by changes in river water oxygen isotope values during the snowmelt period (May-June), or in δ 18 O wr * values from shell portions formed during the 290 same time interval (Table 3). During the four years under study (1975,1976,1978,1979), measured and reconstructed δ 18 O w values were nearly identical during the main growing season (δ 18 O w = -12.46 ‰; δ 18 O wr * = -12.57 ‰ and -12.46 ‰) and during summer (δ 18 O w = -12.39 ‰; δ 18 O wr * = -12.44 ‰ and -12.43 ‰) (Table 3). In contrast, isotopes in precipitation and river water showed larger discrepancies (see above and Fig. 6b, Table 3).

Shell stable oxygen isotope data and synoptic circulation patterns (NAO) 295
Site-specific annual δ 18 O wr * (and δ 18 O wr(SGI) *) chronologies (computed as arithmetic averages of all chronologies at a given river) were significantly (Bonferroni-adjusted p < 0.05) positively correlated to NAO indices ( Fig. 7; Table 4). In NAO + years, the δ 18 O wr * (and δ 18 O wr(SGI) *) values were higher than during NAOyears. The strongest correlation existed between the winter (Dec-Mar) NAO and δ 18 O wr * (and δ 18 O wr(SGI) *) values at NJB (44 to 49 % explained variability). At GTB, the explained variability ranged between 24 and 27 %, whereas at GJ only 16 to 18 % of the inter-annual δ 18 O wr * (and δ 18 O wr(SGI) *) variability 300 was explained by the winter NAO (wNAO) index. Between 1947 and 1991, the time interval for which isotope data were available for all sites, R 2 values were more similar to each other and ranged between 0.27 and 0.46 (Table 4). All sites reflected well-known features of the instrumental NAO index series such as the recent  positive shift toward a more dominant wNAO, which delivered isotopically more positive (less depleted in 18 O) winter precipitation to our region of interest ( Fig. 7a-c). The correlation between δ 18 O wr * (and δ 18 O wr(SGI) *) values and the summer (June-Aug) NAO index was much lower 305 than for the wNAO, but likewise positive and sometimes significant at p < 0.05 (Table 4). Between 1947 and1991, 7 to 43 % of the inter-annual oxygen isotope variability was explained by the summer NAO index.
We have also computed an average δ 18 O wr(SGI) * curve for the entire study region (Fig. 8a-c
The detrended and standardized annual shell stable carbon isotope (δ 13 C s(d) ) curves showed no statistically significant 335 (Bonferroni-adjusted p < 0.05) agreement with the NAO indices or shell growth rate (SGI values) ( Fig. 7; Table 4). A weak negative correlation (10 % explained variability) only existed between δ 13 C s(d) * values and the wNAO at NJB. Some visual agreement was apparent between δ 13 C s(d) values and SGI in the low frequency realm. For example, at NJB, faster growth during the mid-1950s, 1970s, 1980s and 1990s fell together with lower δ 13 C s(d) values (Fig. 7g). Likewise at GTB, faster shell growth seemed to be inversely linked to δ 13 C s(d) values (Fig. 7h). 340

Advantages and disadvantages of using bivalve shells for stream water δ 18 O reconstruction; comparison to sedimentary archives
Our results have shown that shells of freshwater pearl mussels from rivers in northern Scandinavia (fed predominantly by small, open lakes and precipitation) can serve as a long-term, high-resolution archive of the stable oxygen isotope signature of 345 the water in which they lived. Since δ 18 O w values have a much lower seasonal amplitude than δ 18 O p values (i.e., δ 18 O w signals are damped relative to δ 18 O p data as a result of the water transit times through the catchment of the stream), the observed and reconstructed stream water isotope signal mirror the seasonal and inter-annual variability in δ 18 O p values. The NAO and subsequent atmospheric circulation patterns determine the origin of air masses and subsequently the δ 18 O signal in precipitation. 350 Compared to lake sediments, which have traditionally been used for similar reconstructions at nearby localities (e.g., Hammarlund et al., 2002;Andersson et al., 2010;Rosqvist et al., 2004Rosqvist et al., , 2013, this new shell-based archive has a number of advantages.

355
(1) The effect of temperature-dependent oxygen isotope fractionation can be removed from δ 18 O s values so that the stable oxygen isotope signature of the water in which the bivalves lived can be computed. This is possible by solving the paleothermometry equation of Grossman & Ku (1986) for δ 18 O wr * (Eq. 2) and compute the oxygen isotope values of the water from those of the shells and river water temperature. Stream water temperature during shell growth can be reconstructed from shell growth rate data (Eq. 3; Schöne et al., 2004a+b, 2005 or 360 instrumental air temperature (Eq. 1; Morrill et al., 2005;Chen & Fang, 2015). Similar studies in which the oxygen isotope composition of microfossils or authigenic carbonate obtained from lake sediments were used to infer the oxygen isotope value of the water, however, merely relied on estimates of the temperature variability during the formation of the diatoms, ostracods, abiogenic carbonates etc., and how these temperature changes affected reconstructions of δ 18 O w values (e.g., Rosqvist et al., 2013). In such studies, it was impossible to reconstruct the 365 actual water temperatures from other proxy archives. Moreover, at least in some of these archives, such as diatoms, the effect of temperature on the fractionation of oxygen isotopes between the skeleton and the ambient water is still debated (Leng, 2006).
(2) M. margaritifera precipitates its shell near oxygen isotopic equilibrium with the ambient water, and shell δ 18 O 370 values reflect stream water δ 18 O data. This may not be the case in all archives, which have previously been used.
(3) The shells can provide seasonally to inter-annually resolved data. In the present study, each sample typically represented as little as one week up to one full growing season (= one 'year'; mid-May to mid-October; Dunca 375 et al., 2005). In very slow growing shell portions of ontogenetically old specimens, individual samples occasionally covered two or, in exceptional cases, three years of growth which resulted in a reduction of variance.
If required, a refined sampling strategy and computer-controlled micromilling could ensure that time-averaging consistently remains below one year. Such high-resolution isotope data can be used for a more detailed analysis of changes in the precipitation-runoff transformation across different seasons. Furthermore, the specific sampling 380 method based on micromilling produced uninterrupted isotope chronologies, i.e., no shell portion of the outer shell layer remained un-sampled. Due to the high temporal resolution, bivalve shell-based isotope chronologies can provide insights into inter-annual and decadal-scale paleoclimatic variability. With the new, precisely calendar-aligned data, it becomes possible to test hypotheses brought forward in previous studies according to which δ 18 O signatures of meteoric water are controlled by the winter and/or summer NAO (e.g., Rosqvist et al., 385 2007Rosqvist et al., 385 , 2013. (4) Each sample taken from the shells can be placed in a precise temporal context. The very season and exact calendar year during which the respective shell portion formed can be determined in shells of specimens with known dates of death based on the seasonal growth curve and annual increment counts. Existing studies suffer from the 390 disadvantage that time cannot be precisely constrained, neither at seasonal nor annual time-scales (unless varved sediments are available). However, isotope results can be biased toward a particular season of the year or a specific years within a decade. Such biases can be avoided with subannual data provided by bivalve shells.
In summary, bivalve shells can provide uninterrupted, seasonally to annually resolved, precisely temporally constrained 395 records of past stream water isotope data that enable a direct comparison to climate indices and instrumental environmental data. In contrast to bivalve shells, sedimentary archives come with a much coarser temporal resolution. Each sample taken from sediments typically represents the average of several years, and the specific season and calendar year during which the ostracods, diatoms, authigenic carbonates etc. grew remains unknown. On the other hand, the time intervals covered by sedimentary archives are much larger and can reveal century and millennial-scale variations with much less effort than 400 sclerochronology-based records. As such, the two types of archives could complement each other perfectly and increase the understanding of past climatic variability. For example, once the low-frequency variations have been reconstructed from sedimentary archives, a more detailed insight into seasonal to inter-annual climate variability can be obtained from bivalve shells. As long as the date of death of the bivalves is known, such records can be placed in absolute temporal context (calendar year). Although the same is currently impossible with fossil shells, each absolutely dated (radiocarbon, amino acid 405 razemization dating) shell of a long-lived bivalve species can open a seasonally to annually resolved window into the climatic and hydrological past of a region of interest.

M. margaritifera shell δ 18 O values reflect stream water δ 18 O values
Unfortunately, complete, high-resolution and long-term records of δ 18 O w values of the studied rivers were not available. Such data is required for a direct comparison with those reconstructed from shells (δ 18 O wr * or δ 18 O wr(SGI) * values) and to determine 410 if the bivalves precipitated their shells near oxygen isotopic equilibrium with the ambient water. However, one of the study sites (NJB) is located close to the river Skellefteälven, where δ 18 O w values were irregularly analyzed between 1973 and 1980 ( Fig. 6a) by the Water Resources Program (GNIR dataset). It should be noted that the δ 18 O w data of GNIR merely reflect temporal snapshots, but not actual monthly averages. In fact, the isotope signature of meteoric water can vary significantly on short time-scales (e.g., Darling, 2004;Leng and Marshall 2004;Rodgers et al., 2005). In addition, for some months, no GNIR 415 data were available. In contrast, shell isotope data represent changes in the isotope composition of the water over coherent time intervals ranging from one week to one year (and in few cases two or three years). Due to the micromilling sampling technique, uninterrupted δ 18 O s time-series were available. It is thus compelling how well the ranges of intra-annual δ 18 O wr * data compared to instrumental oxygen isotope data of the river Skellefteälven (Fig. 6a), and that summer averages as well as growing season averages of shells and GNIR were nearly identical (Table 3). Furthermore, in each studied river, individual 420 δ 18 O wr * series agreed strongly with each other (Fig. 5). All these aspects strongly suggest that shell formation occurred near equilibrium with the oxygen isotope composition of the ambient water, and M. margaritifera recorded changes of stream water δ 18 O values. Our conclusions are in agreement with previously published results from various different freshwater mussels (e.g., Dettman et al., 1999;Kaandorp et al. 2003;Versteegh et al., 2009) and numerous marine bivalves (e.g., Epstein et al., 1953;Mook and Vogel, 1968;Killingley and Berger, 1979). 425

Site-specific and synoptic information recorded in shell oxygen isotopes
Although individual chronologies from a given river compared well to each other with respect to absolute values, the three studied sites differed by almost 2.00 ‰ (average δ 18 O wr * values between 1947 and 1992: NJB = -12.51 ‰; GTB = -12.21 ‰; GJ = -14.16 ‰; Fig. 5+7). If our interpretation is correct and δ 18 O s values of the studied margaritiferids reflect the oxygen isotope signature of the water in which they lived, then these numbers reflect hydrological differences in the upstream 430 catchment controlled by a complex set of physiographic characteristics: catchment size and elevation, transit times, upstream lake size and depth controlling the potential for evaporative depletion in 16 O, river flux rates, river width and depth, humidity, wind speed, groundwater influx, differences in meltwater influx etc. (Peralta-Tapia et al., 2014;Geris et al., 2017;Pfister et al., 2017). However, detailed monitoring would be required to identify and quantify the actual reason(s) for the observed hydrological differences. Thus, we refrain from speculations. 435 Despite site-specific differences described above, the δ 18 O wr * chronologies of the three rivers were significantly positively correlated to each other, suggesting that common environmental forcings controlled isotopic changes in the entire study region.
Previous studies suggest that these environmental forcings may include changes in the isotope composition of precipitation, specifically the amount, origin and air-mass trajectory of winter snow and summer rain, timing of snowmelt as well as 440 condensation temperature (Rosqvist et al., 2013). The latter is probably the most difficult to assess, because no records are available documenting the temperature, height and latitude at which the respective clouds formed. Likewise, we cannot confidently assess the link between the isotope signature of precipitation and stream water because only limited and incoherent datasets are available from the study region. Besides that, data on precipitation amounts were taken from another locality and another time interval. Yet, it is well known that precipitation in northern Scandinavia, particularly during winter, originates 445 from two different sources, the Atlantic and arctic/polar regions (Rosqvist et al., 2013), and the moisture of these air masses is isotopically distinct (Araguás-Araguás et al., 2000;Bowen and Wilkinson, 2002). During NAO + years, the sea level pressure difference between the Azores High and the Iceland Low is particularly large resulting in mild, wet winters in central and northern Europe, with strong westerlies carrying heat and moisture across the Atlantic Ocean toward higher latitudes (Hurrell et al., 2003). During NAOyears, however, westerlies are weaker and the Polar Front is shifted southward, allowing Artic air 450 masses to reach northern Scandinavia. Precipitation originating from the North Atlantic is isotopically heavier (δ 18 O p = -5.00 to -10.00 ‰) than precipitation from subarctic and polar regions (δ 18 O p = -10.00 to -15.00 ‰). Furthermore, changes in air mass properties over northern Europe are controlled by atmospheric pressure patterns in the North Atlantic, particularly the NAO during winter (Hurrell, 1995;Hurrell et al., 2003). The positive correlation between δ 18 O wr * chronologies of the three studied rivers and the wNAO index (Table 4; Fig. 7a-c, 8a) suggests that the shell isotopes recorded a winter precipitation 455 signal, and this can be explained as follows. A larger proportion of arctic air masses carried to northern Scandinavia during winter resulted in lower δ 18 O p values, whereas the predominance of North Atlantic air masses caused the opposite. In NAO + years, strong westerlies carried North Atlantic air masses far northward so that winter precipitation in northern Sweden had significantly higher δ 18 O p values than during NAOyears. When the NAO was in its negative state, precipitation originated predominantly from moisture of polar regions which are depleted in 18 O, and hence have lower δ 18 O p values. The specific 460 isotope signatures in the rivers were controlled by the snowmelt in spring. Essentially, the bivalves recorded the (damped) isotope signal of the last winter precipitation -occasionally mixed to spring and summer precipitation -in their shells. This hypothesis is supported by the correlation of the few available GNIP and GNIR data with the wNAO index (Fig. 8d+e). Rosqvist et al. (2007) hypothesized that the summer NAO strongly influences δ 18 O p values and thus, the δ 18 O w signature of the open, through-flow lakes in northern Scandinavia. However, our data did not support a profound influence of the summer 465 NAO index on δ 18 O wr * values (Fig. 7d-f). This conclusion is consistent with other studies suggesting that the summer NAO has a much weaker influence on European climate than the NAO during winter (e.g., Hurrell, 1995).
Following Baldini et al. (2008) and Comas-Bru et al. (2016), northern Sweden is not the ideal place to conduct oxygen isotopebased wNAO reconstructions. Their models predicted only a weak negative or no correlation between δ 18 O p values and the 470 wNAO index in our study region (Baldini et al., 2008: Fig. 1;Comas-Bru et al. 2016: Fig. 3a). One possible explanation for this weak correlation is the limited and temporally incoherent GNIP dataset in northern Sweden, from which these authors extracted the δ 18 O p data that were used to construct the numerical models. In contrast, δ 18 O data of diatoms from open lakes in northern Sweden revealed a strong link to the amount of precipitation and δ 18 O p values, which reportedly are both controlled by the predominant state of the NAO (Hammarlund et al., 2002;Andersson et al., 2010;Rosqvist et al., 2004Rosqvist et al., , 2007Rosqvist et al., , 2013. 475 Findings of the present study substantiated these proxy-based interpretations. Furthermore, we presented, for the first time, oxygen isotope time-series with sufficient temporal resolution (annual) and precise temporal control (calendar years) required for a year-to-year comparison with the NAO index time-series.
As Comas-Bru et al. (2016) further suggested, the relationship between δ 18 O p values and the wNAO index is subject to spatial 480 non-stationarities, because the southern pole of the NAO migrates along a NE-SW axis in response to the state of another major atmospheric circulation mode in the North Atlantic realm, known as the East Atlantic Oscillation or East Atlantic Pattern (EA) (Moore and Renfrew, 2012;Moore et al., 2013;Comas-Bru and McDermott, 2014). Like the NAO, the EA is most distinct during winter and describes atmospheric pressure anomalies between the North Atlantic west of Ireland (low) and the subtropical North Atlantic (high). Through interaction of these circulation patterns, the correlation between the wNAO and 485 δ 18 O p values can weaken at times in certain regions. For example, when both indices are in their positive state, the jet stream shifts poleward (Woolings and Blackburn, 2012) and the storm trajectories that enter Europe in winter take a more northerly route (Comas-Bru et al., 2016). The δ 18 O p values will then be lower than during NAO + /EAyears. To identify whether this applies to the study region, we followed Comas-Bru et al. (2016) and tested if the relationship between the wNAO and reconstructed river water oxygen isotope data remain significant during years when the signs of both indices are the same (EQ) 490 and during years when they are opposite (OP) (note, the EA index is only available from 1950 onward). As demonstrated in

Damped river water oxygen isotope signals
Compared to the large isotopic difference between winter precipitation sourced from SW or N air masses, the huge seasonal  (Fig. 6) and inter-annual time-scales (Fig. 5a-c). This figure agrees well with seasonal amplitudes determined in other rivers at higher latitudes of the Northern Hemisphere (Halder et al. 2015) and can broadly be explained by catchment damping effects due to water collection, mixing, storage and release processes in upstream lakes and groundwater from which these rivers were fed. Catchment mean transit time (MTT), determined via a simple precipitation vs. stream flow 505 isotope signal amplitude damping approach (as per de Walle et al., 1997), is approximately 6 months -corroborating the hypothesis of a mixed snowmelt and precipitation contribution to the stream water δ 18 O signal during the growing season.
The attenuated variance on inter-annual time-scales can possibly be explained -amongst others -by inter-annual changes in the amount of winter precipitation and the timing of snowmelt. Colder spring temperatures typically resulted in a delayed 510 snowmelt so that lower oxygen isotope signatures still prevailed in the river water when the main growing season of the bivalves started. However, winter precipitation amounts remained below average in NAOyears, so that the net effect on δ 18 O w values in spring was less severe than the isotopic shift in δ 18 O p values. In contrast, the amount of snow precipitated during NAO + years was larger, but milder spring temperatures resulted in an earlier and faster snowmelt, so that the effect on the isotope signature of river water at the beginning of the growing season of the mussels likely remained moderate. 515

Sub-annual dating precision and relative changes of seasonal shell growth rate
The precision by which the time can be determined that is represented by individual isotope samples depends on the validity of the seasonal growth model. We assumed that the timing of seasonal shell growth was similar to published data of M. margaritifera and remained the same in each year and each specimen. This may not be entirely correct, because the timing and rate of seasonal shell growth can potentially vary between localities, among years and among individuals, yet in M. 520 margaritifera the seasonal timing of shell growth is remarkably invariant across large distances (Dunca et al., 2005). A major dating error exceeding four weeks, however, seems unlikely because the oxygen isotope series of individual specimens at each site were in good agreement. Presumably, the timing of seasonal shell growth is controlled by genetically determined biological clocks, which serve to maintain a consistent duration of the growing season (Schöne, 2008). Although shells grew faster in some years and slower in others, the relative seasonal changes in shell growth rates likely remained similar and consisted of a 525 gradual increase as the water warmed and more food became available in spring and summer followed by a gradual decline as temperatures dropped in fall. It was further assumed that the timing of shell growth has not significantly changed through the lifetime of the studied specimens. In fact, if ontogenetic changes of seasonal growth traits had occurred, it would be impossible to crossdate growth curves from young and old individuals and construct master chronologies (Schöne et al., 2004a(Schöne et al., , b, 2005Helama et al., 2006;Black et al., 2010). Based on these arguments, seasonal dating errors were likely minor. 530

Shell stable carbon isotopes
Our results are consistent with previous studies using long-lived bivalves (Beirne et al., 2012;Schöne et al., 2005cSchöne et al., , 2011 where δ 13 C s chronologies of M. margaritifera did not show consistent ontogenetic trends, but rather oscillated around an average value (ca. -12.00 to -13.00 ‰). The time-series of NJB were too short to reject the hypothesis of directed trends through lifetime, but we propose here that the δ 13 C s values of shells from that river would also average out at ca. -12.50 ‰ like 535 at the other two studied sites, if longer chronologies were available. If a contribution of metabolic CO 2 to the shell carbonate exists in this species (which we cannot preclude because no δ 13 C values of the dissolved inorganic carbon (DIC) data are available for the studied rivers), it likely remains nearly constant through lifetime as it does in other long-lived bivalve mollusks (Schöne et al., 2005c(Schöne et al., , 2011Butler et al., 2011;Reynolds et al., 2017). Observed stable carbon isotope signatures in the mussel shells are in the range of those expected and observed in stream waters of northern Europe (-10.00 to -15.00 ‰; Leng and 540 Marshall, 2004) Seasonal and inter-annual changes of δ 13 C s values could be indicative of changes in primary production, food composition, respiration and influx of terrestrial detritus. However, in the absence of information on how the environment of the studied rivers changed through time, we can only speculate about possible causes of temporal δ 13 C DIC variations. For example, 545 increased primary production in the water would not only have propelled shell growth rate, but also resulted in a depletion of 12 C in the DIC pool and thus, higher δ 13 C DIC and δ 13 C s values. However, just the opposite was observed on seasonal and interannual time-scales. Highest δ 13 C s values often occurred near the annual growth lines, i.e., during times of slow growth, and, although not statistically significant, annual δ 13 C s(d) * values at NJB and GTB were inversely related to shell growth rate ( Fig.   7g+h; Table 4). Accordingly, δ 13 C s(d) * values do not seem to reflect phytoplankton dynamics. Another possibility is that a 550 change in the composition of mussel food occurred which changed the shell stable carbon isotope values without a statistically significant effect on shell growth rate. Since the isotope signatures of potential food sources differ from each other (e.g., Gladyshev, 2009), a change in relative proportions of phytoplankton, decomposing plant litter from the surrounding catchment vegetation, bacteria, POM deriving from higher organisms etc. could have left a footprint in the δ 13 C s(d) * values. Furthermore, seasonal and inter-annual changes in respiration or influx of terrestrial detritus may have changed the isotope signature of the 555 DIC pool and thus the shells. Support for the latter comes from the weak negative correlation between δ 13 C s(d) * values and the wNAO (Tab. 4; without Bonferroni correction, p-values remained below 0.05). After wet (snow-rich) winters (NAO + years), stronger terrestrial runoff may have flushed increased amounts of light carbon into the rivers which lowered δ 13 C DIC values.
To test these hypotheses, data on the stable carbon isotope signature of digested food and DIC would be required, a task for subsequent studies. 560

Error analysis and sensitivity tests
To test the robustness of the findings presented in Tables 3 and 4 as well as their interpretation, we have propagated all uncertainties associated with measurements and modeled data and randomly generated δ 18 O wr *, δ 18 O wr(SGI) *, δ 18 O wr(Norrland) * and δ 13 C s(d) * chronologies (via Monte Carlo simulation). A brief overview of the errors and simulation procedures are provided here below. 565 Water temperature estimates (Eq. 1) were associated with an error (one standard deviation) of ± 2.07 °C. Amongst others, this large uncertainty results from the combination of temperature data of four different rivers, which all varied with respect to average temperature and year-to-year variability. The error exceeds the inter-annual variance (one standard deviation of ± 0.90 °C) of the instrumental water temperature average (8.64 °C) more than two times. Instead of reconstructing T w from T a with 570 an uncertainty of ± 2.07 °C, we could have assumed a constant water temperature value of 8.64 °C with an uncertainty of only ± 0.90 °C. However, our goal was to improve the δ 18 O wr * reconstructions by taking into account the actual year-to-year temperature variability. To simulate the effect of different temperature uncertainties, we randomly generated 1,000 T w * chronologies with an error of ± 0.90 °C as well as 1,000 chronologies with ± 2.07 °C. Both sets of simulated T w * time-series were used in subsequent calculations. Errors involved with shell growth patterns include the measurement error (± 1 µm 575 equivalent to an SGI error of ± 0.06 units) and the variance of crossdated SGI data. In different calendar years, the standard error of the mean of the 25 SGI chronologies ranged between ± 0.03 and ± 0.66 SGI units. The measurement and crossdating uncertainties were propagated, 1,000 new SGI chronologies randomly generated and regressed against simulated T w * chronologies. The uncertainty of the new SGI vs. T w * model (standard error of ± 1.35 °C) was propagated in subsequent calculations of δ 18 O wr(SGI) * values using Eq. 2. A third set of uncertainties was associated with isotope measurements (analytical 580 precision error, one standard deviation = ± 0.06 ‰), the calculation of site-specific annual averages from contemporaneous specimens (for δ 18 O, on average: ± 0.11 ‰ to ± 0.15 ‰; for δ 13 C, on average: ± 0.37 ‰ to ± 0.42 ‰), and the calculation of the Norrland average. All errors were propagated and new δ 18 O wr *, δ 18 O wr(SGI) *, δ 18 O wr(Norrland) * and δ 13 C s(d) * chronologies simulated (1,000 representations each). The simulated chronologies were regressed against NAO and SGI chronologies (results of sensitivity tests for regressions of δ 18 O wr(SGI) * and δ 18 O wr(Norrland) * values vs. wNAO indices are given in Table 5). 585 According to the complex simulation experiments, the observed links between reconstructed stream water oxygen isotope values and the wNAO largely remained statistically robust disregarding of which T w * error was assumed (Table 5). This outcome is not particularly surprising given that even the annual δ 18 O s chronologies of the studies specimens were strongly coherent and values fluctuated at time-scales similar to that of the wNAO (Fig. 4). Apparently, decadal-scale atmospheric 590 circulation patterns exert indeed a strong control over the isotope signature of river water in the study area. However, none of the correlations between shell isotope data and sNAO attained significance level. The importance of summer rainfall seems much less important for the isotope value of river water than winter snow. As before, the relationship between stable carbon isotope data of the shells and climate indices as well as shell growth rate remained weak and were statistically not significant.
Inevitably, the propagated errors, specifically the uncertainty associated with the reconstruction of the stream water 595 temperature from air temperature resulted in a notable drop of explained variability and statistical probability (Tab. 5). The use of instrumental water temperatures could greatly improve the reconstruction of δ 18 O wr * values, because the measurement error would be on the order of 0.1 °C instead of 2.07 °C or 0.90 °C. Future calibration studies should therefore be conducted in monitored rivers.

Summary and conclusions 600
Stable oxygen isotope values in shells of freshwater pearl mussels, M. margaritifera, from rivers in northern Sweden mirror stream water stable oxygen isotope values. Despite a well-known damping of the precipitation signal in stream water isotope records, these mollusks archive local precipitation and synoptic atmospheric circulation signals, specifically the NAO during winter. Stable carbon isotope data of the shells are more challenging to interpret, but they seem to record local environmental conditions such as changes in DIC and/or food composition. Future studies should be conducted in rivers in which temperature, 605 DIC and food levels are closely monitored to further improve the reconstruction of stream water δ 18 O values from δ 18 O s data and better understand the meaning of δ 13 C s fluctuations.
The bivalve shell oxygen isotope record presented here extends back to AD 1819, but there is potential for developing longer isotope chronologies through the use of fossil shells of M. margaritifera collected in the field or taken from museum 610 collections. With suitable material and by applying the crossdating technique, the existing chronology could probably be extended by several centuries back in time. Stream water isotope records may shed new light on pressing questions related to climate change impacts on river systems, mechanistic understanding of water flow and quality controlling processes, calibration and validation of flow and transport models, climate and earth system modeling, time variant catchment travel time modeling, etc. Longer and coherent chronologies are essential to reliably identify multidecadal and century-scale climate 615 dynamics. Even individual radiocarbon dated fossil shells that do not overlap with the existing master chronology can provide valuable paleoclimate information, because each specimen of M. margaritifera can open a seasonally to annually resolved, multi-year window into the history of rivers.
Acknowledgements. We thank Denis Scholz and Erika Pietroniro for constructive discussions. We are grateful for comments 620 and suggestions provided by two anonymous reviewers that greatly improved the quality of this article.
Data and code availability. All data and code used in this study are available from the authors on request. Additional supplementary files are available at https://www.paleontology.uni-mainz.de/datasets/HESS_2019_337_supplements.zip.

625
Sample availability. Bivalve shell samples are archived and stored in the paleontological collection of the University of Mainz.
Author contributions. BRS designed the study, performed the analyses and wrote the manuscript; AEM and SMB conducted the field work and collected samples; SMB sampled the shells and temporally aligned the isotope data; JF isotopically analyzed 630 the shell powder; LP conducted MTT calculations. All authors jointly contributed to the discussion and interpretation of the data.
Competing interests. The authors declare that they have no conflict of interest.

635
Financial support. This study has been made possible through a research grant by the Deutsche Forschungsgemeinschaft DFG to BRS (SCHO793/1). T a -instrumental air temperature T w -river water temperature reconstructed from T a T w * -weighted (considering variations in seasonal shell growth rate) river water temperature reconstructed from SGI and T w ; annual T w * = growing season mean value 680 Table 1: Shells of Margaritifera margaritifera from three rivers in northern Sweden used in present study for isotope and growth 955 pattern analysis. Last set of digits denotes whether bivalves were collected alive (a) or dead (d), specimen number and which valve was used (R = right, L = left). § = minimum estimate of lifespan. i = last sampled year incomplete. A = Add ten years to these values to obtain approximate ontogenetic years.   Table 3. Relationship between stable oxygen isotope values in precipitation (amount-corrected δ 18 Op), river water and shells of Margaritifera margaritifera from Nuortejaurbäcken during different portions of the year (during four years for which data from shells, water and precipitation were available : 1975, 1976, 1978, 1979; hence N = 4). The arithmetic mean δ 18 O values for each portion of the year are also given. The rationale behind the comparison of δ 18 O values of winter precipitation and spring (May-June) river water or shell carbonate is that the isotope signature of meltwater may have left a signal in the water.  Table 4. Site-specific annual isotope chronologies of Margaritifera margaritifera shells linearly regressed against winter and summer NAO (wNAO, sNAO) as well as detrended and standardized shell growth rate (SGI). δ 18 Owr* data were computed from shell oxygen isotope data and temperature data computed from instrumental air temperatures, whereas in the case of δ 18 Owr(SGI)* data, temperatures were estimated from a growth-temperature model.    Tables 3 and 4, uncertainties (one of them the error associated with the reconstruction of stream water temperatures, Tw, from air temperatures, 990 Ta) were propagated and used to randomly generate δ 18 Owr(SGI)* chronologies, which were subsequently regressed against winter North Atlantic Oscillation (wNAO) indices. Simulations were computed with propagated Tw* of 2.07 °C and 0.90 °C. See text for details. Statistically significant values (Bonferroni-adjusted p < 0.05) are marked in bold.   Nuortejaurbäcken and Grundträsktjärnbäcken that were sampled with very high spatial resolution and from which the majority of isotope data were obtained (Table 1)