Distribution of Butyltin Compounds in the Coastal Environment of the Bahía Blanca Estuary, Argentina

This study evaluates for the first time the distribution and accumulation of butyltin compounds (BTs) in different compartments such as seawater, sediments, suspended particulate matter (SPM), and mussels (Brachidontes rodriguezii) in the Bahía Blanca estuary. The samples were collected from six sampling sites with different anthropogenic impacts. A better visualization and interpretation of data was achieved using chemometric tools (Tucker4 model), which made it possible to reveal the main relationships among the variables. This analysis showed the presence of BTs in all the estuarine environmental compartments, even in sites with low human intervention. The relationships found among BTs levels, seasons, and environmental matrices show the importance of biological processes such as phytoplankton blooms and remobilization of sediments (by tidal dynamics and/or periodic dredging) in BTs distribution and degradation. In addition, partition coefficients showed that mussels mainly bioaccumulate tributyltin from sediment, water and, to a lesser extent, SPM.

In the past, organotin compounds (OTCs) were traditionally used in a wide range of applications (Chen et al. 2017). However, in the mid-70 s, OTCs gained notoriety due to the outstanding biocidal properties of triorganotin derivatives. Tributyltin (TBT), in particular, has been used as an active component in antifouling paint formulations to prevent biofilm formation on marine structures and ship hulls (van Gessellen et al. 2018).
TBT reaches the marine environment by leaching from these paints (Okoro et al. 2016). Once in the environment, TBT and its degradation products-dibutyltin (DBT), monobutyltin (MBT) and inorganic tin-can be involved in different processes, such as diffusion, sorption, lixiviation, sedimentation, biodegradation, and resuspension (dos Santos et al. 2016). Degradation of TBT may occur by photolysis, chemical and thermal cleavage, or through microbiological activity (El Hadj et al. 2016). In turn, butyltin (BT) bioavailability is influenced by numerous environmental factors, physicochemical variables, seasonal variations, and geographic conditions (Fang et al. 2012). Moreover, BTs can be found adsorbed to sediment and particulate matter, or dissolved in both pore water and water column. It can even be readily incorporated into the tissues of filter-feeding zooplankton, grazing invertebrates and, eventually, into higher organisms such as fish, water birds, and mammals (Ohji et al. 2007).

3
Owing to the ecotoxicological effects of TBT (Filipkowska and Lubecki 2016;Martínez et al. 2017), the International Maritime Organization (IMO 2001) adopted the International Convention on the Control of Harmful Antifouling Systems, in force since September 17, 2008, which banned the use of TBT as an active component in antifouling systems for ships around the world. However, until now, this restriction has not been implemented for DBT and MBT (Fang et al. 2017). In spite of this regulation, in recent years, a large number of studies have reported that TBT is still present in the environment (Abraham et al. 2017;Batista et al. 2016;Deng et al. 2018;Maciel et al. 2018). This is particularly important in estuaries with restricted circulation hosting harbors, marinas, and shipyards (Castro and Fillmann 2012). This is the case of the Bahía Blanca estuary (BBE), an area under increasing anthropogenic pressure and the second largest estuary in Argentina, including one of the most important industrial hubs in South America.
Considering the possible incidence of BTs in the whole ecosystem, it is necessary to monitor these compounds in different environmental matrices to assess their distribution and bioavailability (dos Santos et al. 2016). In this sense, the aim of the present study was to analyze the levels of TBT, DBT, and MBT in several samples from different environmental compartments (sediment, suspended particulate matter (SPM), water column, and mussel tissues) that were seasonally collected (temporal dimension) from several sampling points (spatial dimension) in the BBE. Data visualization and interpretation were carried out by chemometric multivariate analysis, particularly using the Tucker4 model (Bro 1998;Marinho et al. 2019). The current study is the first comprehensive analysis performed on BTs in different environmental matrices with the assistance of N-way chemometric tools.

Study Area
The Bahía Blanca estuary is located in the southwest of Buenos Aires Province, Argentina (38°45′-39°40′ S and 61°45′-62°20′ W) (Fig. 1). The estuarine area comprises a dense arrangement of meandering channels and islands surrounded by extensive intertidal mudflats and marshes (Dutto et al. 2017). Approximately 85% (250,000 ha) of the total area of this system is under legal protection through four ecological reserves (Bravo et al. 2018). Despite this, Fig. 1 Map showing the study area and sampling points. S1: Rosales Harbor (ship-repair yard); S2 and S3: Ingeniero White Harbor (S2: near thermoelectric plant, S3: fishing and recreational zone); S4: Gal-ván Harbor (industrial harbor, gas, and oil loading buoy); S5: Cuatreros Harbor (fishing and recreational port); and S6: Villarino Viejo (rural zone) the BBE undergoes strong anthropogenic pressure from urban and industrial activities on its northern boundary and also houses the most important deep-water port system in Argentina. In particular, deepening dredging activities carried out in the Main Navigation Channel promote resuspension, transport, and redistribution of sediments (Grecco et al. 2011).
Six sites, affected by different intensities of human activities, were selected from the outer to the inner zone of the estuary (Fig. 1). Sampling station S1 is located in Rosales Harbor, a coastal town in the middle of the Main Navigation Channel. This place, frequented by oil tankers and port service vessels, is located near Belgrano Harbor, the main military port in Argentina. This is where the hulls of many commercial and military ships are continuously repaired and painted. S2 is located in Ingeniero White Harbor, near a thermoelectric generation plant facility. This site is one of the major commercial and general cargo ports in Argentina, where fishing as well as leisure activities normally take place. The S3 site, also located in Ingeniero White Harbor, is characterized by the presence of numerous fishing boats. S4 lies near Galvan Harbor, in the vicinity of a petroleum derivatives loading buoy. Cuatreros Harbor (S5) is located in the inner part of the estuary, in a small recreational/fishing port. Finally, the S6 sampling point (Villarino Viejo) is at the head of the estuary, a rural area with a relatively low anthropogenic impact.

