Effects of uranium mining and milling on benthic invertebrate communities in the Athabasca Basin of Northern Saskatchewan

Abstract

Aquatic environments downstream of uranium operations (mining/milling) in Northern Saskatchewan are exposed to a variety of chemical and physical disturbances. There is extensive regulatory documentation of monitoring and effects assessments for modern uranium operations, though little is available in the scientific literature. The reference condition approach was used to predict expected benthic invertebrate community metrics for unexposed and exposed lakes. This approach was used to identify impacted communities that differed significantly from expected natural regional compositions. The number of taxa at the lowest practical level (taxon richness) downstream of approximately half of the uranium operations’ effluent release points was lower compared to the number of taxa observed in reference conditions. The taxon richness residuals were strongly correlated with about half of the measured concentrations of sediment contaminants and always negatively, such that richness was lower in lakes with higher concentrations of metals or radionuclide activities in the sediments. This was consistent for both reference and exposure lakes, implying a natural background influence of some metals and radionuclides on taxon richness. This exercise produced models that can be used by practitioners in the north of Saskatchewan in the design and delivery of monitoring programs not only for uranium mines, but for regional cumulative effects studies in support of sustainable resource planning and management.

Les milieux aquatiques en aval d’exploitations uranifères (minière/extraction) dans le nord de la Saskatchewan sont exposés à diverses perturbations chimiques et physiques. Il existe une documentation réglementaire exhaustive sur la surveillance et les évaluations des effets des opérations modernes uranifères, bien qu’il y ait peu d’information dans la littérature scientifique. L’approche fondée sur les conditions de référence (ACR) a été utilisée pour prédire ce qui est attendu des communautés d’invertébrés benthiques dans des lacs exposés et non exposés. Le nombre de taxons identifié (richesse taxonomique) en aval de la moitié des points de rejet des effluents des opérations minières et d’extraction d’uranium était inférieur au nombre de taxons observé dans les conditions de référence. Il y avait une forte corrélation négative entre les résidus de la richesse taxonomique et la concentration de la moitié des contaminants mesurés dans les sédiments, de sorte que la richesse était plus faible dans les lacs ayant des concentrations élevées de métaux ou de radionucléides dans les sédiments. Cette relation a été observée dans les lacs de références ainsi que les lacs exposés, ce qui implique une influence naturelle de certains métaux et radionucléides sur la richesse taxonomique. Cet exercice a produit des modèles qui peuvent être utilisés pour la conception de programmes de surveillance dans le nord de la Saskatchewan, non seulement pour les mines et usines d’extraction d’uranium mais aussi pour des études régionales des effets cumulatifs servant à la planification et à la gestion de l’exploitation durable des ressources.

Introduction

Aquatic environments downstream of uranium mining and milling operations (hereafter mining) are exposed to a variety of both chemical and physical disturbances. The individual and combined effects of released contaminants on freshwater organisms have been tabulated in the gray literature, but have rarely been analyzed in the scientific literature (Peterson et al. 2002; Mihok and Thompson 2011) with the exception of selenium (Muscatello et al. 2008; Muscatello and Liber 2009; Wiramanden et al. 2010; Gallego-Gallegos et al. 2012; Tse et al. 2012; Franz et al. 2013). There is a wealth of underutilized data on benthic invertebrate community composition in this context from metal mines in Canada. There have, however, been no substantive peer-reviewed examinations of mine-related data exploring the influences of uranium mining on lake or stream benthic communities.

All active Canadian uranium mines and mills are in the Athabasca Basin in Northern Saskatchewan. The geology of the Athabasca Basin results in natural concentrations of some metals and radionuclides in surface water and sediments that are higher than in other regions. Typical aquatic environment monitoring programs include the sampling of water, sediment and fish tissue chemistry, sentinel fish health and benthic invertebrate community composition. Benthic invertebrate communities are monitored as they are an important part of the aquatic food chain and possess many of the characteristics of an ‘ideal’ sentinel organism. These include being sedentary in nature and sufficiently long lived to allow simultaneous sampling of multiple age classes and developmental stages, thereby reflecting site-specific exposure to both water-borne and sediment-accumulated contaminants (Mousavi et al. 2003; Bonada et al. 2006).

The monitoring of lakes is a key aspect of the regulatory compliance monitoring of active uranium mines and mills. There is, however, a challenge to monitoring benthic invertebrates in lakes as it is often difficult to find suitable reference lakes (unexposed to mining operations) that are good matches for lakes exposed to mining influences (Bonada et al. 2006; Mackey et al. 2013). The composition of benthic invertebrate communities among lakes varies predictably with average lake depth, depth of sample collection and substrate type, among a variety of other subtler influences (Brinkhurst 1974; Saether 1979). Results of monitoring programs that describe differences in composition of fauna between reference and exposed lakes are often confounded with natural differences in physical, chemical and biological characteristics.

