Tracking the spatiotemporal variations of statistically independent components involving enrichment of rare-earth elements in deep-sea sediments

Deep-sea sediments have attracted much attention as a promising resource for rare-earth elements and yttrium (REY). In this study, we show statistically independent components characterising REY-enrichment in the abyssal ocean that are decoded by Independent Component Analysis of a multi-elemental dataset of 3,968 bulk sediment samples from 101 sites in the Pacific and Indian oceans. This study for the first time reconstructs the spatiotemporal variations of the geochemical signatures, including hydrothermal, hydrogenous, and biogenic calcium phosphate components that were closely involved in the formation of REY-rich mud over the past 65 million years. An underlying key factor of significant REY-enrichment is a sufficiently low sedimentation rate that enables the mud to accumulate REY from seawater. In the early Cenozoic, a remarkably small supply of aeolian dust, compared with any other time and region, facilitated the deposition of very high-grade REY-rich mud in the South Pacific. This indicates an important link between the genesis of the seafloor mineral resources and Earth’s dynamic phenomena such as climate change and plate tectonics.


Mathematical formulation of the Independent Component Analysis model
In the literature of Independent Component Analysis (ICA; Ref. S1), the mathematical model is usually formulated as where x is the column vector whose i-th element is an observed random variable x i (i = 1, 2, 3, ..., m, where m is the number of observed variables), s is the column vector whose j-th element is an independent random variable s j (i.e. independent component, IC; j = 1, 2, 3, ..., r, where r is the number of ICs), and A is the m × r mixing matrix. In general, m and r may not be equal, and r should not be larger than m. Equation (S1) indicates that the observed data can be expressed by a linear combination of the independent source components, or the compositional data can be plotted in subspaces spanned by the new base vector s.
Practically, however, x contains the index of the time or sample (Ref. S1) and thus can be written as x i (t) (t = 1, 2, 3..., n, where n is the number of samples). Therefore, in our problem, x is equivalent to the m × n observed data matrix whose rows and columns correspond to the m variables and n samples, respectively, which can be denoted by the transposed matrix X T . The superscript T denotes the transpose, and X is the same as the n × m observed data matrix defined in 'Fundamentals of ICA' in the Methods section. In the same manner, s can be written as the r × n independent source matrix S T because t is the common index between x and s (Ref. S1). In addition, if we arbitrarily rewrite the mixing matrix A in equation (S1) as A T , where A T is also an m × r matrix because this is a merely notational change, the formulation of equation (S1) can be written as

X T = A T S T (S2)
By transposing both sides of equation (S2), we obtain where X is the n × m observed data matrix, S is the n × r independent source matrix, and A is the r × m mixing matrix (Refs. S2,S3).
In the open R package that we employed to perform ICA (Ref. S2), the observed data matrix X is assumed to have the samples in its rows and the variables in its columns (i.e. an n × m matrix). Indeed, actual data are often organized as a matrix whose rows correspond to individual samples, and the ICA model for such data was represented by the same equation as (S3) (e.g. Ref. S3). Therefore, for the convenience of the potential users of the package based on this paper, we denoted our ICA model by equation (S3), as discussed in 'Fundamentals of ICA' in the Methods section.

Definitions of independent component scores and independent component loadings
We defined IC scores here as the coordinate values of the sample data in the IC space, which represent the intensities of each independent source signal in each sample. The IC score of the j-th element in the i-th sample corresponds to the (i, j) component of the estimated independent source matrix S.
We defined IC loadings here as the relative contributions of each element to the individual ICs, which correspond to the slope of each IC vector projected in the compositional space. In our ICA model, as given in the preceding section, the loadings of the i-th IC correspond to the i-th row of the estimated mixing matrix A. However, the absolute values of the contents and their variations differed significantly among elements. The CaO content varied from 0 to ~56%, whereas the Ce content varied from 0 to ~660 ppm, which makes it difficult to evaluate the contribution of each element to each IC at a glance. Therefore, in the present paper, to clearly visualize the influence of each element to each IC, we show the relative loading values rather than the raw loading values ( Supplementary Fig. S6). Because the ICs do not have specific norms, the values of the relative loadings can be expressed in arbitrary units. A relative loading of the j-th element for the i-th IC can be obtained by dividing a ij , or the (i, j) component of A, by the sample standard deviation of the j-th element.

