Total alkalinity (TA) was determined on the same samples as dissolved nutrients (Mori et al., 2021). The respective samples were taken daily after the hourly mixing via ultrapure water (purified on the Milli-Q system Millipore, 18.2 MΩ x cm) and sample rinsed syringes attached to low density polyethylene (LDPE) tubing. Samples were filtered directly through 0.45 µm pore size Surfactant-Free Cellulose Acetate (SFCA) membrane filters, collected into pre-cleaned 5mL Eppendorf tubes and toxified with HgCl₂ to a final concentration of 0.01 vol%. Until further analysis, samples were stored in the cool and dark. For TA analysis we used a multiscan GO microplate spectrophotometer (Thermo Scientific) and followed a modified method by Sarazin et al. (1999).

Samples for dREE analyses were taken in three replicates (P2 to P4) at intervals of 1-5 days after the hourly mixing from a depth of about 30-40 cm below surface. Sampling was performed using LDPE tubes and syringes (25 mL), pre-rinsed with ultrapure water and sample. Samples were filtered directly through 0.45 µm pore size SFCA membrane filters, collected into acid pre-cleaned LDPE bottles and acidified with ultra-clean concentrated HNO₃ to 2 vol. % and a final pH of ≤2.

For the analysis of dREEs, we diluted 3 mL of each sample with 8 mL ultrapure water, which was acidified to pH 2 by adding ultraclean HNO_3. The diluted sample was then further adjusted to pH 2, if necessary, using ultraclean HNO₃ and/or NH₄OH The sample handling was carried out in a clean lab. Nitric acid was purified by single sub-boiling distillation and NH₄OH was purified by evaporation and absorption into ultrapure water under ultraclean conditions. Preconcentration, analysis and quantification of dREEs followed the method described in Behrens et al. (2016). In brief, all samples were spiked with a multi-element REE isotope spike (DKM, prepared by Gilbert N. Handson of SUNY Stony Brook). After an equilibration time of at least 48 hrs, the samples were purified and preconcentrated using the seaFast-pico™.system (Elemental Scientific) in offline mode. Elemental analyses were carried out via inductively coupled plasma mass spectrometry (ICP-MS, Thermo Finnigan Element II). We introduced the samples into the ICP-MS using an autosampler (CETAC ASX-100) and a desolvation introduction system (CETAC Aridus II) to enhance sensitivity and lower oxide formation rates (<0.05% for Ce). We therefore did not apply any correction for oxide formation to the measured values. Concentrations of dREEs were quantified via isotope dilution and additional measurements of external standards (Coral Sea Standard in concentrations of 40, 55, 65, 80 and 100 ppt). Seawater standard SAFe 3,000 m (North Pacific) was purified and analysed regularly (about 2-3 times) in each seaFAST and ICP-MS session. Procedural blanks were prepared and analysed in the same manner as samples using 11 mL of ultrapure water at pH 2.

External reproducibility, that we determined by repeated preconcentration and analysis of seawater samples from North Pacific station SAFe 3,000 m, was <3% for all elements except Gd (3.5%) and Ce (8.1%). Apart from La, all REEs were within the confidence interval of the long-term laboratory mean of SAFe 3,000 m measurements published by Behrens et al. (2016), even though our values for Ce were much lower than previously reported ( Table 1 ). Procedural blanks accounted for <0.5% of sample concentrations ( Table 1 ), sample concentrations were therefore not blank-corrected.

The group of REEs was subdivided into three individual subgroups: LREEs (La, Ce, Pr, Nd), MREEs (Eu, Gd, Tb, Dy, Ho) and HREEs (Er, Tm, Yb, Lu). This is possible without losing valuable information, since the behaviour of individual REEs within the subgroups is quite similar, which was confirmed within this study. The element Sm was excluded from all calculations and discussions since its concentration over time was significantly different compared to the other REEs (see also results below). In this study, we want to emphasise the common observations for REEs and focus on more general processes between REEs and OM.

In order to illustrate changes throughout the incubation period, we normalized the dREE concentrations of each time point to its initial concentration (T0) in the respective replicate. Fractionation was calculated via the ratio of HREE/LREE = (Lu+Yb+Tm)/(La+Pr+Nd) using values normalized to PAAS (Rudnick and Gao, 2003). The uncertainties for HREE/LREE, as calculated from the SAFe 3,000m replicates, was on average 0.06

