Environ. Sci. Technol. 1997, 31, 1411-1415
U.S. Mussel Watch Data from 1986 to 1994: Temporal Trend Detection at Large Spatial Scales B . B E L I A E F F , †,| T H O M A S P . O ’ C O N N O R , * ,‡ D. K. DASKALAKIS,‡ AND P. J. SMITH§ IFREMER, B.P. 1105, 44311 Nantes Cedex 03, France, NOAA NORCA21, 1305 East West Highway, Silver Spring, Maryland 20910, and Department of Mathematics, University of Maryland, College Park, Maryland 20742
National Mussel Watch Programs collect bivalves at many sites along their respective coastlines in order to assess chemical contamination of coastal waters. This paper aims at detecting temporal trends at large spatial scales using contaminant data from 1986 to 1994 of the U.S. National Status and Trends Mussel Watch Program. U.S. coasts are divided into 10 large areas. After a logarithmic transformation of the data, linear models are used and include a ‘sampling site’ effect, a linear trend estimated for each area, weather covariates, and an ‘analysis’ effect reflecting changes of analytical procedures or even laboratories. Except for the Californian coast, chlorinated organic compounds show rather uniform decreases, which can be attributed to bans put on these contaminants. A few increasing trends were detected for trace elements, which also showed few decreases. Additions of weather covariates and possible analysis effects were intended to model apparently random between-year variations in contaminant time-series. Although low in terms of percentage of variability explained, the addition of these effects might have avoided the detection of spurious trends for some contaminants.
Introduction National Mussel Watch Programs collect mussels and oysters at many sites in order to cover coastlines in a roughly systematic way. For example, the 250 sites used by the National Status & Trends (NS&T) Program allow a sensible spatial grid for the U.S. coastlines (1). For the same objectives, samples are collected at about 100 sites for the French Monitoring Network (2). Monotonic trends were detected in annually measured concentrations at numerous NS&T sites for the same contaminant (1), e.g., ∑chlordane, but the approach used does not allow a spatial analysis of trends. Plots of NS&T data have shown time-series behavior for sites located in the same geographical area or for sites that are not spatially related. These similarities might have two origins. It is well known that contaminant concentration in biota exhibits seasonal variations. Mollusks are collected at the * Corresponding author fax: 301-713-4388; e-mail address:
[email protected]. † IFREMER. | Present address: ARA, Environment Department, P.O. Box 21, Aqaba, Jordan. ‡ NOAA NORCA21. § University of Maryland.
S0013-936X(96)00658-X CCC: $14.00
1997 American Chemical Society
same time each year, but the physiological cycles of the bivalves are subject to shifts within 1 year due to changes in environmental conditions and in particular temperature (3). Another possible explanation would correspond to biases in analytical procedures. This hypothesis would favor similar temporal patterns for sites located on different coasts and analyzed by the same laboratory. The aim of this paper is to look for contaminant trends at large spatial scales using the U.S. Mussel Watch data from 1986 to 1994. In this study, an analysis effect and weather covariates are introduced in linear models.
Materials and Methods NS&T Data. A fully detailed description of the U.S. NS&T Mussel Watch Project sampling strategy is given in ref 4. The following contaminant data have been used: As, Cd, Cu, Hg, Ni, Se, and Zn for the trace elements and total chlordane (∑chlordane), total dieldrin (∑dieldrin), total DDT (∑DDT), PCBs (∑PCB), and PAHs (∑PAH) for the organic compounds. As defined in ref 1, ∑PCB is twice the sum of concentrations of 18 PCB congeners; ∑PAH is the sum of concentrations of 24 compounds; ∑DDT is the sum of the concentrations of the 2,4′- and 4,4′-forms of DDT, DDD, and DDE; ∑chlordane is the sum of concentrations of R-chlordane, trans-nonachlor, heptachlor, and heptachlorepoxide; and ∑dieldrin is the sum of concentrations of dieldrin and aldrin. We selected 150 NS&T sites with at least six observations (Figure 1). We did not use data from the four sites in Hawaii and Alaska. Ten geographical areas were defined (Table 1) on the basis of medians and ranges of air temperature with three additional constraints. First, the areas correspond to sets of geographically contiguous sites. Second, in each year, samples from an area are analyzed by the same laboratory. Third, all samples in an area are either mussels or oysters. There are sites in Long Island Sound (within area E2) where both mussels and oysters have been collected, but only mussels have been collected consistently and only mussels have been used for analysis of trends. There is a West Coast site where Mytilus edulis is routinely collected but where, in 1 year, two species of mussels were collected. Both M. edulis and Mytilus californianus had very similar concentrations for all chemicals (1), and this allows us to include both M. edulis and M. californianus sites in areas W1 and W2. Information about laboratories are given in ref 5: the East Coast samples (areas E1-E4) have been analyzed by Battelle Ocean Sciences (Duxbury, MA); the Gulf samples (areas G1G3) have been analyzed by the Geochemical and Environmental Research Group of Texas A&M University (College Station, TX); the California samples (areas W1 and W2) have been analyzed by Science Applications International Corporation Inc. (La Jolla, CA) until 1989 and by Battelle Ocean Sciences (Sequim, WA) from 1990; the Oregon and Washington samples (area W3) have always been analyzed by Battelle. Weather Data. Weather data have been extracted from the Summary of the day Cd-ROMs, provided by the National Climatic Data Center (Asheville, NC). Eighty-seven weather stations have been selected on the basis of their distance from NS&T sites, their records including the period 19861994, and their records having less than 10% missing values. The Model. A logarithmic (base 10) transformation of the concentration data was performed, in particular to ensure an approximate normal distribution for the residuals. Therefore, in the following, concentrations are log10 concentrations. From 1986 through 1991, three separate composites were collected
VOL. 31, NO. 5, 1997 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
1411
FIGURE 1. Map of NS&T site locations and areas.
TABLE 1. Area Definition coast
area
no. of sites
East
E1 E2 E3 E4 G1 G2 G3 W1 W2 W3
9 21 11 10 8 28 19 12 17 15
Gulf
West
TABLE 2. Occurence of Analytical Changes species
Tmedian
Trange
1.7 6.4 8.3 15 18.9 15 17.2 13.3 8.3 6.7
0-5 3.3-10.6 6.4-10.6 11.7-17.8 17.2-20 12.2-16.7 11.7-19.4 8.3-15 5-11.7 3.9-8.3
Mytilus edulis Crassostrea virginica Crassostrea virginica M. edulis and M. californianus
at each site in each year. In this case, an arithmetic mean was computed on the log-transformed values. This is strictly equivalent to a geometric mean computed on the raw data. A multiple linear regression is used for trend estimation. A detailed description of this statistical tool can be found in ref 6. Briefly, the concentration (Y) measured at the jth site in the ith area for the kth year is described by the model l)1
Yijk ) µ + δi + τj(i) + γitk + ξik +
∑R X
l lk
+ ijk
(1)
q
where µ is a constant corresponding to the unknown general mean, δi corresponds to the effect of the ith area on the general mean, τj(i) corresponds to the effect of the jth site nested in the ith area, γi corresponds to the slope of the interaction cross-product between area j and the continuous variable year t, and Rl is the coefficient of a weather covariate Xl (e.g., air temperature or precipitation), with q as the number of covariates. ξik corresponds to a ‘dummy’ variable (ref 6, pp 134-142) allowing us to model a possible vertical shift in γi due to an analysis effect (either a change of laboratory or a change in analytical methods, or both). In this model, the total number of parameters to be estimated ranged from 164 to 166, the difference being due to the number of ξik coefficients for a given contaminant. It should be clear that the area effect, δi, includes concentration differences among areas due to differences in temperature and species.
1412
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 31, NO. 5, 1997
contaminant
year
areas under concern
all all trace elements Ni, Pb, Cd Cd Cu all organics ∑PAH ∑PAH ∑PAH
1990 1992 1992 1987 1988 1989 1987 1988 1989
W1, W2 Gulf Coast East and West Coasts W1, W2 W1, W2 Gulf Coast East Coast and W3 East Coast and W3 W1, W2
The residuals ijk are assumed to be independent and normally distributed, with mean zero and constant variance. Normality of residuals was checked using normal probability plots. The model parameters are simultaneously estimated using standard regression techniques to minimize the squared sum of differences between observations and estimations. Then their significance is statistically assessed: the basic principle is to check whether the decrease of the residual sum of squares induced by the addition of a new parameter in the model is significant or not. For details on this extra sum of squares, see ref 6 (pp 67-69). Table 2 gives the years for which there were changes either in laboratory and/or analytical procedures. Details in these changes can be found in refs 1 and 5. For example, a change of laboratory, already mentionned in NS&T Data, occurred for the California Coast in 1990. Adding the ξi step function does not mean that it is certain that contaminant results are biased for these samples, but it allows us to take into account this possible analysis effect. The weather covariates used are the daily minimum temperature averaged over the 2 months preceeding the sample collection month (T) and the precipitations cumulated over the 2 months preceeding the sample collection month (P).
Results and Discussion Detailed Example of Model Results. To provide an example on how the model works, we selected selenium data for the
TABLE 3. Percentage of Total Variability Explained by Different Effectsa As Cd Cu Hg Ni Pb Se Zn ∑chlordane ∑dieldrin ∑DDT ∑PCB ∑PAH
p
N
site
trend
W+L
R2
164 167 165 164 167 167 164 164 164 164 164 164 167
1146 1145 1146 1146 1145 1143 1143 1146 1111 1053 1131 1126 1065
82.2 76.6 96.3 71.5 68.6 88.6 55.1 96.8 64.5 61.6 81.9 81.8 73.0
1.6 2.7 0.3 0.7 1.2 0.9 2.5 0.1 11.9 6.7 2.4 4.0 0.9
0.8 0.8 0.1 -0.2 0.7 0.2 1.1 0.1 0.5 0.4 0.1 0.4 2.4
84.5 80.1 96.7 72.0 70.3 89.7 60.7 97.0 76.9 68.9 84.4 86.2 76.3
a W + L ) weather effects [T and P] + analysis effect; R 2 ) final R 2 value (%). p is the number of estimated parameters for the final model, and N is the number of observations.
TABLE 4. Significant (p < 0.005) Trend Resultsa E1
FIGURE 2. Model result example for site j ) NBWJ (area i ) W1). 1 2 3 Yijk are the intial log-concentration values. Y ˆ ijk ,Y ˆ ijk , and Y ˆ ijk are the residuals after having added respectively and successively the site effect, the weather covariates, and the temporal trend into the 4 ˆ ijk is the residuals after having modeled the analysis model. Y 4 effect; finally Y ˆ ijk ) Eˆ ijk.
As Cd Cu Hg Ni Pb Se Zn ∑chlordane ∑dieldrin ∑DDT ∑PCB ∑PAH a
NS&T site NBWJ in area W1. Selenium contamination at year k can be estimated as
Yˆ W1NBWJk ) 0.9231 + 0.2050 - 0.0952 - 0.0028TNBWJk 0.0027RNBWJk - 0.0055tk + ξˆ W1k (2) where 0.9231 is the national grand mean for selenium (in logarithmic scale); +0.2050 is the adjustment to the national mean required for area W1; -0.0952 is the adjustment required for site NBWJ in area W1; -0.0055 is the slope corresponding to area W1; -0.0028Tk and -0.0027Pk are the effects of time (year), temperature (°F), and precipitation (in.), respectively, on site NBWJ mean; and the effect of laboratory ξˆ W1k is 0.0853 in 1986-1989 and 0 thereafter, i.e., if k > 89. An estimation of the residual error W1NBWJk is
ˆW1NBWJk ) YW1NBWJk - Yˆ W1NBWJk
(3)
For example, 0.4015 log10 [µg (g dry weight)-1] is the selenium concentration observed in January 1986 (tk ) 86), where TNBWJ86 and PNBWJ86 were respectively 70.5 °F (21.4 °C) and 4.15 in. (0.11 m) (respectively mean and sum of November and December 1985 monthly data). Finally, the residual value, i.e., the part of the concentration value not explained by the model, is W1NBWJ86 ) 0.4015 - 0.4343 ) -0.0328 log10 [µg (g dry weight)-1]. Figure 2 shows model results with sequential deletions of possible effects. It allows one to consider residuals obtained adding effects sequentially in the model. As an illustration, the impact of adding an analysis effect can be seen when 3 4 moving from Yˆ ijk to Yˆ ijk , which is the final residual value. In this particular case, the trend estimate goes from a highly significant decrease (p < 10-4) to non-significance once the
E2
E3
E4
I D D
G1
G2
G3
W1
W2
W3
D D D
D D D D
D D
D D D D
D D D
D D D
D
I I D
D D
D D
D D D D
D D D
D
D, significant linear decrease; I, significant linear increase.
step function ξˆ W1k is added. Note that ˆW1NBWJ86 is equivalent 4 to the 1986 value of Yˆ ijk . Table 3 gives the percentage of total variability explained by the different effects added to the model. Most of the variability is explained by the site effect, i.e., δi + τj(i). The relative contribution of the temporal trend, once the model is adjusted for the site effect, is the ratio of the percentage given by Table 3 for the trend to the remaining variability. For example, this contribution is 9% for arsenic [100 × 1.6/ (100-82.2)]. The highest contributions for trace elements are Cd (12%) and As (9%) and for organics are ∑chlordane (34%), ∑PCB (22%), ∑dieldrin (17%), and ∑DDT (13%). The percentage of the combination of weather and analysis effects is very low. The highest contributions are obtained for selenium, arsenic, cadmium, and nickel for the trace elements and are obtained for ∑PAH for the organics. The loss of observations due to missing values for weather covariates induced a reduction in R 2. Thus, the mercury case shows a loss (-0.2) in the percentage of variability due to negligible influences of weather covariates and of the analysis effect. Trends. Trend results per area are presented in Table 4. We only considered trends that are significant (p < 0.005) using the final model (1) and that remained significant in the case where the analysis effect is discarded from the model. On the contrary, ‘new’ trendsstrends appearing once the analysis effect is finally added to the modelsare not included. These trends might be true, but given the model’s sensitivity to the addition of new effects due to the scarcity of data for each site, it is safer to consider them as possible numerical artifacts. Among the trends we accept as real, three are increases: arsenic in E2 and selenium and zinc in G2. For the trace
VOL. 31, NO. 5, 1997 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
1413
TABLE 5. Significant Parameters for Analysis Effects with Areas under Concern and for Weather Effects (rˆ T and rˆ P) ξˆ ik As Cd Hg Se Zn ∑chlordane ∑dieldrin ∑PCB
areas
-0.092 for k < 1990 W1, W2 W1, W2 0.155 for k < 1990 W1, W2 0.085 for k < 1990 -0.089 for k < 1990 W1, W2
rˆ T
-0.0032 -0.0026 -0.0029 -0.0028 -0.0027 -0.0034
0.260 for k < 1990 0.145 for k > 1988
rˆ P
0.0041
W1, W2 Gulf Coast
elements, the highest number of decreasing trends is obtained for copper (E2, G3, and W3) and cadmium (E2, G3, and W2). A previous study has shown a decrease of cadmium on a national scale (1), which does not appear here. In that study, two approaches were used to assess trends on a national scale. First, the number of decreases detected at individual sites was tested against a randomness hypothesis considering the total number of sites. This shows a significant number of decreases, but it should be noted that a high number of decreases were obtained for neighboring sites. Second, using a Spearman test, it was shown that a significant correlation existed between year and annual geometric mean based on data for all sites. However, we should be cautious on the national interpretation of this trend. It is not necessary for cadmium concentrations to change at all sites to alter the national geometric mean. Except for ∑PAH, there is a strong tendency for organics to decrease as already stated and interpretated in ref 1. Although still a rather crude spatial analysis, dividing the U.S. coasts into areas reveals some regional differences in trends. For example, we did not detect any significant decreases for the W1 and W2 areas, constituting the California Coast, for ∑chlordane, ∑dieldrin, and ∑DDT. This appears consistent with data in ref 7, showing sharp decreases in chlordane, DDT, or PCB in the early 1980s, but with generally more stable levels from 1985 to 1986, corresponding to the start of the NS&T Mussel Watch Program. ∑chlordane and zinc display two extremes in the model resuls. The site effect explains 96.8% of the total variability for zinc, showing low between-year variabilities, but it accounts for only 64.5% for ∑chlordane (Table 3). The remaining 35.5% corresponds to between-year variability, i.e., either random or systematic changes. Among these possible changes, modeling a linear trend per area extracted another 11.9% from the residual variability, while doing so for zinc added a negligible 0.1% for zinc. Zinc or copper are essential elements whose concentrations are to some extent regulated by mollusks (8). This might explain the lower between-year changes and the ability to detect small trends. Thus trends were detected for copper in areas (E2, G3, and W3) with slopes 5 times lower than the slopes obtained for significantly decreasing ∑chlordane. Random Between-Year Changes. Significant contributions for weather and analysis effects are described in Table 5. In all but one of those cases, the added analysis effect is the change of laboratory for samples from the California Coast, from 1990. We found a few significant contributions of weather covariates and always a negative effect, except for ∑chlordane, for precipitation and temperature (Table 5), e.g., for cadmium, as found in ref 3. The reader must be aware that most of the weather influence is already captured by the site effect. Higher mean levels for some sets of sites can, of course, be due to differences in contamination but could also be due to differences in bioaccumulation between mussels and oysters. What weather covariates intend to model are among-year variations due to the possible influence of environmental factors. Precipitation can affect the transfer of contaminants from land to the coastal zone, and it can
1414
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 31, NO. 5, 1997
affect salinity which itself influences the availability of contaminants. Temperature affects bivalve spawning period, which in turn affects the masses of tissue, lipid, and chemicals within mollusks (8). Although significant, the contribution of weather covariates to the explanation of the residual variability is very low and here does not exceed 1%. An obvious cause for this is that temperature and precipitation, by themselves, cannot quantitatively account for complex indirect relationships among environmental, physiological, and biochemical processes involved in contaminant bioaccumulation. Possible model refinements could be brought by the addition of covariates, if available, more directly related to physiological processes in bivalves. The addition of weather covariates and analysis step functions was intended to model between-year changes that might spuriously appear as trends. The detailed example given in Results stresses the sensitivity of trend results to the addition of this set of continuous (weather) or categorical (analysis effect) variables. In the example, a possible spurious trend disappeared when modeling a change of laboratory for California sites. In other cases, new trends appeared after modeling the analytical changes or the influence of meteorological factors. In this case, the noise reduction allows the trend to be detected. These new trends were not considered as true ones because we cannot distinguish beween a true trend and a spurious trend due to a signal distortion, induced by the rather crude way of modeling these changes. With the NOAA Mussel Watch data for 1986-1989, Wilson et al. (9) sought concordant shifts over large regions of the Gulf of Mexico in concentrations in oysters of trace elements and PAHs. They found more than random concordance among selenium concentrations throughout the Gulf. Arsenic concentrations also varied in close to unison except along the Louisiana Coast. Mercury and cadmium behaved concordantly among sites in the Western Gulf. Copper, zinc, and PAH concentrations varied together in the eastern and western Gulf. These large-scale concordances were attributed to similarly large-scale influences that the authors speculated could include El Nino/southern oscillation events. With so few years of data, a monotonic trend could not have been estimated by these authors. However, if an actual trend exists in the data, it is likely to modify concordance results. Concordance among annual shifts in concentration do not necessarily correspond to linear trends, but linear trends might be seen as successive concordance shifts. Our results (Table 4) based on 9 years of data include some trends for As, Cd, Cu, Se, and Zn in large regions of the Gulf. The trends for As, Cd, Cu, and Se were in regions where Wilson et al. (9) found concordant behavior. Linear models were already used to explain variability in chemical contamination in bivalves, but with quarterly data and with a spatial scale restricted to an estuary (3). The data discussed here are based on a different sampling strategy. But again, a linear model appears to be a flexible and powerful statistical tool, allowing inclusion of continuous covariates and of ‘dummy’ variables like those used to model analysis effects.
Literature Cited (1) (2) (3) (4)
O’Connor, T. P. Mar. Environ. Res. 1996, 41, 183-200. Claisse, D. Mar. Pollut. Bull. 1989, 20, 595-603. Beliaeff, B.; Smith, P. J. Water Resour. Bull. 1996, 32, 1-9. O’Connor, T. P. In Environmental Statistics, Assessment, and Forecasting; Cothern, C. R., Ross, N. P., Eds.; Lewis Publishers: Boca Raton, 1994; pp 331-349. (5) Lauenstein, G. G.; Cantillo A. Y. Sampling and analytical methods of the NS&T Program National Benthic Surveillance and Mussel Watch Projects. Volume I. Overview and summary of methods; NOAA Technical Memorandum NOS ORCA 71; NOAA: Silver Spring, MD, 1993. (6) Draper, N. R.; Smith H. Applied Regression Analysis; Wiley: New York, 1981.
(7) Stephenson, M. D.; Martin, M.; Tjeerdema R. S. Arch. Environ. Contam. Toxicol. 1995, 28, 443-450. (8) Phillips, D. J. H.; Rainbow, P. S. Biomonitoring of Trace Aquatic Contaminants; Elsevier Applied Science: New York, 1993. (9) Wilson, E. A.; Powell, E. N.; Wade, T. L.; Taylor, R. J.; Presley, B. J.; Brooks, J. M. P. S. Helgol. Meersunters. 1992, 46, 201-235.
Received for review July 30, 1996. Revised manuscript received January 17, 1997. Accepted January 23, 1997.X ES9606586 X
Abstract published in Advance ACS Abstracts, April 1, 1997.
VOL. 31, NO. 5, 1997 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
1415