Extensive environmental monitoring is undertaken in the Athabasca Basin of Northern Saskatchewan due to the association of uranium deposits with this geological formation. Cameco Corporation (Cameco) and Areva Resources Canada Inc. (Areva) undertake the monitoring at the five operating uranium mines and/or mills concentrated along the eastern margin of the basin (Figure 1). Collectively they have developed a large database of benthic invertebrate communities collected from dozens of reference and exposed lakes using the same methods. These data provide an opportunity to develop a regional reference condition approach (RCA) model (Bailey et al. 2004; Bailey et al. 2014) for benthic communities from reference lakes. RCA models use ‘template’ variables as defined by Imhof et al. (1996) to predict an expected reference condition in the absence of human influences. RCA models are most typically applied to assessments of benthic invertebrate communities, but are also developed and applied to surveys of fish communities (Stanfield and Kilgour 2006). RCA-style models for lakes are not common, and have been generally restricted to models predicting lake chemistry or algal community composition (Hawkins et al. 2010). That said, Reynoldson et al. (1995) used RCA models to assess the benthic invertebrates of areas of concern in the Great Lakes, and Parsons et al. (2010) developed an RCA model for benthic invertebrates in small, shallow (generally <1 m deep) lakes in the vicinity of oil sands operators near Fort McMurray, Alberta. There is every expectation that RCA models for lake benthic communities could be used to predict expected benthic communities for lakes exposed to mine and mill operations and for which there are no nearby and highly similar reference lakes, and will thus improve assessments of benthic invertebrates from exposed lakes.

Figure 1. Map of Northern Saskatchewan showing locations of reference lakes, and lakes assessed using a reference condition approach model. Samples of benthos and sediment were collected from these lakes at various times from 2002 to 2009, from various depths, and used in the models developed herein.

Cameco and Areva monitor benthic invertebrate communities in lakes exposed to treated effluent to meet their regulatory requirements under the Nuclear Safety and Control Act and the Metal Mining Effluent Regulations (MMER). Each of the five active uranium mines and/or mills in Northern Saskatchewan has its own independent program predominantly consisting of one or two reference areas and multiple downstream exposure areas. Some mining/milling operations have adopted the use of multiple reference areas in more recent years. Combining data from all the reference areas collected from the Athabasca Basin presents an opportunity to explore causes of variation in benthic community composition. In this study, conventional and multivariate indices were used to summarize variation in benthic invertebrate community compositions among lakes considered to be in a reference condition on the basis that they are not exposed to influences of uranium mining operations, and because other potential influences are minimal. There is no logging in the region, so forests are generally intact except for natural disturbance due to fire. There are minimal road networks, and no use of salts on roads for winter de-icing. There is no agriculture. There is some angling and netting of fish by local community members but because the area is remote, the impact is expected to be minimal.

Reference condition models were constructed using data from reference monitoring areas. The reference models were used to quantify the degree of deviation from reference for benthic communities in lakes influenced by uranium mining operations. Deviations from the reference condition, for specific indices of benthic community composition, were examined for correlation with sediment chemistry variables, specifically concentrations of inorganic constituents (metals, metalloids), to identify potential causes of variation in exposed benthos community indices. The models developed here can be used by industry to enhance existing monitoring activities by providing a modeled prediction for lakes for which an exact local reference lake may not be available, and through the comparison of local benthic community variations to regional variations that provide a broader ecological perspective (Kilgour et al. 1998, 2017).

Methods

The lakes from which benthic invertebrates were collected and incorporated into this database/analysis are illustrated in Figure 1. The data were collected by Cameco and Areva in support of their Environmental Effects Monitoring (EEM) reporting requirements under the MMER and monitoring required for their Canadian Nuclear Safety Commission and Province of Saskatchewan licenses. The sampling period included the years 2002 through 2009. All samples were collected using a standard Ekman Grab. Samples consisted of multiple grabs (three, five or 10 subsamples) pooled to produce a composite sample.