Losses and gains of dREEs were calculated by dividing the minimal (maximal) concentration by the maximum (minimum) of each phase. A phase is defined as a period within the total incubation time, where the dREE concentration trend is consistent, i.e., no shift from decreasing (increasing) concentrations to the opposite occurred (see section 4.2).

The intercorrelations of the various biogeochemical parameters presented in the context of this study were determined by a three-step approach combining a Spearman correlation analysis with a cluster and a non-metric multidimensional scaling (NMDS). All statistical analyses were carried out using R version 4.1.0 (R Core Team, 2016) and performed on a combined dataset for all three biological replicates (P2-P4).

In a first step, we performed a correlation analysis, where we correlated each biogeochemical parameter (n) with each other, creating a n x n similarity matrix, based on Spearman’s rho (rho), using the rcorr function of the Hmisc Rpackage (Harrel and Dupont, 2019). The corresponding correlation matrix was visualized using the corrplot function by corrplot R-Package (Wei and Simko, 2007). The parameter rho indicates correlation based on ranks, where +1 indicates a strong positive and -1 indicates a strong negative correlation (Spearman, 1904). As significant, we considered correlations with a p value of <0.1.

In a second step we combined a hierarchical clustering analysis with a NMDS using the complete linkage algorithm. Clusters are classically shown as dendrograms a NMDS, however, is able to visualize intercorrelations of parameters in a two-dimensional space. We used a 1-|rho| distance matrix, so that parameters with strong positive and negative correlations plot close to each other. To determine the optimal number of clusters (k_optimal), we used the maximal average silhouette value (Kaufman and Rousseeuw, 2005; Borcard et al., 2018). Setting a maximum of k_max=9 clusters, the optimal number of clusters for the present dataset was k_optimal=4 ( Supplementary Figure 2 and 3). We visualized the resulting dendrogram using the dendextend package (Galili, 2015). The NDMS was calculated using the metaMDS function of the vegan package (Oksanen et al., 2019).

The interplay between dREEs and dissolved ligands, with special emphasis on organic complexation, was determined based on a model approach using the aqueous geochemical modelling program PHREEQC version 3 (Parkhurst, 1995). This program can simulate aqueous species of a specified solution composition by using equilibrium chemistry. The needed thermodynamic data is defined in an associated database. Since no predefined database contained all the stability constants needed for the modelling of chemical speciation of dREEs under differing elemental concentrations, we constructed a new database following Schijf et al. (2015), including desferrioxamine B (DFOB) as the representative of a strong organic ligand. A detailed description of the database as well as the applied model script is outlined in the supplementary material B. Many studies have modelled the organic complexation of dREEs with terrestrial humic and fulvic acids (HFAs) (Tang and Johannesson, 2003; Sonke and Salters, 2006; Pourret et al., 2007). Catrouillet et al. (2020) showed that microbiologically produced HFAs, in contrast to those of terrestrial origin, preferentially complex HREEs and thus show different complexation patterns with dREEs. Since the existence of terrestrial HFAs can be excluded in our mesocosm approach and the DOC pool in general is probably less complex than in natural coastal waters due to the absence of large algae and a refractory DOC pool, we decided to use only one well described microbially produced organic ligand, DFOB, in our model. This one molecule is not able to reproduce the complexity of DOC, but Schijf et al. (2015) mention that its REE complexation pattern is likely to be typical for various organic ligands. Therefore, the model is expected to be able to show general trends in the complexation of REEs with microbially produced DOC.

The initial modelling step is to characterize the artificial seawater, based on concentrations of the main seawater ions, as well as TA and pH. Each dREE data point is considered in the model as a separate reaction step and the measured dREE concentrations with it entered. The parameters TA and pH are re-defined for each individual reaction step due to their high importance in chemical REE speciation (Byrne, 2002; Liu and Millero, 2002; Luo and Byrne, 2004). In order to focus on the complexation of dREEs, precipitation of particulate phases as well as gas exchange and interfering elements (i.e., Fe and Mn) were not considered in this approach. With entry of the measured dREE concentrations via the individual reaction steps, inorganic REE-complexes with the seawater components and free dREE ions were calculated. To account for and calculate additional REE-DOC complexes, we entered DOC-concentrations to the model by equating parts of the DOC pool with the strong organic ligand DFOB. The final output of a model simulation run were concentrations of inorganic and organic complexes, and free ions, as well as the total concentration for each respective dREE for each reaction step and thereby each data point.