Samples
Sampling was performed every three months in 2014 on board the IADO IV research vessel, at low tide. Sediment samples were collected from the superficial layer (1-5 cm) using a stainless steel grab sampler and an acrylic corer. Native mussels (Brachidontes rodriguezii) were collected, according to their availability, from natural banks, dock columns, platforms, surface sediments, or surface water (0-1 m). On average, sixty mussels from each sampling point were used to prepare composite samples. Next, approximately 4 L of seawater was collected in amber glass bottles at each sampling location. Immediately after, the samples were refrigerated in iceboxes, stored in containers to avoid exposure to light, and rapidly transported to the laboratory. Afterward, the water samples were filtered through a polycarbonate membrane with a pore size of 0.4 μm (Mil-lipore® HTTP 04700) using a vacuum pump to prevent BTs adsorbed onto SPM from interfering with the analytical determination of these pollutants in the water samples (Berg et al. 2001;Michel et al. 2001). Thus, both filtered estuarine water samples and SPM retained in the filter were obtained. Before using them, the membranes were treated with 0.7% nitric acid, rinsed with deionized water, and dried at 50 °C to a constant weight. The filters with the retained material, sediments, and mussel samples were frozen at − 18°C, lyophilized, and stored at 4°C until analysis. The water samples were stored in the dark in amber glass bottles at 4°C until further use.

BTs Determination
At the beginning of the procedure, tripropyltin chloride (TPrT) was added as a surrogate recovery standard together with the sample and ultrapure deionized water in a centrifuge tube. Next, BTs determination in seawater samples was performed by liquid-liquid extraction and derivatization with sodium borohydride in hexane (Quintas et al. 2019). For sediments, SPM, and mussels, the freeze-dried samples were placed in a methanolic NaOH solution under ultrasound energy (Quintas et al. 2018). The same derivatization reaction was then carried out by adding sodium borohydride to the hexane. In all cases, one microliter of the hexane layer was injected in splitless mode (220 °C) into a gas chromatograph (Agilent 7890 B) coupled to a mass detector (Agilent 5977A). Before injection into the gas chromatograph, the tetrabutyltin (TeBT) standard solution was added as an internal standard.
The GC column was an HP-5MS fused silica column (30 m; 0.25 mm i.d.; 0.25 μm film thickness) with helium as a carrier gas (flow rate of 1 mL min −1 ). The mass spectrometer was operated in the electron impact mode (70 eV). The temperature program was 55 °C for 1 min, followed by 20 °C min −1 up to 200 °C, maintaining that temperature for 5 min. Organotin hydrides were determined in the selected ion monitoring (SIM) mode. The selected ions were 119,121,177,179,233, and 235 m/z for TBT; 121, 177, and 179 m/z for DBT, and the same ions selected for DBT plus ion 119 m/z for MBT. For TeBT, the ions were the same as for TBT plus 289 and 291 m/z. TPrT was detected using ions 119,121,163,165,205, and 207 m/z. Each BT was confirmed by the retention time and relative abundance of ions against certified BTs standards. Chromatograms of the standard solution with MBT, DBT, TBT, TrPT, and TeBT and of the mussel and water sample extract are illustrated in Fig. S3 as examples.

Quality Control
The procedure was validated in terms of trueness and precision, using PACS-2 reference material for sediments. The other samples (i.e., water, SPM, and mussels) were spiked at three concentration levels for each BTs and then analyzed six times. The recoveries used for assessing trueness ranged between 82 and 97% for the reference material and between 94 and 100% for the other samples. Precision was evaluated under repeatability conditions and estimated as RSD (%) values. They varied between 1.1 and 1.9% for the 1 3 determinations of PACS-2 reference material and between 0.9 and 7.1% for the spiked samples. The results were in good agreement with sediment certified values (no statistical differences were obtained, α = 0.05) and with EURA-CHEM analytical validation recommendations, i.e., recovery of 70-120% and RSD (%) below 20% (EURACHEM 1998;Thompson et al. 2002). Limits of detection (LODs) were calculated considering S/N = 3 (signal-to-noise ratio) and taking the chromatogram baseline as an estimation of noise (Vidal et al. 2003). For sediment and SPM samples, the LODs obtained were 0.8, 1.1, and 3.5 ng Sn g −1 dry wt. for TBT, DBT, and MBT, respectively. The LODs for mussel samples were: TBT 1.2, DBT 1.8, and MBT 3.0 (ng Sn g −1 dry wt.). Finally, the LODs for water samples were 2.1, 5.4, and 8.4 ng Sn L −1 for TBT, DBT, and MBT, respectively.
To assess extraction efficiency, TeBT was used as a surrogate standard added to the samples at known concentrations. The samples were then subjected to the analytical procedure using TPrT as an internal standard. Thus, only the extraction step was evaluated because TeBT does not undergo the derivatization process. The recoveries obtained were 100%.