The meaning of positive and negative signs in Independent Component Analysis results
The positive and negative loadings represent the elements that increase and decrease accompanying the unit increase of each IC score. In fact, ICA contains the ambiguity of the sign of each IC. In the ICA model, because both the independent source matrix S and the mixing matrix A are estimated at the same time, we can multiply the i-th independent component (i.e. the i-th column of S) and the corresponding i-th row of A by -1 simultaneously without affecting the analysis (Refs. S1 and S4). Therefore, we can arbitrarily choose the signs of the IC scores and loadings. Here, for an intuitive description of the geochemical ICs, we adjusted the signs of the IC scores and loadings as the characteristic elements increase when the IC score positively increases. This adjustment never affects the essential meaning of the ICA result.
Checking the statistical robustness of the independent components

Effect of outliers in the dataset
In this study, to minimize the undesirable influence of outliers or specific sites/lithologies/samples, we chose the exponential function as the evaluation function of non-Gaussianity (i.e. G(y) in 'Fundamentals of ICA' in the Methods section) because it reflects the characteristic structure around the centre (mean) of the data more intensely than other proposed functions such as kurtosis, and it theoretically gives the most robust results theoretically (Refs. S1 and S5).
Here, we checked the robustness of our ICA result by using a limited dataset, excluding the samples with high values for each IC (greater than 3) in the original result. Because each IC is constrained to have unit variance in the fast ICA algorithm (Ref. S1), an IC score of 1 corresponds to 1 standard deviation (1σ) along each IC. Thus, samples having IC scores greater than 3 constitute the long tails outside the 3σ range in the IC score distributions that have strong non-Gaussianity such as large skewness and kurtosis ( Supplementary Fig. S5). IC3 ranges from -1.7 to 1.3 and does not show a long-tailed distribution like other ICs ( Supplementary Fig. S5). Therefore, no outlier was considered for IC3.
The ICA result with no outliers (n = 3,492) was fundamentally the same as the original result ( Supplementary Fig. S8 and S9). To confirm this, we further calculated the projection of the excluded 'outlier' samples into the IC spaces by using the whitening matrix and the un-mixing matrix estimated from the dataset excluding high-IC samples ( Supplementary Fig.   S10). The projected data distributions of all samples were almost identical to the original ICA result derived from all samples ( Supplementary Fig. S4), although only IC6 showed a slightly oblique feature compared with the original result. In general, this test confirmed that the geochemical ICs derived from all data certainly reflect the fundamental structure in the dataset and are not artefacts generated by a small number of outliers.

Effect of the difference in sample subsets between the Pacific and Indian oceans
Here, we evaluated the change in ICA results change when each basin is treated as separate populations in the ICA. The ICA result using only the Pacific Ocean samples extracted almost the same components as the original result (Supplementary Figs. S11 and S12), although the loading values changed slightly. Thus, the slopes of IC vectors projected in the compositional space differed slightly from the original result.
On the contrary, the ICA result using only the Indian Ocean samples showed partly different features (Supplementary Figs. S13 and S14). The ICs corresponding to the original IC1, IC3, IC5, IC6, and IC7 were extracted, although the loading of Mn in IC7 showed a somewhat smaller value than those in the hydrothermal IC in the original and Pacific Ocean results. The IC4 showing Ca-phosphate-enrichment was not extracted from the Indian Ocean subset because the rare-earth elements and yttrium (REY)-rich mud samples in the Indian Ocean are predominantly IC1-type mud in our dataset. Alternatively, another IC assigned as new IC4 in the Indian Ocean result was extracted which showed a concurrent increase in Si, Ti, Al, Fe, Mg, and K ( Supplementary Fig. S14). This IC signal appeared to be strong in the sites near the continent (e.g. Site 241 near the eastern coast of Africa, Site 223 off the Arabian peninsula, and Site 1171 south of Australia; Fig. 1 in the main text). This finding suggests that the IC corresponds to a terrigenous detrital component, although its contribution appears to be very minor in the entire dataset. The remaining ambiguous IC was assigned as IC2.
Although partly different components were extracted, the main structures of the data appeared to be common to data subsets from the two ocean basins. Thus, we discussed the general mechanism of REY-enrichment on the basis of the entire sample set in the main text.