Since the specific fraction of the complex DOC pool, which contains strong organic ligands for REEs and therefore equals DFOB, is not known, we used two different assumptions following Schijf et al. (2015). Thereafter, the definition of the ‘High-DOC’ and ‘Low-DOC’ approach follows maximal simulated REE-DOC proportions to the dREE pool. The ‘High-DOC’ approach results in a maximum of 40% Lu-DOC of total dLu by assuming that strong organic ligands occur with a concentration of c = 10^-6 * c_m (c_m = measured DOC concentration at specific sampling point, that is reaction step). The ‘Low-DOC’ approach results in a maximum of 10% Lu-DOC by assuming concentrations of c = 2*10^-7 * c_m for strong organic ligands. These assumptions and associated calculation of concentrations were not changed throughout the incubation period, assuming that the parts of the DOC pool acting as strong REE ligands were stable over the entire incubation period.

Values for TA fluctuated between 2.6 and 2.9 mmol L^-1 with a peak around day 8 to 15 in all biological replicates ( Supplementary Figure 1B ).

In comparison to spring conditions in the coastal North Sea (Paffrath et al., 2020), the mesocosm matrix was up to 4-times higher in dLREE concentrations, while dHREEs were half the natural concentrations. Overall, however, the dREE concentrations in the mesocosms were in line with natural relative abundances with enrichment of PAAS-normalized dHREEs relative to dLREEs ( Supplementary Figure 4 ).

In general, the concentrations of dLREEs, dMREEs and dHREEs showed a similar time course, with the largest changes over time with respect to T0 for dLREEs and smallest for dHREEs. The temporal concentration course of dSm differed distinctly from the other dREEs, with concentrations increasing constantly from around day 20 onwards ( Supplementary Figure 5 ). In previous studies dissolved Sm anomalies have been attributed to contamination either during analysis (Möller et al., 2003; Zhu et al., 2006; Lawrence and Kamber, 2007; Takata et al., 2009) or to anthropogenic sources (Kulaksız and Bau, 2013). In the present study, external anthropogenic contamination is, based on the closed incubation set-up, unlikely. Contamination from the set-up itself and the analytical procedures can be excluded since neither the control experiments nor the REE standards show positive Sm anomalies. Thus based on the available dataset, any statement regarding the reason for high or rather increasing Sm concentrations in the biological replicates are highly speculative and are thus not further discussed within the framework of this study. The focus of the present study, is solely on the common REE trends and the drivers responsible for those. Thus, as mentioned above, Sm was not included in the LREE subgroup and will not be discussed further in the context of this study. The general pattern for dREEs can be described as decreasing concentrations during the diatom bloom followed by increasing concentrations during the subsequent bacterioplankton maximum. Concentrations decreased again during the Phaeocystis sp. bloom. During its stationary phase, dREE concentrations increased again in P2 and P4, while concentrations continued to decrease in P3 ( Figure 3A ).

Within the biota-free controls (C1 and C2), dREE concentrations were slightly higher compared to the biological replicates at T0. At the end of the incubation in the controls, we observed a 40% and 21% loss of dLREEs and dHREEs, respectively ( Supplementary Figure 5 ).

Based on the overall course of dREE concentrations in P2-P4, the incubation time is subdivided into four characteristic phases ( Figure 2A ):

Phase I: dREE concentrations decrease in P2, P3, P4; day 1-11

Phase II: dREE concentrations increase in P2, P3, P4; day 11-23 (P2, P4), day 11-26 (P3)

Phase III: dREE concentrations decrease in P2, P3, P4; day 23-30(P2, P4); day 26-38 (P3)

Phase IV: dREE concentrations increase in P2 and P4; day 30-38

