In this step, the spectra were first calibrated, normalized by median peak intensity, aligned onto common mass-to-charge (m/z) ratios detected by bin-wise KDE, and cleaned by a two-step peak picking method to extract robust and informative signals in the spectra, including a peak prominence filter and a geochemical-context-based filter based on the varying deposits visualized by the X-Ray image of the sediment.

The performances of NMF at varying ranks between 3 and 20 are compared in Supplementary Text S2. As a result, an NMF rank of 10 was chosen for both dataset A and dataset B, resulting in the extraction of in total 17 and 14 distinct spatial signatures, respectively, over 30 NMF runs (Supplementary Figures S6, S7). From these, 5 and 10 spatial signatures (Figure 3) are stable because they were reached in at least 27 out of 30 (i.e., >90 percent) NMF runs. The unstable signatures likely represent the local outliers among these stable signatures, but their differences are too trivial to be stably separated. For example, in dataset B, the unstable signature B-NMF12 is in fact similar to the stable signature B-NMF6. Although they are not well-separated from the rank selected in this study, they are in fact separated from each other in the consensus-matrix-based molecular network.

In the following, we compared the spatial signatures with the varved structures of the sediment revealed by its X-ray photograph, and their quarterly averages after applying an age model with the environmental parameters determined in the water column to show that these mathematically derived signatures are indeed relevant to biogeochemical processes. However, one should keep in mind that the inferred ecological and geochemical implications hereafter are speculative and could be oversimplified due to the numerous stereochemical possibilities of the molecular formulas detected, along with other factors. The spatial signatures are coded in the following way: assuming dataset V ∈ { A ,   B } , given V ∼ W H for dataset V , we refer to the columns in W as signatures and the k t h column as a signature V -NMFk, whose ion images are denoted by its entries. The contribution of each molecule to a signature V -NMFk is denoted by the entries in the k t h row of H . Although the unstable signatures are worth investigating, for example, B-NMF10 and B-NMF12 show almost complementary ion images (Supplementary Figure S7) that likely indicate distinct sources, only the stable signatures are discussed in the following for succinctly demonstrating the capability of MSI in conjunction with untargeted data mining in extracting biogeochemically relevant molecular signatures from the sediments. For stable signature V -NMFk, we show

(1) its ion image that represents the cluster center (the weighted average of all ion images in the cluster; Figure 3, left), constructed using the entries in the k t h column of W ;

Among all signatures discovered in the two datasets, A-NMF3 shows the most distinct laminations (Figure 3). Figure 6 further shows that the light laminae of A-NMF3 match well with the light diatom ooze laminae of the SBB sediments, which result from enhanced primary production (Schimmelmann and Lange, 1996; Bull et al., 2000). The corresponding molecular network in Figure 4 indicates that A-NMF3 consists of chlorin-like compounds, and the most abundant molecule is putatively pyropheophorbide a. We speculated that pyropheophorbide a and the other chlorin-like compounds in A-NMF3 derive from zooplankton and zoobenthos grazing on diatoms (Head et al., 1994; Szymczak-Żyła et al., 2011), and the distinct light/dark laminae of A-NMF3 may result from the opportunistic bloom-and-bust life cycle of diatoms (Butterfield, 1997). In fact, by comparing the quarterly averaged abundance of A-NMF3 to environmental parameters (Figure 5), we found that the abundance of A-NMF3 is positively correlated with salinity and nutrient concentrations including the concentrations of nitrite, phosphate, and nitrate, while negatively correlated with bottom water oxygen concentrations. The positive correlation could indicate the increase in diatom populations and consequently A-NMF3 associated with increases in nutrient supply (Lange et al., 1997), while the negative correlation could suggest the depletion of oxygen in bottom waters, resulting from the subsequent organic matter remineralization (Reimers et al., 1990; Alfken et al., 2021).

A-NMF6 consists of alkenones produced by the ubiquitous coccolithophores (Volkman et al., 1980b, 1980a; Marlowe et al., 1984a, Marlowe et al., 1984b), i.e., C₃₇- and C₃₈-, di- and tri-unsaturated alkenones, together with their corresponding isotopologues and a compound with a tentative chemical formula of C₃₇H₇₀O₂, possibly a coccolith derived di-unsaturated alkenoate, e.g., methyl hexatriacontadienoate (Marlowe et al., 1984a; Figure 4). Compared to A-NMF3, A-NMF6 shows less distinct laminations (Figure 6), and its light laminae do not match very well with each other. This agrees with previous observations that coccolithophores usually do not bloom along with diatoms (Zhao et al., 2000). In addition, the quarterly averaged A-NMF6 has a slightly different correlation with environmental parameters compared to A-NMF3 (Figure 5). It is positively correlated with nutrient concentrations, including the concentrations of silicate, phosphate, and nitrate, while negatively correlated with the temperature of the water column. Moreover, although both A-NMF6 and A-NMF3 are negatively correlated with oxygen concentrations in the water column, A-NMF6 is best matched with the oxygen concentrations in the shallower water column, while A-NMF3 is best matched with the oxygen concentrations in the bottom water column.