Effect of the recalculation of raw data of major elements
In our dataset, the major element contents were determined by X-ray fluorescence (XRF) analysis, as discussed in the Methods section. Then, the total contents of the major elements (i.e. SiO 2 , TiO 2 , Al 2 O 3 , Fe 2 O 3 , MnO, MgO, CaO, Na 2 O, K 2 O, and P 2 O 5 ) and loss on ignition (LOI) were normalized to 100 wt.% (Supplementary Data S1 and S2), although Na 2 O and LOI were not included in the dataset analysed by ICA because of the effect of sea salt and the ambiguity of the origin, respectively.
To check the influence of the recalculation of raw data of the major elements, we performed ICA on the dataset by using XRF raw data of the nine major elements before normalizing to 100%, including LOI (Supplementary Figs. S15 and S16). The test resulted in essentially the same figures as the original ICA result discussed in the main text (Supplementary Figs. S4 and S6).

Additional ICA test on a sample subset without significant dilutions by using selective variables
We performed ICA on the sample subset, excluding high-Ca (CaO > 10 wt.%) and high-Si (SiO 2 > 70 wt.%), with the selective elements of Fe 2 O 3 , MnO, CaO, and P 2 O 5 and all of the separated REY elements. As described below, the additional ICA resulted in extraction of three key ICs. The new IC1, IC2, and IC3 corresponded well to the original IC7, IC4, and IC1, respectively. It should be noted that even if the number of extracted ICs was increased, the new three key ICs were always extracted and essentially never changed. Therefore, we consider that these are inherent features of REY-rich mud characterized by the original IC1, IC4, and IC7.
The result of the additional ICA is shown in Supplementary Figure S17, which indicates the relative loadings of each IC in the cases in which the numbers of extracted ICs were 3 and 4. These results led us to consider that the new IC3 corresponds to the original IC1.
When the number of ICs was increased to 4, the above three ICs were extracted again with no essential changes. This confirms that these three ICs extract the fundamental features inherent in the data structure. The additionally extracted IC4 has relatively large loadings of Fe 2 O 3 and MnO, suggesting the influence of the hydrothermal component. This IC also largely contributes to P 2 O 5 content, which implies the influence of biogenic Ca-phosphate. In addition, the relatively large Y loading compared with the other REEs, together with negative Ce loading, is characteristic of this new IC4. The samples with high IC scores regarding the new IC4 are from Sites 75, 319, and 597 in the eastern South Pacific, which showed intense signals in both original IC4 (i.e. biogenic Ca-phosphate) and IC7 (i.e. hydrothermal). The additional IC4 might reflect the influence of a relatively minor effect superimposed on the fundamental features, such as a diagenetic process. However, for more detailed and reliable discussion, various trace elements should be incorporated in the dataset.