A normalization to T0 displays enrichments (gains) and depletions (losses) over the total incubation time relative to the starting point, with values of under 1 showing a general loss of dREEs during the experiment. Values of above 1 show a gain of dREEs and were observed for dLREEs and dMREEs from phase II onwards in all biological replicates ( Figure 2A ). The losses and gains of dREEs within each phase throughout the incubation time show differences within and between replicates. One notable observation is the variations in dHREE behaviour in P2. During phase III (Phaeocystis sp. bloom), the loss of dHREEs in P2 was distinctly lower (-4%) than during phase I (-17%), while the loss of dLREEs (~ -37%) was similar in both phases. In contrast, in phase IV (stationary Phaeocystis sp. bloom) dHREEs increased strongly (+9%), while dLREEs were roughly stagnant (-2%) ( Figure 2B ).

The HREE/LREE ratios mirrored the overall temporal dREE concentrations, with increasing HREE/LREE ratios for increasing dREE concentrations and vice versa ( Figure 3A ). Comparison of HREE/LREE ratios with dNd concentrations, as the representative for the LREEs, revealed a linear relationship for the data points with DOC concentrations <200 µmol L^-1 (R² = 0.98, p=1.1x10^-6 (P2); R^{2 =} 0.93, p= 4.9x10^-6 (P3); R² = 0.09, p= 0.24 (P4)). At times with DOC concentrations >200 µmol L^-1, reached during the end of phase III and during phase IV in P2 and P4, the HREE/LREE ratios were above the linear regression line ( Figure 3B ).

The correlation analysis based on Spearman’s rho ( Figure 4 ) revealed weak positive correlations (0.1<rho <0.4) between the dREEs (especially dHREEs) and directly spring bloom-associated parameters, like POM, DOC and Chl-a. Additionally, positive correlations (rho > 0.4) of the dREE pools with the dFe pool were identified. The combinatorial statistical approach of cluster and NMDS analyses showed an intercorrelation between the three individual dREE subgroups and the dFe pool, that are grouped in one cluster and plotted close to each other and far from other parameters in the NMDS ( Figure 5 ).

The simulations of the applied PHREEQC model showed variations in the proportions of organic complexes (DOC as ligand, REE-DOC complexes) to the dREE pool, that seem to be driven by DOC concentrations and the atomic mass of REEs ( Figure 6 ). In general, the proportions of REE-DOC complexes increased with the onset of the Phaeocystis sp. bloom, in line with increasing DOC concentrations. Depending on the DOC concentration, the proportions of the REE-DOC complexes differed among the three individual biological replicates and over the temporal course. Highest proportions were simulated for the incubation end in P2, which showed highest DOC concentration. While the proportion of REE-DOC were of minor importance for the dLREE pool, it increased for the dMREEs and dHREEs, with highest proportions simulated for ytterbium, the second heaviest element. The ‘Low-DOC’ ( Figure 6B ) and ‘High-DOC’ approaches ( Figure 6A ) showed the same temporal and element-specific variations in organic complexation, only differing in the overall proportion of REE-DOC complexes to the total dREE pool. Organic complexes displaced carbonate complexes, which were identified to be the dominant inorganic species.

The chemical speciation of trace metals, including REEs, has a strong influence on their behaviour in aquatic systems. The speciation itself is highly dependent on temperature and pH as well as salinity and TA, as they influence the stability constants of solution complexes as well as the availability of ligands in the respective water mass (Henderson, 1984; Byrne, 2002; Liu and Millero, 2002; Luo and Byrne, 2004; Quinn et al., 2006; Quinn et al., 2007). For the speciation of dREEs, TA is particularly important, as it determines the availability of carbonate complexes, the dominating inorganic dREE complexes. In the present study, the initial salinity (at T0) as well as the temporal course of pH and TA showed small temporal changes that were similar in all biological replicates ( Supplementary Figure 1 ). Our simulations of the chemical speciation of dREEs using the PHREEQC model showed, that the small changes of the measured TA of ±0.3 mmol L^-1 did not influence the existence of REE-carbonate complexes. A model approach with the measured TA entered and fixed for each individual sampling point/reaction step simulated the exact same results as one using a constant mean TA throughout the simulations. The pH values increased from <8.5 up to 9.5 throughout the biological experiment. For dREEs, studies show that the sorption ability on iron hydroxides, one of the main inorganic scavengers for REEs in seawater, increases with higher pH (Quinn et al., 2006). Solubility of iron oxides appear to stay constant at the pH range observed in our experiments (Liu and Millero, 2002). Moreover, contrary to an expected rise in dREE with decreasing pH, we observe a concentration increase in dREE with increasing pH (phase III in P2 and P4). These observations therefore argue against a major role of pH in dREE behaviour in this study. ( Supplementary Figure 1 ), even though we do not rule out any influence of pH variability.