Although B-NMF6 also shows distinct laminations in its ion image, in contrast to A-NMF3 and A-NMF6, it is more likely to be a terrigenous signal as its light laminae match with the dark laminae of the sediment (Figure 6) formed by seasonal runoff during heavy winter rains (Schimmelmann and Lange, 1996; Bull et al., 2000). MSI also reveals scattered hotspots in B-NMF6 in addition to the laminae couplets, which could indicate the presence of terrigenous particles. The corresponding molecular network of B-NMF2 shows that it consists of steranes and unspecified oxygen-containing compounds (Figure 4). As steranes are not formed from autochthonous organic matter in immature sediments, they are likely derived from eroded source rocks of the Monterey Formation or related oil seeps (Hinrichs et al., 1995); the oxygen-containing compounds with even-numbered carbons, i.e., C₂₄H₄₈O₂, C₂₆H₅₂O₂, and C₂₈H₅₆O₂ could be derived from higher plants (Ratnayake et al., 2005; Kusch et al., 2010). Moreover, the quarterly averaged B-NMF6 is positively correlated with the bottom water oxygen concentrations, while negatively correlated to salinity and nutrient concentrations (Figure 5), which is consistent with its assignment to wintertime terrigenous input during periods of weakened upwelling.

The ion image of B-NMF2 is rather complicated: it is ubiquitous in the sediment, while some of its light laminae (more abundant) match quite well with the dark laminae of the sediment (Figure 6). Hierarchical clustering of the most representative compounds in B-NMF2 displayed in Supplementary Figure S9 further reveals the internal structure of the molecular cluster, which consists of sterols and compounds that have one to three oxygen atoms and double bond equivalents typically found in steroids. Sterols and steroids are known to be produced by diverse species of eukaryotes such as zoo- and phytoplankton and terrestrial plants (Volkman et al., 1998; Volkman, 2003). The mechanism causing the colocalization of these compounds in the sediment deserves further investigation. It is possible that B-NMF2 captures the superimposed signal of steroid-like compounds from varying sources. For example, one of the most abundant molecules in B-NMF2, C₃₀H₅₂O, can be assigned to dinosterol, which is mainly produced by dinoflagellates (Volkman et al., 1993), the algae which often blooms during relaxed upwelling (Kennedy and Brassell, 1992; Smayda and Trainer, 2010). However, the formula is also consistent with terrigenous biomarkers such as lupanol (Pancost et al., 2002) or the bacterial biomarker hopanol (Pearson et al., 2001). Alternatively, B-NMF2 could reflect benthic biogeochemical processes because of its strong positive correlation to bottom water oxygen concentration and strong negative correlation to bottom water salinity and phosphate concentrations (Figure 5).

The other signatures exhibit distinctive spatial patterns as well, although not all chemical formulas of the associated molecules can be assigned with certainty, partly due to the absence of isotopologue peaks. For example, signatures A-NMF5 and B-NMF7 are both characterized by the “intrusion” structure in the uppermost ∼5 cm sediment in addition to the dark and light laminae (Figure 3); such a structure is also visible on the sediment slide. A-NMF7 is omnipresent in the sediment, with distinct laminations only visible in the topmost sediment, and it consists of a large number of molecules, with the most abundant one likely to be bacteriohopanetetrol (Figure 4). Quarterly averaged A-NMF7 is positively correlated with silicate concentration but negatively correlated with temperature and nitrite concentration (Figure 5). A-NMF9 is of low abundance in the topmost surface sediment, and it is associated with three compounds that are apparently homologs, together with their ¹³C isotopologues. Notably, A-NMF9 does not show a strong correlation with any environmental parameters. In addition, B-NMF1 is characterized by curved dark/light boundaries in the sediment (Figure 3). We speculated that they are polyethylene glycol contaminants (PEG; H⁺ and Na⁺ adducts) because of the characteristic mass difference of 44 Da between individual molecules (Figure 4). Moreover, it should be noted that all the signatures discovered in this study spatially overlap with each other to some extent, causing image segmentation to be difficult. This is due to either a blurred zonation of biomarkers or the relatively small mass window analyzed, which limits the number of diagnostic biomarkers.