Remeasurement of major element contents in previous data by X-ray fluorescence
In our inductively coupled plasma mass spectrometry (ICP-MS) analysis, bulk sediment samples were dissolved by HNO 3 -HF-HClO 4 digestion, and the sample solutions were finally evaporated until dry. During this process, silicate minerals were digested by HF, and Si was lost as SiF 4 through evaporation. Therefore, ICP-MS analysis of acid-digested samples cannot determine the Si content. In addition, the LOI content was not analysed. The SiO 2 content in the previously published data by Ref. S6 was estimated by using the analysed contents of other major element oxides and was assumed to have a constant value of H 2 O + CO 2 . In this study, to construct a complete dataset of uniform quality, we reanalysed all of the major elemental contents by using XRF, and we measured the LOI for all samples used in the ICA. pelagic sediments independently, resulting in a spurious correlation between REY content and phillipsite abundance (Refs. S9 and S10). Indeed, typical REY-rich mud with relatively high REY content is generally described as zeolitic clay; thus, the REY-enriched IC that showed an Al-rich trend in the previous study should have extracted the composite feature of pelagic clay abundant in phillipsite that also contains other aluminosilicates, Fe-Mn (oxyhydr)oxides, and REY-enriched biogenic Ca-phosphate.

Age models
Age models were constructed for each site on the basis of previously published biostratigraphic data. We assigned numerical ages to our samples by linear interpolation between the closest age-diagnostic events such as first/last occurrence of key fossils (e.g. calcareous nannofossils, foraminifera, and radiolaria) by using the geologic time scale (Ref.

Reconstruction of spatiotemporal distribution of independent component scores
Reconstruction of palaeogeography and the tracks of each site from the present day to 65 Ma was achieved by using GPlates (Refs. S15 and S16) in steps of 5 Myr. On the palaeomap of each time slice, the IC scores were plotted on the palaeopositions of each site as age-weighted mean values. We calculated the weighted mean of IC scores s j at time slice T i as where t i is an assigned age of sample j, d j is the age deviation between the time slice i and sample j, w j is the weight of sample j, and s j is the IC score of sample j. Equations (S4) to (S6) indicate that the closer in age a sample is to each time slice, the larger the weight of the sample.    Note that the reference materials were not included in a data set analyzed by ICA, but were merely projected into IC spaces by applying the same linear transformation as that estimated from samples by ICA. The approximate IC scores of the reference materials S rm are calculated by multiplying their compositional data matrix X rm by the whitening matrix K and the un-mixing matrix W estimated from the samples (i.e. S rm = X rm KW). Supplementary Figure S5. Frequency distribution of the seven ICs extracted in this study. The horizontal axis corresponds to each IC score. The vertical solid and dashed lines represent mean and median values for each IC, respectively. Because each IC is constrained to have unit variance in the fastICA algorithm, an IC score of 1 corresponds to 1 standard deviation along each IC. A Gaussian distribution of zero mean and unit variance is also shown for comparison.

Relative loadings (arbitrary unit)
Supplementary Figure S6. Relative loadings of each element for IC1 to IC7. The relative loading of the j-th element for the i-th IC can be obtained by dividing a ij , or the (i, j) component of the mixing matrix A, by a sample standard deviation of the j-th element.  Fig. 2a. (b) Results of 100 computational runs for evaluating the effects of random initial values given in the fastICA algorithm to estimate the matrix W. All ICs showed well-converged slopes, indicating that the ICs have robustness to built-in random valuables of the algorithm. (c) Results of other 100 runs using synthetic perturbed datasets to evaluate the effects of data uncertainties such as analytical error. The slopes of ICs still converge well except for IC2 (dark blue), which shows slopes that different from b and much less convergence than other ICs. This suggests that IC2 is not robust to data uncertainty.  Figure S5). Therefore, no outlier is considered for IC3.

Relative loadings (arbitrary unit)
Supplementary Figure S16. Relative loadings of each element for IC1 to IC7 in the result of ICA using XRF raw data of nine major elements before being normalized to 100% including Na 2 O content and loss on ignition.  Figure S19. Flux of Nd required to form REY-rich mud as a function of the sedimentation rate. The assumed water depth is 5,000 m, which is the typical depth for the occurrence of REY-rich mud.
The slopes of the lines changed with the bulk ΣREY content. The horizontal blue bars indicate available Nd supplies from the overlying seawater. The pink areas represent REY-rich mud with relatively high resource potential (bulk ΣREY = 1,000-2,000 ppm).