The main inorganic scavengers for REEs in seawater are Fe and Mn oxides (Elderfield and Greaves, 1982; Nozaki, 2001) due to their large surface area. They can adsorb REEs, acting as a transfer shuttle between the dissolved and particulate REE pool. To investigate the transfer processes associated exclusively with the “abiotic shuttles” and their effects on the dREE pool, in addition to the biotic incubations, two biota-free controls were studied in parallel. The absence of OM in the mesocosms makes it possible to study the influence of abiotic drivers on dREE behaviour. The controls further served to identify possible contamination processes of the stainless steel tanks regarding trace elements. No contaminating input of dREEs or dFe was observed. Iron as well as the dLREEs decreased over the control experiment by 40%, while dHREEs decreased by 20%. Dissolved Mn, on the other hand, increased by 7%. For Fe it was suggested that the loss is driven by adsorption directly onto the mesocosm walls or by precipitation as oxides in the water column (Mori et al., 2021). Similar processes or co-precipitation might also influence the dREEs. The observed dREE fractionation is in line with the higher particle reactivity of LREEs compared to HREEs (e.g. Elderfield and Greaves, 1982; Nozaki, 2001). The biota-containing replicates, in contrast, did not show overall decreasing dREE concentrations, which would indicate adsorption processes on mesocosm walls or on oxides.

These metal oxides were likely formed prior to T0 (see section 2 above) and partly dissolved during phase II. Release of dFe during the same time frame was also attributed to oxide solution (Mori et al., 2021). The behaviour of dMn and dFe in the biological replicates do not support the production of significant amounts of oxides (Mori et al., 2021, see also section 2 above) and thereby limits the possibility of adsorption on inorganic particles as a main transfer-shuttle of the REEs within the mesocosms.

In summary, we conclude that the observed dREE concentration trends are only marginally influenced by abiotic processes/shuttles and that other processes must be more important for driving the overall dREE changes over the incubation period.

Given the dominance of biological cycling and exclusion of other abiotic processes and external influences in the mesocosms, we have to consider bio- or rather bloom-associated processes/transfer shuttles to explain the behaviour of REEs. To determine the nature and effects of these bio-associated shuttles in the mesocosms, we discuss the biological processes which have a possible influence on REEs in the following. Since individual bloom-phases coincided with the identified phases of REE behaviour (phase I-IV), each phase is discussed separately.

The organic complexation of dREEs with parts of the complex DOC pool and the potential influence on dREE abundances and fractionation patterns were evaluated using the dREE speciation model in PHREEQC. The model results show that organic REE-DOC complexes can make up a large proportion of the dREE pool. The simulated proportion of REE-DOC complexes significantly increased over the course of the incubation, in line with increasing DOC concentrations, as well as with atomic mass of the REEs. Modelled organic complexation of LREEs thereby seemed irrelevant for the total dLREE pool, contrasting the high relevance of organic complexation for dHREEs. This is in good accordance with previous model studies (Christenson and Schijf, 2011; Schijf et al., 2015) as well as field studies (Johannesson et al., 2004), in which organic complexation increases in relevance with atomic number of REEs. It confirms that OM-cycling and associated existence of DOC influences the distribution and fractionation of dREEs in marine environments.

Stability constants of organic complexes with DFOB, the organic component that was used as a representative for strong organic REE ligands in the applied model, exceed the ones for inorganic carbonate complexes by a factor of 10⁸ (Christenson and Schijf, 2011). As a result, REE-DOC complexes are highly stable and buffered from adsorption processes. Since the HREEs show highest potential for organic complexation, dHREEs are strongly stabilized while LREEs are not influenced by the presence of organic ligands. This explains the gain of dHREEs in P2 and P4 during phase IV, where high DOC concentrations (>200 µmol L^-1) were reached, and the reduced loss of dHREEs in P2 during phase III with high DOC concentrations (mean of 215 µmol L^-1) compared to phase I (mean of 40 µmol L^-1).