The proposed data processing workflow in this study is explorative and mostly unsupervised as except in the geochemical-context-based peak picking step, it only utilizes the information carried by sedimentary MSI datasets themselves. In order to demonstrate the potential of supervised learning approaches for novel molecular proxy discovery on the basis of data mining results, we trained linear models for indicating SST and sediment-water interface redox conditions, using the abundances of co-occurring molecules associated with clusters A-NMF6 (alkenones and derivatives) and B-NMF2 (incl. steroid-like compounds), respectively, based on the assumption that residuals among these co-occurring molecules are not random and result from biogeochemical influences such as changing environmental conditions or selectivity in biogeochemical processes. ElasticNet, a popular regularized linear regression algorithm implemented in the python package scikit-learn (Pedregosa et al., 2011), was employed to tackle multicollinearity (Gunst and Webster, 1975; Alin, 2010). Time-based cross-validation implemented in scikit-learn (Pedregosa et al., 2011) was employed to reduce overfitting and autocorrelation; whether the derived proxies are applicable in longer time series or at other sites is beyond the scope of this study.

Since it has been established that in dataset A, the U 37 K ′ index has a strong correlation with SST (Alfken et al., 2020), we first validated the use of ElasticNet on the alkenone cluster A-NMF6 and its relationship to SST. Although A-NMF6 contains four ¹³C isotopologues and an alkenoate in addition to the monoisotopic alkenones, the whole cluster is used for training without removing any peaks in order to show the reliability of ElasticNet with redundant variables. Figure 7A compares the SST measured in the water column and the SST indicated by the trained model and the U 37 K ′ index, suggesting that the model trained in this study indicates changes in SST. Supplementary Table 3 shows the resulting coefficients assigned to the nine molecules in A-NMF6 after cross-validated ElasticNet training and indicates that the two di-unsaturated alkenones and their isotopologues are positively correlated with SST, while the two tri-unsaturated alkenones and their isotopologues are negatively correlated with SST. This is in agreement with the definition of the U 37 K ′ index proposed in the previous study (Prahl and Wakeham, 1987). In addition, the putative alkenoate was assigned with a negative coefficient after training, although it weighed less than the two tri-unsaturated alkenones. This may suggest that the molecule is also regulated by coccolithophores to adapt to changing SST, which agrees with the observation by Conte et al. (1998) that the cellular ratio of alkenoates to alkenones in coccolithophores decreases with increasing water temperature.

B-NMF2 consists of diverse steroid-like compounds with one or more oxygen-containing functional groups, which might be sensitive to redox conditions. We further applied ElasticNet to train a model for indicating bottom water redox conditions using the abundances of co-localized molecules in B-NMF2 against bottom water oxygen concentrations. Figure 7B compares the bottom water oxygen concentrations with the trained model and the redox-sensitive C₂₉ stanol/stenol ratio. In contrast to the model for indicating SST (cf. in Supplementary Table S3), the resulting coefficients for indicating redox conditions are all equal or greater than 0, with smaller absolute values (Supplementary Table S4). The non-negative coefficients of the compounds suggest that variations in the model could be driven by 1) the relative abundances amongst these steroid-like compounds in B-NMF2 or 2) the relative abundance of these steroid-like compounds as a group compared to all other compounds in the spectrum. Although all non-negative coefficients could make the model less favorable as a paleoenvironmental proxy and more sophisticated nonlinear models are likely required to characterize the link between these compounds and the redox condition, such rather complicated supervised learning approaches are beyond the scope of this study and will be evaluated in our future work. Since standard scaling was applied before model training, the magnitudes of resulting coefficients depend largely on the number of dependent variables. The smaller absolute coefficients assigned in B-NMF2 compared to those assigned in A-NMF6 results from the much larger numbers of variables (molecules) in B-NMF2. Nevertheless, the resulting coefficients agree with the observation in the aforementioned section that B-NMF2 is positively correlated with bottom water oxygen concentrations, and the ElasticNet model, here, further extracted the most redox-related biomarkers from B-NMF2. It is interesting that the top three molecules are all C₃₀ pentacyclic triterpenoids with varying unsaturation degrees. This may suggest that they are the most recalcitrant components of B-NMF2 in the sediment. However, other factors cannot be ruled out, for example, it is also possible that the three putative C₃₀ pentacyclic triterpenoids have the strongest contribution from terrestrial input during the winter season as compounds with the same chemical formulas have been reported in higher plants (e.g., Escobedo et al., 2012; Kennedy, 2012).