TBT Degradation
As mentioned above, TBT may undergo degradation via successive debutylation (Hoch 2001). This degradation can be evaluated through the butyltin degradation index (BDI) (Díez et al. 2002), which is calculated as the ratio between the concentration of TBT and the sum of its breakdown products, i.e., DBT and MBT: A BDI > 1 indicates recent inputs of TBT, while BDI < 1 implies that the parent compounds have been degraded and, therefore, the contamination should be considered 'old'. However, many factors could be responsible for the degradation of these compounds in the marine environment (Diez et al. 2005;Sousa et al. 2009), whose recentness is difficult to assess.
Nevertheless, in view of the fact that the samples were collected at the same time in a relatively small area under similar environmental conditions, using BDI seems to be adequate. The use of BDI is based on the assumption that the decomposition of TBT is the only source of DBT and MBT (Díez et al. 2002). The BDI data are shown in Supplementary Table S1.

Additional Analysis and Measurements
Temperature, pH, turbidity, and conductivity were measured at each sampling point. Salinity was estimated by conductivity and expressed as practical salinity units. The

BDI =
TBT DBT + MBT measurements, which are provided as Supplementary Material (Table S2), were performed in situ using a Horiba U-10 multi sensor probe. Sediment granulometry was analyzed by laser diffractometry (Malvern Mastersizer 2000), and results were expressed as the percentage of sediment particles below 63 µm in diameter (Supplementary Table S3).
Ignition loss was used to calculate organic matter as a percentage (%OM). A subsample of 10-20 g of dry sediment (105 °C until constant weight) was heated in a muffle furnace at 500 °C for 3 h (Commendatore and Esteves 2004). Particulate Organic Matter (POM) concentration was determined following the Strickland and Parsons method (1968) (range 10-4000 mg C m −3 ) using a Jenway 6715 UV-vis spectrophotometer. Before performing the analysis, ~ 250 mL of seawater was filtered through a muffled (450-500 °C, 1 h) glass fiber membrane (47 mm in diameter and 0.7 µm pore size) and frozen at − 20 °C in plastic bottles (APHA 1998;Grasshoff 1976). Both parameters are shown in Supplementary Table S4.
In addition, to evaluate the physiological status of mussels, the condition index (CI) was determined (Lau and Leung 2004) according to the following equation (Orban et al. 2002): This index may exhibit natural fluctuations as a result of food availability, spawning, reproduction, and the degree of stress of the organism, among others. The increase in CI represents a higher concentration of organic constituents in individuals, which is associated with growth. A decrease in CI, on the other hand, indicates a period of stress and implies a greater use of the organism's reserves. The condition index for the current study can be found in Fig. S1.

Data Treatment
The experimental data were arranged in a four-way array X (I × J × K × L). The values of some variables were beneath the detection limits. Therefore, they were substituted for a value equal to half of the detection limit for the corresponding variable (Álvarez et al. 2008). In this arrangement, the I-mode is related to the number of studied samples, the J-mode corresponds to BTs, the K-mode is related to the seasons of the year, and the L-mode to the different environmental compartments analyzed, i.e., sediments, mussels, SPM, and water. All data were scaled within the J-mode to allow BTs the same possibility of contributing to the model (Singh et al. 2007). Data array X was decomposed by applying a Tucker4 model (Bro 1998) using the N-way Toolbox (Bro and Anderson 2013) in Matlab 8.0 (The MathWorks).
Four loading matrices were obtained according to the following equation: where P, Q, R, and S are the number of factors in each mode, e ijkl are the elements of residual array E, and g pqrs are the elements of core array G (P × Q × R × S). This core array gives information about the interactions between the loading vectors in the different modes (Bro 1998). The optimal model was chosen on the basis of the residual sum of squares and degrees of freedom (named st-criterion) using the so-called convex hull method (Ceulemans and Kiers 2006;Stanimirova et al. 2009). The model in the convex hull with the largest st-coefficient was selected. The data were also analyzed using INFOSTAT software version 2008. The normal distribution was evaluated with the Shapiro-Wilk test (p > 0.05), while the Kruskal-Wallis test was used for the comparison of data from different sites and seasons. Correlations between BTs levels and physicochemical variables were evaluated using Spearman's test (p < 0.05 was regarded as significant).

Butyltin Levels and Physicochemical Variables
Butyltin concentrations are shown in Table 1. No statistically significant differences were found among the sampling sites (Kruskal Wallis, p > > 0.05) in all compartments, indicating a quite homogeneous distribution throughout the study area. This fact could reflect widespread diffusion due to the strong vertical mixing and tidal currents that traditionally led to defining the estuary as vertically homogeneous (Freije et al. 2008).
The data collected for BT concentrations could broadly be described as follows:

Sediments
TBT was quantified in 15 out of the 24 sediment samples. Levels ranged from < LOD (0.78) to 259 ng Sn g −1 dry wt. (summer, S3). Most of the samples in which TBT was not quantified belong to the winter and spring sampling. DBT was quantified in 21 out of 24 samples. Its concentration ranged from < LOD (1.08) to 782 ng Sn g −1 dry wt. (spring, S5). The lowest levels of DBT were found in winter (according to the Kruskal-Wallis test, winter significantly differed from the other seasons (p = 0.01), and the highest values a ip b jq c kr d ls g pqrs + e ijkl were observed in spring. MBT was quantified in 23 of the total samples, and its concentration ranged from < LOD (3.5) to 831 ng Sn g −1 dry wt. (autumn, S3). No statistically significant differences were found between seasons (Kruskal Wallis, p > > 0.05).
Regarding sediment grain size, all the sampling stations were classified as silty/clay. In general, the silt/clay ratios ranged from 0.77 to 1.60, growing from site S1 to the innermost part of the estuary. Even though BTs would be preferably adsorbed in the fine fraction (Hoch et al. 2003), no significant correlation was found between BT concentrations and grain size (p > 0.05). This lack of correlation can be attributed to the high variability of the sediment deposition dynamics caused by the high transport rate of sediments in the study area (Perillo and Piccolo 1991) and the remobilization of sediments that occurs as a consequence of the periodic dredging of the Main Navigation Channel. Finally, the %OM concentration ranged from 2.1 (S6, autumn) to 7.9% (S5, winter) (Table S4). Despite the fact that organic matter has been reported as the main adsorbing component of BTs (Berg et al. 2001), no correlation was found between OM% and BTs concentrations (p > 0.05) at the different sites or with respect to the seasons.

SPM
TBT and DBT were quantified in all SPM samples. TBT levels ranged from 27.4 (spring, S1) to 351 ng Sn L −1 (summer, S6), and DBT concentration ranged from 30.3 (spring, S1) to 461 ng Sn L −1 (summer, S6). On the other hand, MBT levels ranged from < LOD to 790 ng Sn L −1 (autumn, S1) and were below the LOD in all the samples collected in spring. POM was also studied, and the values were between 490 (S6, spring) and 1728 mg C m −3 (S1, summer) (Table S4). Given the affinity of MBT for organic matter (Hongxia et al. 1996), it makes sense to find low concentrations of this compound in SPM in the spring. Despite this fact, even though the content of organic matter plays an important role in BT destination, mobility, and availability in marine-coastal ecosystems, no significant correlations (p > 0.05) were found between BTs concentrations and POM.

Water Column
TBT was quantified in all the samples with concentrations ranging from 21.4 (summer, S6) to 387 ng Sn L −1 (autumn, S6). In turn, DBT was only quantified in spring and ranged between 260 (S5) and 510 ng Sn L −1 (S3). The MBT was between 289 (summer, S6) and 647 ng Sn L −1 (autumn, S2). We have recorded BTs in water samples of the Bahía Blanca estuary for the first time. As there is no previously published data from the study area, the concentrations from the present study were compared with data from other parts of the   (9) 297 (22) 374 (62) (13) 294 (6) world. Okoro et al. (2016) reported TBT levels of between 67 and 111,290 ng Sn L −1 in Cape Town harbor in South Africa, and Al-shatri et al. (2015) found TBT concentrations ranging from 140 to 1900 ng Sn L −1 in Saudi Arabia. In both research works, the maximum TBT levels described were much higher than those reported in this study. At the same time, the concentrations found in BBE were considerably higher than those reported on the Korean coast, with maximum TBT levels of 23.9 and 10.2 ng Sn L −1 (Kim et al. 2014;Lam et al. 2017) and in Taiwan by Liu et al. (2011) (28.8 ng Sn L −1 ).
The mean values of the physicochemical variables measured in situ in the water column are summarized in Table S2. All the parameters are within the previously reported ranges in the study area (Spetter et al. 2015). According to the Kruskal Wallis test, no significant differences were found among the sites (p > > 0.05), which would indicate that the variables evaluated were homogeneous regardless of the site, probably as a result of the mixing force of tidal currents. No correlation was found (p > 0.05) between the different physicochemical parameters evaluated -that is, particulate organic matter, temperature, turbidity, pH, and salinity-and the levels of TBT, DBT, and MBT found in the water column. Possibly, the presence of microorganisms, dredging activities, and characteristic tidal currents, among other factors, could hide a clear correlation between physicochemical variables and BTs levels (Basheeru et al. 2020;Simonetti et al. 2017).

Mussels (Brachidontes Rodriguezzi)
Brachidontes rodriguezii is a bivalve mollusk that belongs to the family of Mytilidae (d'Orbigny 1842). It is the dominant organism on middle and low intertidal rocky substrates on the warm-temperate shores of Argentina, acting as a key species in the ecosystem (Adami et al. 2008). This mussel species can be used as a bioindicator to monitor pollutants due to its wide geographical distribution, sessile lifestyle, easy sampling, tolerance to a considerable range of salinity, resistance to stress, and high accumulation of a wide range of chemicals (Baumard et al. 1998;Farrington and Tripp 1995). Furthermore, mussels are not able to depurate pollutants because they lack an efficient system to metabolize them (Burgeot et al. 2001), which underlines the fact that they can be used as good bioindicators. They also bioaccumulate pollutants throughout their lives. An analysis of their tissue will thus provide an indication of the bioavailable fraction of TB environmental contamination.
TBT and DBT were quantified in 18 out of 24 samples. The lowest TBT and DBT levels were found in autumn and spring (Table 1). TBT levels ranged from < LOD (0.8) to 118 ng Sn g −1 dry wt., with the highest concentration detected at sampling site S5 in the summer. DBT levels ranged from < LOD to 100 ng Sn g −1 dry wt. (winter, S2). MBT values ranged from < LOD to 51.35 ng Sn g −1 dry wt. (summer, S4). No seasonal differences were observed in BT levels for native mussels.
Regarding the condition index (CI), it ranged between 3.8 and 10.7 (Supplementary Fig. S1). No significant differences were found among the different sampling sites (Kruskal Wallis, p > 0.05), suggesting that the location in the estuary did not affect mussel characteristics. Moreover, the CI began to increase in the winter, reaching its maximum in the spring. However, no correlation between CI and BT concentrations was found. Rather, this behavior could be related to the availability of food during the so-called "phytoplankton bloom" (Guinder et al. 2015), which seems to be a key factor regulating CI (Arrieche et al. 2002).

Loading Plots
Sometimes, inspection of data tables is not enough to discover relationships between the different variables. In this sense, chemometric tools can help to interpret and clarify these hidden relationships. The experimental data inherently have a four-way structure, i.e., three BTs were analyzed in the four seasons in several compartments of different kinds for various sampling sites. Thus, Tucker4 was the selected multi-way method to extract useful information from the data set because of its capability to deal with different numbers of factors in each mode (i.e., sampling site, BTs, seasons, and environmental compartments). Supplementary Fig. S2 shows the deviance plot with the different model complexity (P × Q × R × S). Nine models were located on the lower boundary (the so-called convex hull). Ceulemans-Kiers st-coefficients were then calculated for these nine models, as shown in the table in Fig. S2 (an explanation of their calculation can be found in Stanimirova et al. 2009). The convex hull method makes it possible to identify suitable models in terms of low deviance and a high number of degrees of freedom. However, some points, representing the optimal models, lie on the line defined by the convex hull (some of them have a low deviance and number of degrees of freedom, while others have a high deviance and number of degrees of freedom), as shown in Fig. S2. To choose the optimal model among the candidates selected by the convex hull, the st criterion was used (Ceulemans and Kiers 2006;Stanimirova et al. 2009). The model with the highest stcoefficient was selected as the optimal choice. Therefore, the model with complexity (2 × 2 × 4 × 4) was chosen for the interpretation of BT distribution in BBE (see the table in the inset of Fig. S2). It is important to note that the low explained variance values (SS Fit /SS total ) obtained for the selected model are not uncommon when dealing with environmental data, which are affected by high noise due to weather conditions and to the relatively high experimental error in some variables (Leardi et al. 2000). Figure 2a shows the loading plot for the I-mode (sampling sites). All sampling sites show similar negative values for Factor 1 but present a greater variability for Factor 2. In this way, site S3 appears separately from the other sampling points. A possible explanation for this fact is that BT concentrations at this site were, in general, the highest recorded for most of the seasons in the different environmental compartments. Indeed, S3 is characterized by the presence of numerous fishing and leisure boats throughout the year, which could lead to high concentrations of TBT. On the other hand, the extreme sampling sites in the estuary (S1 and S6) appear together, while the other sites are grouped in the second quadrant of the Factor 1 versus Factor 2 plot. This behavior could indicate that these sampling sites have similar characteristics regarding BT concentration, although they are located at the extreme points of the estuary. In addition, sites S4 and S5 appear close to each other, which would indicate similar overall BT concentrations. Again, tidal currents seem to play an important role in the transport of pollutants through the BBE.
In the butyltin mode (J-mode), all BTs have negative values for Factor 1 (Fig. 2b) since this factor is negative and higher as the number of butyl substituents in the molecule is higher. Regarding Factor 2, TBT and MBT show negative values, whereas DBT has a positive one, which could be related to the different affinity of BTs for organic matter (Hongxia et al. 1996). In fact, the order according to Factor 2 is in agreement with the order of desorption reported by Dowson et al. (1996), i.e., DBT > TBT > MBT.
The K-mode (season) involves four factors. When Factors 1, 2, and 3 are plotted, each season appears in a different quadrant, which reflects the contrasts in the behavior of BTs at different times of the year (Fig. 2c shows, as an example, the plot of F2 vs. F1). On the contrary, all the seasons have negative values for Factor 4 (Fig. 2d only shows F4 vs. F3 for simplicity).
The L-mode also involves four factors and contains information about the environmental compartments (Figs. 2e, f show some combinations of factors to provide only substantial information and avoid repetition). In the graphs obtained for the compartment mode, the different sampling types have been grouped in a particular way. When Factor 1 is included in the graph, the mussels appear to be separated from the other compartments. Factor 2 then distinguishes between the water column and all other types of samples. Factor 3 groups water, mussels, and SPMs together, leaving sediment aside. Factor 4 separates sediment and SPMs from water and mussels.

The Core Elements
The distribution interpretation of the different variables in each mode is not always easy, and the relationship among these variables could be revealed by studying the elements of the core matrix. The core array obtained after Tucker4 decomposition (Eq. 1) contains information about the relationships among the different elements of the loading vectors in each mode. A better understanding of the system could be achieved by inspecting these relationships. The core elements arise from the product of the loading values, both in magnitude and sign (i.e., aip × bjq × ckr × dls). In this sense, many combinations of loading values are possible, taking into account that there are two factors for the first and second modes and four factors for the third and fourth modes. The most important core values are those that collect a high percentage of the variability observed in the core matrix (Stanimirova et al. 2009). Thus, Table 2 reports the four core values with the largest contributions.

The first Core Element
The first core element [1, 2, 1, 1] has a value of − 6.76, which accounts for the relationships between the first factors in modes I, K, and L and the second factor in the J-mode. Since the core value is negative, the loadings of each mode must be combined in such a way that the product among them yields a negative result. The loadings for Factor 1 in the I-mode are always negative (Fig. 2a), whereas Factor 2 in the J-mode, and Factor 1 in both K-and L-modes have both positive and negative values. Therefore, four combinations are possible to obtain a negative core value: The four possible relationships observed for the first core value are explained below: Fig. 2 a Loading plots for the I-mode (sampling sites), b loading plots for the J-mode (BTs), c-h loading plots for the K-mode (season) and i-n loading plots for the L-mode (compartments) obtained by the Tucker4 model

Combination 1A
As mentioned above, the loadings for Factor 1 in the I-mode are always negative and all sampling points contribute to this factor in a similar way. Therefore, no particular appreciation can be made of the sampling sites, which coincides with the Kruskal-Wallis test, for which no statistical differences were found among the sampling sites, as mentioned in Butyltin Levels and Physicochemical Variables Section.
The negative loadings for Factor 2 in the J-Mode correspond to MBT and TBT (Fig. 2b); the positive loadings in the K-mode are related to autumn (spring to a lesser extent); and the negative loadings for the L-mode are linked to sediment and water samples (SPM to a lesser extent). These relationships could be explained as follows: As a whole, the amounts of TBT and MBT were higher in autumn for both sediment and water samples. In this season, there is a notable increase in maritime traffic due to cereal export activities in the Bahia Blanca Port (from April to June), which may explain the increase in TBT levels. Furthermore, a recent study showed that recreational boats continue to be sources of BTs released into the environment due to their presence in historic paint layers on leisure vessels (Lagerström et al. 2017). The authors indicated that the old TBT-based coatings were not properly removed from the vessels, but only covered with a new coat of non-TBT paint. Therefore, this could be one of the causes of the presence of BTs in the BBE, in addition to the removal of industrial waste and sediment in the sampling area. As mentioned above, the tidal currents that characterize the study area affect the transport and recirculation of pollutants and could contribute to the dispersion of TBT in the estuary. These factors may have caused the BDI value (Table S1) calculated in autumn for the more internal sampling point (S6), which has a relatively low anthropogenic impact, to exceed the unit in both sediment and water samples.

Combination 1B
Another possible relationship established by the first core element links the concentration of DBT (positive contribution of Factor 2 to the J-mode) with autumn (more positive contribution of Factor 1 to the K-mode) and mussel samples (the only positive for Factor 1 in the L-mode). This could be due to the fact that the only season in which DBT was not detected in any sample of Brachidontes rodriguezii was autumn.
The levels below the LOD (for both TBT and DBT, see Table 1) could be related to the BTs elimination rate constant in mussels, which could be accelerated when the environment reaches high concentrations of these pollutants (Furdek et al. 2012;Tang et al. 2010). Some authors found in other bivalve species that the steady state was achieved earlier and that the elimination rate was higher when mollusks were subjected to larger amounts of TBT (Gomez-Ariza et al. 1999). This unusual behavior -low TBT levels in mussels versus high levels in the water column-could be explained by the accumulation potential of this particular mollusk species.

Combination 1C
The concentration of DBT (positive contribution to Factor 2 in the J-mode) is related to spring (the most negative contribution to Factor 1 in the K-mode) and to water and sediment samples (similar negative contribution of Factor 1 to the L-mode). This could be associated with the high DBT concentrations generally found in seawater and sediment samples in spring, which were higher than those determined for this compound in the other seasons. As stated in Butyltin Levels and Physicochemical Variables Section, in seawater, DBT was only determined in the spring. Furthermore, sediment samples collected this season at sites S2, S3, S4, and S5 showed a high concentration of DBT (Table 1).
Mounting evidence suggests that biotic degradation is one of the major pathways for the removal of TBT from the environment (Gadd 2000), turning it sequentially to DBT and MBT (Bernat and Długoński 2006). The increase in DBT concentrations both in the sediments and the water column in spring could be related to a degradation process caused by diatoms (Fang et al. 2017). In fact, the internal zone of the BBE is characterized by the flowering of diatoms during the winter-spring season (June-September) (Guinder et al. 2015), known as the "phytoplankton bloom." High DBT concentrations in spring in the water column could also be due to lower amounts of DBT adsorbed onto SPM (in other seasons, DBT in the SPM is much higher).
Furthermore, Dowson et al. (1996) have reported that DBT in sediments, once formed from TBT degradation, could undergo two different processes: on the one hand, it could be degraded to MBT and, on the other, it could be desorbed into the surrounding water column. The high rate of biological degradation in spring possibly generates a high amount of DBT, part of which is redissolved in the water column (Seligman et al. 1986(Seligman et al. , 1996. Likewise, Hongxia et al. (1996) have reported that DBT shows the lowest affinity for organic matter, so it tends to remain in the water column.

Combination 1D
Finally, MBT and TBT (the negative loadings for Factor 2 in the J-Mode) are also related to spring and mussel samples because these BTs show the lowest concentrations of these organisms in this season. An interesting finding is the relatively high correlation observed between TBT and DBT in spring (Spearman coefficient: r = 0.89, p < 0.05), which would indicate an active metabolism in mussels that results in the biotransformation of TBT into DBT (Furdek et al. 2012). Although these authors focused on a different mussel species, their physiology and reproduction strategies are very similar, as Mytilus galloprovincialis and Brachidontes rodrigezzi are both members of the Mytilidae family. Brachidontes rodriguezii (and Mytilidae species in general) enter the breeding season when the water temperature and food availability increase (Gosling 2003). These conditions are reached in spring and summer at the BBE, as well as during the big phytoplankton bloom (Domingues et al. 2005;Guinder et al. 2012). Furthermore, according to Adami et al. (2008), the Brachidontes rodriguezii population shows a high percentage of recruits throughout the year that peaks during the summer. This could also be related to the highest CI observed for mussels in spring, possibly due to an increase in nutrients during the phytoplankton bloom (Fig.  S1). Finally, the correlation between DBT and MBT was not significant, which may indicate that part of the DBT could be eliminated from the mussel (depuration) before being degraded to MBT or that the degradation metabolism of DBT is slower (Furdek et al. 2012).

The Second Core Element
The second element of the core array is [1, 1, 2, 2] with a value of − 6.15 (Table 2). As mentioned above, all the elements of the I-mode (sampling sites) for Factor 1 have negative values, a situation that is also observed for the elements of the J-mode (butyltins). The possible combinations in this case are the following: The explanation of these combinations is detailed below:

Combination 2A
This combination links the butyltins of the J-mode (particularly those with a higher, negative contribution, i.e., TBT and DBT) with the winter season (the most positive contribution of Factor 2 to the K-mode) and also with mussel samples (the most negative contribution of Factor 2 to the L-mode). It was observed that, in general, the mussel samples collected in winter presented the highest concentrations of the three BTs at all the sampling points. A rationale for this fact is that bivalves could have a limited capacity to metabolize BTs due to the low temperatures recorded in winter, thus generating a bioaccumulation effect (Furdek et al. 2012;Tang et al. 2010). Accordingly, the BDI values calculated for Brachidontes rodriguezii in winter were the highest in the study period, suggesting TBT accumulation. Previously, similar behavior was observed at the same sampling sites in 2013 (Quintas et al. 2017). Again, the higher amounts of TBT seem to be related to its availability for mussels due to the dredging carried out in the autumn of 2014 or to increased navigation in the Main Channel.

Combination 2B
The other relationship expressed by the second core element linked BTs with autumn and spring (negative loading values for Factor 2 in the K-mode) and with the water samples (the only one with positive loading for Factor 2 in the L-mode). Indeed, the overall concentration of BTs (as the sum of TBT, DBT, and MBT) in the water samples was higher at all sampling sites in the autumn and spring. In autumn, the high total concentration of BTs is mainly due to the contribution of TBT, which possibly appears as a consequence of the greater maritime traffic, as stated by other authors in different study areas (Champ et al. 2001;Hoch et al. 2001). In spring, the large amount of BTs in water is due to DBT, which was determined in high concentrations only in this season, as already discussed. On the other hand, the contribution of MBT to the total BT concentration is relatively high and shows little seasonal variation, probably due to the low degradation rate in comparison with the other BTs (Furdek et al. 2012;Wang et al. 2008).

The Third Core Element
The third element of the core is [1,1,3,4], with a positive value of 4.10. Again, two possible combinations could be found: A Since this core element has a smaller magnitude than the elements previously discussed, they involve the relationships between the elements with moderate to low contributions to the respective modes.

Combination 3A
This combination links the spring season (positive and lower contribution of Factor 3 to the K-mode) with SPM samples (positive and lower values for Factor 4 in the L-mode). It is noteworthy that all the SPM samples collected in the spring had MBT concentrations under the LOD, and the other BTs also showed their lowest concentration for this season. As mentioned before, the phytoplankton bloom takes place in winter/spring, and it has been well established that the growing biomass results not only in the adsorption of pollutants onto freshly produced particles, but also in higher biodegradation rates (Witt et al. 2002). This phenomenon has already been described in several coastal zones around the world, such as the Baltic Sea (Thorsson et al. 2008), the Swedish coasts (Gunnarsson et al. 2000), and the American Great Lakes (Taylor et al. 1991). Probably for this reason, the highest BDI values for SPM were recorded in spring in all sampling locations (Table S1), with some values close to the unit due to an accelerated degradation of DBT and MBT.

Combination 3B
The other possible combination in the third core element associates autumn (negative and low contributions of Factor 3 to the K-mode) and mussel samples (negative and low contributions of factor 4 to the L-mode). The relationship between mussels and autumn has already been commented on (first core element) and has to do with the low BT values recorded. In autumn, only MBT was quantified in mussels, with the only exception of TBT in site S3 (Table 1). In fact, the highest BDI value for mussels in all seasons was 2.12 (Table S1), recorded at sampling point S3 in autumn. This would indicate that the presence of TBT is common on this site, where numerous fishing and recreational vessels are observed on a regular basis.

The Fourth Core Element
The fourth core element is [2, 1, 1, 3] with a negative value of − 3.03 and a small magnitude compared with the first and second core elements, so the lowest contributions to the loadings are considered. The possible loading combinations are: The explanation of these relationships is detailed below:

Combination 4A
In this core element, the sampling points (I-mode) with positive values are sites S4 and S5, which are related to MBT (lower and negative for Factor 1 in the J-mode), the summer season (low and positive element for Factor 1 in the K-mode), mussels, and water samples (lower and positive values for Factor 3 in L-mode). In fact, sampling point S4 had the highest MBT value for mussels throughout the year. For water column samples, the S4 point also recorded the highest value for MBT in the summer. The BT degradation process seems to be favored at this sampling site, which is also reflected in the low BDI values (Table S1). In turn, sampling site S5 showed the highest value of TBT throughout the year for mussels. Indeed, S5 in the summer was one of the two points with a BDI value greater than unity. As mentioned above, in the high summer water temperatures, a higher metabolic TBT rate is to be expected. In fact, BDI values for all sampling sites were between 0.02 and 0.68, except for site S5 where the BDI value exceeded unity (1.65). This site, in particular, is located in the estuary's interior, with a small leisure and fishing vessel traffic, primarily during the summer. Therefore, the mussels from this site have probably been more exposed to TBT.

Combination 4B
In this combination, the same sampling points (S4 and S5) appear related to the winter season (the lower and negative values for Factor 1 in the K-mode) and sediment samples (the only negative element for factor 3 in the L-mode). These sites displayed the highest concentration of MBT for the season. This fact, coupled with the low BDI value, would indicate an efficient BT degradation process.

Combination 4C
Another important relationship found in the analysis of the fourth core element is the link of sampling point S3 (negative value for Factor 2 in the I-mode) with MBT concentrations (the lower and negative contribution of Factor 1 to the J-mode), the summer season (positive and lower contribution of Factor 1 to the K-mode), and sediment samples (the negative element of Factor 3 in the L-mode). This is to be expected since this sampling point showed the highest MBT concentration in summer for sediment samples, which could be related to active biological degradation of BTs at this location.

Combination 4D
The last relationship linked sampling point S3 (again, the negative element of Factor 2 in the I-mode), MBT (lower and negative contribution of Factor 1 to the J-mode), the winter season (lower and negative loading for Factor 1 in the K-mode), mussels, and water samples (positive and lower contribution of Factor 3 to the L-mode). The resulting relationship is not clear enough but could be associated with the relatively large BDI value, mainly in water samples, evidenced in the low concentration of MBT (and high concentration of TBT) in this sampling site in winter.

Partitioning of TBT among Different Compartments
Sorption is considered one of the most important processes responsible for the reduction in the concentrations and toxicity of organotin compounds in the marine environment (Fent 1996). Therefore, apparent distribution coefficients for sediment-water k s w , mussel-water k m w , mussel-sediment were calculated as the ratio of TBT concentrations in the different matrices to assess their distribution in seawater, sediments, mussels, and SPM (Table 3). The coefficients k m w ,k m s , and k m SPM can be used to evaluate the ability of mussels to concentrate TBT from water, sediment, and SPM, respectively. In the current study, k m w ranged from 0 to 2.09 (mean: 0.47), k m s from 0 to 1.03 (mean: 0.24), and k m SPM from 0 to 4.47 (mean: 0.89). In the winter, the three coefficients had the highest values. As mentioned above, at low temperatures, the ability to metabolize TBT is reduced, so mussels tend to accumulate this compound (Furdek et al. 2012;Tang et al. 2010). By contrast, the lowest values for these coefficients were found in autumn (Table 3). This is so because the pumping activity of mussels (Brachidontes rodriguezzi) is clearly reduced at high TBT burdens, which also may lead to a decrease in TBT uptake efficiency (Choi et al. 2009;Quintas et al. 2017). In autumn, the levels of TBT in water ranged from 168 to 397 ng Sn L −1 , while its presence could not be detected in mussel samples. The TBT concentrations present in the water column exceed the upper EAC value (0.62 ng Sn L −1 ) (OSPAR 2004). Therefore, certain biological effects are probable for some species located in the area in contact with this water. According to the average values obtained, Brachidontes rodriguezzi seems to bioaccumulate TBT from sediment, water, and, to a lesser extent, SPM. The distribution of TBT between the dissolved and particulate phases is the main factor that determines its bioavailability for degradation by microorganisms, as well as its deposition into sediments (Gao et al. 2019). In this respect, k s w values ranged from 0.12 to 3.46 (mean: 0.87), while k s SPM ranged from 0 to 3.26 (mean: 0.63). These values were very similar, indicating that TBT could be transferred to the sediments from both water and SPM. In turn, the k w SPM coefficient displayed a wide range (from 0.06 to 12.32, average: 1.98), indicating that SPM has a high sorption capacity and does not reach equilibrium with seawater (Choi et al. 2009). The k w SPM values suggest that almost all the dissolved TBT is adsorbed onto suspended particles and could then accumulate in sediments (Ritsema et al. 1991).

Conclusions
This study evaluated for the first time the distribution and accumulation of BTs in the Bahía Blanca Estuary (BBE) through a comprehensive analysis performed in both temporal and spatial domains. Chemometric tools, such as the Tucker4 model, made it possible to deeply analyze, visualize, and interpret the information hidden in the data tables and to discover the relationships among BTs, sampling sites, seasons, and the environmental matrices, including the water column, sediments, SPM, and native mussels. The results pointed out that all compartments were impacted by the presence of the three BTs, which showed a rather homogeneous distribution throughout the estuary, even though the use of TBT has been banned since September 2008. The presence of these compounds in BBE can have different origins; that is, industrial waste from the PVC plant, the resuspension of sediments, or the leaching of TBT from antifouling paints not properly removed before the use of TBT-free paints. Therefore, it is necessary to further monitor and assess BTs in BBE and to continue with the TBT ban to protect the aquatic environment. This fact highlights the importance of pollutant transport by tidal dynamics and other anthropic factors, such as periodic dredging that could remobilize the contaminants trapped in sediments. Moreover, other factors, like the phytoplankton bloom, undoubtedly affect BT degradation and distribution in the different environmental compartments. Furthermore, the apparent distribution coefficients also reveal that TBT bioaccumulates from sediment, water and, to a lesser extent, SPM, which seems to play a critical role as a vector for BT deposition in sediments.