Additionally, stabilization of dHREEs in organic complexes explains the HREE/LREE ratios in P2 and P4 at the end of phase III and in phase IV that are strongly elevated over those in phase II, indicating organic HREE complexation in addition to preferential LREE scavenging. This is also evident in HREE/LREE vs. dNd plots, with Nd representing the LREEs ( Figure 3B ). While in general a linear relationship is observed, indicating that the fractionation is mainly driven by scavenging of LREEs, HREE/LREE ratios for DOC concentrations >200 µmol L^-1 (end of phase III and phase IV) lie above the linear trend line. Based on our results, the concentration of DOC has to exceed a threshold level of 200 µmol L^-1 to cause an alteration in REE patterns. However, we did not observe a fractionation pattern specific to REE-DOC complexation that would allow identification of organic REE complexation in nature.

Many organic ligands may be able to complex other trace metals with higher affinities for the respective binding agent besides the REEs. The organic ligand used in the model approach is DFOB, a siderophore that is bacterially produced to bind iron and specifically mediate its transport inside the cell (Boukhalfa and Crumbliss, 2002). Complexation of REEs by DFOB could occur, if the siderophore is available in excess of free dFe (Christenson and Schijf, 2011; Schijf et al., 2015). In marine environments, this condition seems to be associated with maxima in Chl-a, occurring during phytoplankton blooms (Wagener et al., 2008). The model does not consider the influence of Fe or other trace metals on REE cycling, so this competition effect is not included in the simulations. We assume, that the DOC we enter into the model, is available for dREE complexation without any competition with dFe. Further model limitations are the fact that we only work with hypothetical numbers of strong organic complexes within the DOC pool without being able to identify which parts and structures of the DOC pool actually occur in the mesocosms and which complex dREEs. The two different approaches, ‘High-DOC’ and ‘Low-DOC’, do however show that the concentration of DOC only influences the proportions of REE-DOC, but not the variations within the REE series ( Figure 5 ). Schijf et al. (2015) stated, that the effect of DFOB on REEs and the associated fractionation is likely to be typical for a large variety of organic ligands. Therefore, regardless of the true structure of the organic ligands, the modelled organic complexes can display a natural scenario. We are, however, aware of the fact that DFOB is not able to mirror the complexity of the DOC pool and that some organic ligands could have different effects on the REEs. Lanthanide-binding proteins excreted from bacteria that need LREEs to grow seem, in contrast to the organic ligands described above (Cotruvo, 2019), to preferentially complex LREEs. These are, however, as stated above (section 5.2.4) likely to be of minor importance given the small niche of bacteria producing these lanthanophores (Huang et al., 2019). The model is an important first step to assess the influence of dissolved organic ligands within the DOC pool on REEs without a claim to quantitative conclusions.

The threshold of 200 µmol L^-1 DOC is above the mean DOC concentration in the southern North Sea of around 130 µmol L^-1 (Van Engeland et al., 2010). It seems that only in highly productive waters either at specific locations or days during the year, as for instance during phytoplankton spring blooms, as presented in this study, the required ligand concentrations may be reached to cause significant changes in REE fractionation and temporal concentration patterns. It is important to notice, that the DOC threshold presented here is not a fixed number, which can be directly applied to any natural system, but it is rather specific to the studied mesocosm approach. In general, the DOC concentrations are more likely to act in a continuum, with more DOC available, leading to an increase in dREE complexation. If, however, REEs compete with other metals like Fe for the siderophores, dREE complexation might only be possible if organic ligands are available in excess to dFe (Schijf et al., 2015). Thus, a threshold of DOC, or rather a threshold of organic complexing agents, needs to be exceeded in order to significantly complex dREEs and outcompete inorganic REE binding ligands. This threshold then would depend on trace metal concentrations other than REEs, as well as structure and complexity of the DOC pool. To further address this issue, future modelling approaches should include a molecular analysis of the DOC pool to identify the potential REE-ligands within this pool and their temporal variations and to predict their influence on the respective REE concentration patterns in seawater.