Benthic samples were collected in March, April, June, September and October. Samples from March and April were classified as winter collections as they were collected through the ice; June samples were classified as spring; and September and October samples were classified as fall. Samples were sieved using 500-μm mesh, and then preserved with formaldehyde. Benthic samples were processed by the same taxonomist (Dr. Jack Zloty) and taxa were identified to the lowest practical taxonomic levels according to MMER EEM guidance (Environment Canada 2012). Total invertebrate density (# organisms/m²) was calculated as well as conventional indices of benthic invertebrate community composition (taxon richness, Simpson’s Diversity Index) calculated at the lowest practical level of taxonomic identification. Index computations followed MMER EEM guidance (Environment Canada 2012). In addition to the indices described above, correspondence analysis (CA) was used to summarize (ordinate) variations in community composition among stations. CA was completed with benthic invertebrates identified to the lowest practical level. Taxon abundances were log₁₀(x + 1) transformed prior to analysis. Non-rare taxa (i.e. those accounting for at least 0.5% of the total number of benthic invertebrates in the data set, and occurring in at least 10% of the samples from all reference and exposure lakes) were included in the analysis. The CA was conducted using an Excel add-in (Biplot 1.1 2002).

Template variables (Stanfield and Kilgour 2006) included easting, northing, lake area (ha), lake perimeter (m), year, average lake depth, maximum lake depth and depth (m) of sample collection. Sediments were collected concurrently with most benthic samples. Sediments were collected using Tech-Ops Extruded corers, with the top 2 cm collected and analyzed for total metals (microwave-assisted acid digestion process, EPA Method 3052) and grain size at the Saskatchewan Research Council. Texture of sediments was described as percent coarse sand, percent fine sand, percent silt, percent clay, and loss on ignition (LOI), which were also considered potential template variables. Principal component analysis (PCA) was used to summarize the variation in site conditions using the template variables. Box and probability plots were used to explore the distributions of potential template variables. Lake area, lake perimeter, average lake depth, percent coarse sand, percent fine sand, percent silt, LOI and sample depth were non-normally distributed, and thus log₁₀-transformed prior to carrying out the PCA. PCA axes were retained for consideration when they accounted for >10% of the total variance.

To address the issue of low abundance in these oligotrophic aquatic systems, site study designs included compositing multiple grabs (three, five or 10 grabs depending on study site) to make a single composite sample. The number of grabs pooled to make a composite sample was explored as an additional potential template variable because of the general influence sampling effort has on estimates of taxon richness and diversity.

The year that samples were collected was also a significant predictor of the various indices of benthic invertebrate community composition. Variations from year to year may be natural (i.e. regional climatic variations or spatial heterogeneity), or mine influenced (i.e. changes in effluent volume and/or loadings), although low sedimentation rates in this region mean changes in sediment contaminant concentrations are likely not significant over this 7-year time period. Assuming that climatic variation is somewhat random, the yearly variations were ignored and allowed to contribute to the unexplained variation. Doing so has the potential to somewhat increase the ‘normal range of variation’. Further work could assess the cause–effect relationship between year and variations in indices of benthic invertebrate community composition (e.g. average annual air temperature, average annual precipitation). The causes of the annual variations in benthic invertebrate community composition were considered beyond the scope of this present investigation. As more years of data are added to this data set from Saskatchewan, exploration of causes of annual variations such as they may relate to climatic variables would be beneficial considering the scant literature on this topic (Baulch et al. 2005).

Reference-sample models produced by a general linear model (GLM) procedure using SYSTAT 13.0 (Crane Software International Ltd.) were used to estimate ‘expected’ reference-condition values for all individual samples, including those from both exposure and reference lakes. The difference between the expected and observed values is a ‘residual’. Residuals¹ were standardized to a Z score by dividing the difference by the square root of the mean square error (MSE) term (Kilgour and Stanfield 2006). Residuals that were > ± 2 standard deviations (SDs) from an expected mean value are beyond the range of values for about 95% of reference samples, while samples that are >± 3 SDs from the expected mean are outside the 99% region (Kilgour et al. 1998). Samples producing standardized residuals outside ± 2 SD for the reference models are considered here, like elsewhere (e.g. Kilgour et al. 1998, 2017), to be outside the normal range of reference values.

The logarithm of lake area and logarithm of percent silt were included as template variables on the basis of the PCA that was conducted in the previous step. The logarithm of sample depth was included because it is a known template variable influencing the distribution of benthic invertebrates in lakes (Stanfield and Kilgour 2006) and on the basis of the PCA in the previous step. Sampling year (2002 through 2009) and season (winter, spring, fall) were included as classification variables due to their potential influence on community composition.

SYSTAT’s GLM procedure was used to construct models that predicted each of three EEM endpoints (abundance, taxon richness, Simpsons Diversity Index). The seven template variables were initially included in the model, and were removed in a backward stepwise procedure. Variables were retained if they explained enough variation to be significant at p < 0.05. The robustness of model parameter estimates was evaluated by jackknife-resampling the data and generating a distribution of model parameter estimates. The lower 2.5th and upper 97.5th percentiles of generated parameter estimates were determined and used as an estimate of confidence regions for model parameters. Box plots of residuals were produced to visualize variability in ‘deviations’ from the normal expected data ranges. Box plots included horizontal lines at ± 2 SDs, to illustrate limits of the normal ranges (Kilgour et al. 1998) for residuals for each index of community composition.

Sediment samples were collected within the same general location (stations) as benthic invertebrate community samples. Benthic invertebrate community residuals were averaged by station, as were concentrations of the various measured sediment chemistry parameters (total metals and radionuclides, total organic carbon). Sediment chemistry and radionuclide activity data were log-transformed after being averaged, to increase the potential linearity of association with indices of benthic invertebrate community composition. Spearman rank correlations between sediment variables and indices of benthic invertebrate community composition were computed.

Analysis of Covariance (ANCOVA) was used to determine whether the nature of the observed associations between elements and taxon richness residuals differed between lakes classified as reference and exposure, by testing whether the slope of the negative relationship between element concentration or activity and taxon richness residuals differed for reference and exposure lakes. ANCOVA focused on taxon richness residuals because they tended to produce the strongest and most consistent correlations with metals/radionuclides. The ANCOVA models included a classification variable (Condition: i.e. reference or exposed), the Element Concentration or Activity (e.g. selenium, uranium, radium-226), and the interaction term (Condition × Element Concentration or Activity). The interaction term tested for differences in slopes of the relation between concentration/activity and taxon richness residuals between the two groups, reference and exposure.

Results

The database contained 21 lakes in a reference condition, unexposed to mining/milling operations, and 29 lakes that were exposed to influences from mining/milling operations. Reference lakes varied in size (area) from 1 to 600+ ha; exposure lakes varied from 1 to 2099 ha. Lake depth (average and maximum) and perimeter co-varied with lake area (see discussion of PCA below), with several of the exposure lake study areas (i.e. Hidden Bay and Pow Bay from the Rabbit Lake Operation) being larger than the reference lakes. Benthos sample depths varied between ~1 m and almost 20 m in both reference and exposure lakes, although sample depths were matched within individual mine site study designs. All sediments consisted of fine-grained materials (coarse sands through to silt). LOI of sediments varied considerably, from 0 to 60+%, in both reference and exposure lakes (Table 1).

Table 1. Variation in template variable conditions in reference and exposure lakes in Northern Saskatchewan from which benthic macroinvertebrates were collected between 2002 and 2009. Number of samples refers to the number of lakes (first column) or the number of samples among those lakes (subsequent columns).

Condition	Statistic	Lake area (ha)	Average lake depth (m)	Maximum lake depth (m)	Perimeter (m)	Sample depth (m)	Coarse sand (%)	Fine sand (%)	Silt (%)	Loss on ignition (%)
Reference	Number of samples	21	21	21	21	212	173	173	173	71
	Minimum	2	0.84	1.5	790	0.8	0	0	0	0
	Maximum	632	6	21	41570	19	55	63	78	63
	Average	158	3	7	12363	4	4	16	45	37
	Standard deviation	179	2	6	11844	3	9	17	16	15
Exposure	Number of samples	29	29	29	29	243	113	113	112	103
	Minimum	1	1	1	596	0.4	0	0	16	4
	Maximum	2697	12	42	65500	18	52	59	91	72
	Average	377	5	10	12019	5	4	13	58	22
	Standard deviation	734	4	10	17785	5	9	16	16	13

The first PCA axis described a gradient of lake size (Table 2). Based on the nature of the variable correlations with the axis scores, larger lakes tended to have greater maximum and average lake depth, be sampled in deeper water, have more coarse or fine sand and less silt or clay, and have less organic material in the sediments (i.e. LOI). The second PCA axis reflected the geographic location of the lakes, with the reference and exposure lakes generally being distributed along an axis from northeast to southwest. Average lake depth and northing were selected as representative template variables for model construction on the basis that they were the most strongly associated with independent PCA axes, and thus will generally and potentially represent different aspects of natural environment-induced variation.

Table 2. Results of principal component analysis (PCA) of reference and exposure lake template variables. Data are from 50 lakes in Northern Saskatchewan associated with benthic community samples collected between 2002 and 2009. Values shown are Pearson correlations (r values) between the measured variables and the PCA axis scores. Pearson r values > |0.6| are bolded to indicate a strong association between the variable and the axis. Some variables were log₁₀ transformed, as indicated in the table.

Variable	Axis
	1	2
Easting (m)	0.25	−0.95
Northing (m)	0.06	−0.97
Log10 of lake area (ha)	0.86	0.12
Maximum lake depth (m)	0.93	−0.21
Log10 of sample depth (m)	0.84	0.06
Log10 of average lake depth (m)	0.95	−0.12
Lake perimeter (m)	0.81	0.29
Log10 of coarse sand (%)	0.66	0.4
Log10 of fine sand (%)	0.73	0.46
Clay (%)	−0.84	0.03
Log10 of loss on ignition (%)	−0.76	0.46
Variance explained (%)	56	23

Community composition

The combined data sets across years and studies produced a database with 262 samples from reference lakes unexposed to mining/milling, and 328 samples from lakes exposed to the potential influences of mining/milling (i.e. exposure lakes). Total abundances varied from lows of ≤ 100 individuals per m², to upward of 20,000 per m² in reference lakes and 100,000 per m² in exposure lakes. The number of taxa (at the lowest practical taxonomic level, i.e. taxon richness) per sample varied similarly in reference and exposure samples between ~1–2 and 30–40 per sample. Simpson’s Diversity Index also varied similarly in reference and exposure samples (Table 3).

Table 3. Variations in reference and exposure indices of benthic macroinvertebrate community composition in 21 reference lakes and 29 lakes exposed to uranium mining and milling in Northern Saskatchewan sampled between 2002 and 2009.

Condition	Statistic	Total Abundance (#/m^2)	Taxon richness	Simpson’s Diversity
Reference	Minimum	36	2	0.39
	Maximum	21,480	34	0.92
	Arithmetic mean	2629	15	0.73
	Standard deviation	3158	6	0.12
Exposure	Minimum	44	1	0
	Maximum	103,564	44	0.91
	Arithmetic mean	5622	14	0.67
	Standard deviation	11,280	6	0.17

The dominant benthic taxa used in the multivariate ordination (Figure 2) included a variety of chironomids, fingernail clams (Sphaeriidae), gastropod snails, amphipods (Hyalella azteca, Diporeia hoyi), ceratopogonids (Probezzia), nematodes, tubificid worms, leeches (Helobdella stagnalis), phantom midges (Chaoborus) and water mites (Hydracarina). None of the taxa used in the ordination were unusual or otherwise unexpected; all are typical inhabitants of boreal lake systems.

Figure 2. Plot of scores on CA Axis 1 and CA Axis 2 for taxa. The scores are from a correspondence analysis of logarithms of abundance of non-rare benthos (as defined in the text) from sediments in 50 lakes in Northern Saskatchewan varying in exposure to uranium mining/milling, between 2002 and 2009.

Reference models

The associations between benthic indices of composition and lake area, season, year and sample depth are illustrated in Figures 3 to 7. All model parameters were confirmed by jackknife resampling to be non-zero. The logarithm of lake area explained significant variation in each of the three core indices of benthic invertebrate community composition, as well as the two CA axes (Table 4). Total abundance, taxon richness (number of taxa) and Simpson’s Diversity Index generally increased with lake area and decreased with sample depth. Taxon richness and diversity were generally lower when more benthic invertebrate community grabs were pooled to make a composite sample. Reference models for CA axes reflected the general tendencies anticipated for specific taxa in relation to water depth (Table 4). For both ordination axes, coefficients for lake area were much smaller than the coefficients for sample depth, reflecting the much greater influence of sample depth on community composition. Samples collected from deeper water tended to have negative CA Axis 1 scores and negative CA Axis 2 scores, reflecting a higher relative abundance of phantom midges (Chaoborus), the midge Chironomus, Tubificidae and the deep-water amphipod Diporeia hoyi. Samples collected from shallower water tended to produce more positive CA Axis 1 scores, reflecting a greater diversity of benthic fauna, and including higher relative abundances of a variety of chironomids (including Pseudochironomus, Microtendipes, Heterotanytarsus, Dicrotendipes, etc.), the shallow-water amphipod Hyalella azteca and the snail Valvata sincera (Figure 2).

Figure 3. Scatter and box plots of total abundance (#/m²) of benthic organisms in 50 lakes in Northern Saskatchewan between 2002 and 2009, in relation to lake size (area in ha), sampling season, year of sample and sample depth (m).

Figure 4. Scatter and box plots of taxon richness (number of taxa at the family level) in 50 lakes in Northern Saskatchewan between 2002 and 2009, in relation to lake size (area in ha), sampling season, year of sample and sample depth (m).

Figure 5. Scatter and box plots of Simpson’s Diversity in 50 lakes in Northern Saskatchewan between 2002 and 2009, in relation to lake size (area in ha), sampling season, year of sample and sample depth (m).

Figure 6. Scatter and box plots of CA Axis 1 scores in 50 lakes in Northern Saskatchewan between 2002 and 2009, in relation to lake size (area in ha), sampling season, year of sample and sample depth (m). Scores were derived from a correspondence analysis ordination of logarithms of abundances of non-rare taxa (as defined in the text).

Figure 7. Scatter and box plots of CA Axis 2 scores in 50 lakes in Northern Saskatchewan between 2002 and 2009, in relation to lake size (area in ha), sampling season, year of sample and sample depth (m). Scores were derived from a correspondence analysis ordination of logarithms of abundances of non-rare taxa (as defined in the text).

Table 4. Results of reference-condition model for benthic community indices of composition for samples in lakes in Northern Saskatchewan. Models are based on benthic samples collected between 2002 and 2009 from 21 reference lakes. Model coefficients are provided, as are descriptors of model fit (R² = variance explained; MSE = residual mean squared error) and sample size (N = number of samples in the model).

Variable	Total abundance	Taxon richness	Simpson’s Diversity	CA Axis 1 scores	CA Axis 2 scores
Constant	3.17	20.61	0.789	2.157	0.258
Log10 of lake area (ha)	0.27	2.092	0.032	0.25	0.444
Log10 of sample depth (m)	−0.736	−12.45	−0.156	−3.451	−2.81
Fall	0.092	−0.69	–	0.37	−0.495
Winter	0.108	1.922	–	0.053	0.822
Spring	−0.198	−1.242	–	−0.422	−0.327
# Grabs/sample	–	−0.608	−0.006	−0.155	–
R^2	0.161	0.345	0.127	0.796	0.472
MSE	0.229	23.96	0.009	0.229	0.633
N	200	202	203	203	203

Note: In all cases, ‘jackknifing’ analysis confirmed that coefficients were non-zero.

Assessment of exposure lakes

A box plot of residuals for each index or CA axis was produced for all exposure samples, combining all years, for illustration purposes (Figure 8). Residuals between −2 SD and +2 SD were community index values from exposure lakes that fell within the ‘normal range of variation’ (or approximate 95% region) for predicted reference samples. Residuals <−2 SD were the result of index values that were less than the normal range, while residuals >+2 SD were the result of index values that were greater than the normal range. Of the three EEM metrics, taxon richness was the most consistent in producing large negative residuals, implying that many of the samples from exposure lakes had fewer taxa than the reference model. Total benthic invertebrate community abundance in samples from exposure lakes was rarely above or below the expected range of values for a reference condition, with the exception of samples from East Boomerang Lake, where abundances were at times high. Abundances in samples from McDonald Lake and Vulture Lake were also frequently lower than the predicted reference condition. Simpson’s Diversity Index values in samples from exposure lakes were never higher than the expected reference condition, and were frequently lower than the expected range of reference values, particularly in samples from East Boomerang Lake, Sink Lake and Fox Lake.

Figure 8. Box plot of indices of benthic community composition for each of the lakes exposed to uranium mining effluent from each of the five major uranium mining and milling projects in Saskatchewan: Rabbit = Rabbit Lake Operations; Key = Key Lake Operations; ML = McClean Lake Operations; MR = McArthur River Operations; CL = Cigar Lake Operations. Lakes within project areas are: 1 = Horseshoe Pond; 2 = shallow basin of Hidden Bay; 3 = deep basin of Hidden Bay; 4 = Upper Basin of Upper Link Lake, 5 = Lower Basin of Upper Link Lake; 6 = Lower Link Lake; 7 = Pow Bay; 8 = Fox Lake; 9 = Unknown Lake; 10 = Delta Lake; 11 = Little McDonald Lake; 12 = McDonald Lake, 13 = Sink Reservoir; 14 = Vulture Lake; 15 = East McClean Lake; 16 = Kewen Lake; 17 = East Boomerang Lake; 18 = Little Yallowega Lake; 19 = Aline Lake; and 20 = Seru Bay. Gray shading indicates the regional reference data defined by ± 2 standard deviations. ‘CA Axis 1’ and ‘CA Axis 2’ are scores on axes 1 and 2 from a correspondence analysis of logarithms of abundance of non-rare benthos (as defined in the text) from sediments in 50 lakes in Northern Saskatchewan varying in exposure to uranium mining/milling, between 2002 and 2009.

In contrast to the conventional core indices of community composition, the residuals based on multivariate CA Axes 1 and 2 were more variable for samples from exposure lakes than for samples from reference lakes. In this analysis, large positive residuals on CA Axis 1 implied greater benthic diversity and higher relative abundances of benthic taxa considered more sensitive. Samples producing large negative residuals on CA Axis 1 in contrast could be considered to have had lower benthic diversity and more taxa that were tolerant (e.g. Chironomus). Fox Lake, Horseshoe Pond, Little McDonald Lake and Lower Link Lake each produced a number of residuals on CA Axis 1 that were large and negative, indicative that the benthic communities of those lakes were less diverse, had lower abundances of taxa that are indicative of good water and sediment quality, and had higher abundances of taxa that are tolerant of impaired water and sediment conditions when compared to the normal range of reference conditions.

Correlations between benthos and sediment chemistry

Relationships between benthic index residuals (Z scores) and sediment contaminant concentrations were studied as a means of examining additional sources of variation in indices of benthic community composition. This further assessment amounts to preliminary ‘investigation of cause’ (per the MMER methodology, Environment Canada 2012) of the observed deviations from a ‘normal’ reference condition. Of particular interest were strong negative Spearman rank correlations that implied a decrease in benthic invertebrate community condition with increasing elemental concentration or activity, particularly when the range of concentrations of a measured element include values that exceed a sediment quality guideline. A critical Spearman rank correlation (r_s) value of |0.3| was selected on the basis that there is a low probability of random occurrence (i.e. p~0.001 for n = 72), and correlations of that magnitude are generally visibly apparent. Nineteen of the 120 correlations (i.e. ~16%) were greater than the critical value (r_s > |0.3|), and in a direction that could possibly indicate an impaired benthic community condition with increasing elemental concentration or activity (Table 5).

Table 5. Sediment variable correlations with benthic macroinvertebrate community indices of composition. Data are from 21 reference lakes and 29 lakes exposed to uranium mining and milling in Northern Saskatchewan between 2002 and 2009. Spearman rank correlations (r_S) > |0.3| are in bold, indicating a strong and potentially causal relationship (e.g. when negative for taxon richness or Simpson’s Diversity).

Parameter	Low Effect Level (LEL) Thompson et al. (2005)	Interim Sediment Quality Guideline (ISQG) CCME (2015)	Probable Effect Level (PEL) CCME (2015)	Statistic			Spearman rank correlation (rS) with benthic community index (Z score)
				Minimum	Maximum	Mean	Total abundance	Taxon richness	Simpson’s Diversity	CA Axis 1 scores	CA Axis 2 scores
Antimony				0.2	67.0	0.2	0.29	−0.24	−0.17	−0.22	−0.09
Cobalt							−0.08	−0.17	−0.04	−0.05	0.00
Arsenic	9.8	5.9	17	0.5	2000	12.6	−0.09	−0.35	−0.12	−0.21	−0.11
Aluminum				5.5	41900	6970	−0.07	0.04	0.11	0.27	0.17
Nickel	23.4			0.2	1300	18	−0.06	−0.39	−0.21	−0.17	−0.05
Copper	22.2	35.7	197	0.5	92.0	6.5	−0.13	−0.29	−0.13	−0.05	0.02
Beryllium				0.1	3.0	0.4	0.00	0.18	0.23	0.10	0.04
Barium				1	2300	89	−0.02	−0.26	−0.10	0.14	0.20
Radium-226	0.6			0.01	13	0.09	−0.11	−0.37	−0.13	−0.05	0.01
Thorium-230				0.02	400	0.05	−0.16	−0.07	0.08	0.36	0.22
Zinc		123	315	0.5	480	44.3	−0.18	0.08	0.13	0.09	−0.02
Chromium	47.6			2	1280	16	0.20	0.07	0.03	0.07	0.05
Iron		20,000	40,000	0	332,001	30,609	0.06	0.17	0.21	0.17	0.16
Molybdenum	13.8			0.2	6990	17.0	−0.08	−0.47	−0.15	−0.19	−0.12
Uranium	104.4			0.25	7100	12.2	−0.22	−0.62	−0.27	−0.11	−0.14
Cadmium		0.6	3.5	0.1	22.0	0.5	0.04	−0.07	0.07	−0.02	0.10
Lead	36.7	35	91.3	0.3	86.0	7.1	−0.08	−0.11	−0.10	0.02	0.03
Selenium	1.9			0.1	243	1.4	−0.23	−0.63	−0.32	−0.48	−0.30
Manganese				15.0	16100	403	−0.11	0.18	0.26	0.52	0.36
Boron				1	126.0	4	−0.03	−0.01	−0.03	−0.29	−0.17
Polonium-210				0.02	22.00	0.48	−0.15	−0.46	−0.26	−0.39	−0.28
Loss on ignition				2.0	220	27.5	0.03	−0.24	−0.15	−0.28	−0.23

Abundance residuals did not strongly correlate with any measured elements. Taxon richness residuals were strongly correlated with about half of the measured elements, and usually negatively, such that taxon richness was lower in lakes that had sediments with higher concentrations of measured metals or activity of radionuclides. Taxon richness residuals were most strongly negatively associated with selenium (r_s = −0.63) and uranium (r_s =_-0.62) (Table 5). Concentrations of nickel, arsenic, molybdenum and strontium, and activity of radium-226, polonium-210 and lead-210 were also strongly (and negatively) associated with taxon richness residuals, indicating potential influences of one or more of those elements on benthic invertebrate community composition and the correlation of these elements with each other.

ANCOVA was used to test whether the slope of the negative relationship between element concentrations or activity and taxon richness residuals differed for reference and exposure lakes. The interaction term tested for differences in slopes of the relation between concentration/activity and taxon richness residuals between the two groups, reference and exposure. The p values for the interaction terms were all >0.05, indicating no differences in slopes for reference and exposure areas, or more generally that the relationships between element concentration or activity and taxon richness residuals were similar in reference and exposure lakes. The statistical result (i.e. no differences in slopes) was consistent with what was observed through visual inspection of the data (Figure 9) (i.e. the data for both reference and exposure lakes tended to be well overlapped). Therefore, benthic community taxon richness at natural background concentrations in the Athabasca Basin was negatively influenced by increases in sediment metal concentration or radionuclide activity; however, element concentrations and activities in exposure lakes tended to be higher than in reference lakes.

Figure 9. Scatterplots of benthic taxon richness residuals in relation to various elements with strong associations with richness residuals (see Table 2) for samples from 50 lakes in Northern Saskatchewan between 2002 and 2009.

Residuals for Simpson’s Diversity Index and for scores on CA Axes 1 and 2 did correlate strongly with some element concentrations. Correlations for CA Axis 1 and Axis 2 residuals that were negative were considered evidence of a potential negative influence. Tubificidae and Chironomus had lower negative CA Axis 1 scores. Therefore, a negative relationship between CA Axis 1 residuals and an element would imply that higher concentrations or activities of the element potentially covaried with higher relative abundances of tubificids and/or Chironomus (Figure 2). Generally, the elements that correlated strongly with taxon richness residuals were also strongly (and negatively) associated with residuals for CA Axes 1 and 2 (Table 5).

Discussion

This analysis of lake benthic invertebrate data from Northern Saskatchewan developed reference condition models for assessing the condition of benthic invertebrate communities against regional reference data. General linear models were developed that predicted values for indices of benthic invertebrate community composition on the basis of template variables including season, sampling depth and various measures of lake size. These models have the potential to support the assessment of the effects of uranium mining and milling on aquatic resources in the north of Saskatchewan. Operations that have difficulty finding suitable reference lakes can use these models to predict benthic community indices of composition for the lakes that are exposed to their mining/milling operations. Further, mining/milling operations that have suitable reference lakes can still use the regional models to more fully describe reference variability that can be used to interpret local variations in community composition (Arciszewski and Munkittrick 2015; Kilgour et al. 2017).

There are four additional observations from this analysis that deserve further discussion. First, there were some differences in lake size between reference and exposure lakes, making it impossible to consummately assess benthic communities for some exposure lakes. Second, the associations between richness and metals were developed after adjusting richness for a variety of natural factors, and should be robust. Third, the models used in this analysis did not explore and take into account potential sources of natural annual or year-to-year variations. If climatic conditions, for example, had been taken into account, the reference models might have accounted for more of the background reference variability, and the assessment of mining/milling might have resulted in more significant effects in lakes exposed to mining/milling. Fourth, the various conventional core EEM indices of composition were about as sensitive as the multivariate measures developed using correspondence analysis to exposure to mining/milling operations, while taxon richness (at the lowest practical level) was very sensitive. Fifth, the negative correlations between residuals of indices of benthic invertebrate community composition and several of the metal concentrations and radionuclide activities were consistent with what might be expected on the basis of published sediment quality guidelines.

Abundance, taxon richness and Simpson’s Diversity Index generally decreased with lake size (Figures 3 to 5). The use of lake size as a template predictor variable in the RCA model was, therefore, well supported. Without reference lakes larger than ~600 ha, however, it is impossible to predict whether those decreases in abundance with lake size would continue; and, further, it is not possible to produce a rigorous statistical assessment of the condition of the larger lakes. This analysis did not explore the potential influences of climatic variables, such as mean annual air temperature, precipitation, etc., which can be anticipated to have some degree of influence on annual variations in indices of benthic invertebrate community composition (Ward and Stanford 1982; Hogg and Williams 1996). Future reference condition models could use indices of climatic variation, but it would involve a research exercise to identify the important key variables. Alternatively, and as was the case here, the climatic variations can be assumed to be part of the natural variations that inevitably occur, and could then be used to assist in setting benchmarks for indices of composition; for example, the observed variations naturally occur, and are tolerable by the benthic invertebrate community, and so should be tolerable in a lake exposed to treated uranium mining/milling effluents.

This analysis focused on regulated MMER variables, in addition to axes from CA. Other variables could have been selected, and might have produced different reference models and differences in sensitivity of variables to mining/milling. Multivariate measures of benthic invertebrate community composition are generally the most sensitive indicators of change, and will produce more striking statistical effects compared to more conventional indices (Barton 1996; Kilgour et al. 2004). In this analysis of exposure lakes, total benthic abundances were rarely high or low relative to what was observed in reference lakes, while most exposure lakes (with four exceptions) had at least some samples with low taxon richness compared to the normal range of variations in reference lakes. The multivariate ordination was not as ‘sensitive’ to exposure to mining/milling as was taxon richness, which may in part be due to the ordination process. Only taxa accounting for at least 0.5% of the total number of organisms, and found in at least 10% of samples, were retained in the analysis. These criteria resulted in the CA ordination being carried out with 28 of 187 possible taxa. Clearly, then, the analysis focused only on those taxa that were broadly distributed in the study area, and might not have reflected potential subtle variations in exposure lakes caused by some of the rarer taxa. The data reduction that was used in this analysis was conventional and not out of the ordinary (Gauch 1982). That a simpler metric like taxon richness was more sensitive in this case study than the multivariate indices is not troubling, but does illustrate the necessity of a multi-faceted approach to the analysis of benthic invertebrate community data in order to understand the nature and ecological significance of effects. Now that the data are available in one repository, any number of different approaches (Bailery et al. 2014) to the development of reference models and assessment of mining/milling influences can be carried out for this somewhat unique sector.

The results indicate that the composition of some invertebrate communities has changed due to exposure to uranium mining activities in Northern Saskatchewan. The number of benthic invertebrate taxa downstream at approximately half of the uranium mine/mill effluent release points was lower compared to the number of taxa observed in reference conditions. That observed effect was not consistent across the set of receiving environments in this analysis, and clearly reflects localized effects and conditions. Benthic taxon richness residuals were strongly correlated with about half of the measured sediment elements and always negatively, such that taxon richness was lower in lakes with higher concentrations/activities of measured metals and radionuclides in the sediments. Taxon richness residuals were most strongly associated with selenium and uranium concentrations, both of which have been studied within these systems with accumulated weight of evidence identifying them as contributors to observed aquatic effects (Robertson and Liber 2007; Muscatello et al. 2008; Muscatello and Liber 2009; Wiramanden et al. 2010; Gallego-Gallegos et al. 2012; Tse et al. 2012; Franz et al. 2013). The general decreasing trend in taxon richness residuals with increasing metal concentration or radionuclide activity was consistent for both reference and exposure lakes, implying that there is a natural background influence of metal concentrations and radionuclide activity on benthic taxon richness. Further, the association between richness and metals and radionuclides was found after accounting for variations in richness related to water depth and associated substrate texture and organic carbon content. The association between richness and metals/radionuclides, therefore, is considered reasonably robust in this analysis. These data could be further explored to identify limits to natural background concentrations of inorganic elements in sediments, or the concentrations/activities that correspond with large variations in benthic communities.

Conclusions

The results indicate the value of benthic invertebrate community monitoring programs supported by comprehensive monitoring of sediment chemistry. Taxon richness was the most sensitive benthic metric for monitoring benthic community responses downstream of uranium mine and mill effluents. The application of the reference condition approach for identifying impacted benthic communities that differ significantly from natural regional compositions should be encouraged in programs designed to meet regulatory requirements.

Acknowledgements

Cameco and Areva provided the data for independent analysis and interpretation by the Canadian Nuclear Safety Commission (CNSC) and reviewed an earlier version of the paper. The work of developing the reference models was funded by the CNSC under contract to Kilgour & Associates. All data checking, analysis and interpretation were conducted independently of Cameco and Areva. Comments on the manuscript by Keith Somers are appreciated. We thank Catherine Proulx for translating the abstract, and Dianne Drummond for help with figures.

Notes

1. Hereafter, the term ‘residual’ is used to refer to standardized residuals, or Z scores.