Wetland macrophyte community response to and recovery from direct application of glyphosate-based herbicides

Ecotoxicology and Environmental Safety 183 (2019) 109475

Contents lists available at ScienceDirect

Ecotoxicology and Environmental Safety journal homepage: www.elsevier.com/locate/ecoenv

Wetland macrophyte community response to and recovery from direct application of glyphosate-based herbicides

T

Joseph F. Mudgea, Jeffrey E. Houlahanb,∗ a b

Cooperators General Insurance Company, Guelph, Ontario, N1H 6P8, Canada Department of Biological Sciences, The Canadian Rivers Institute, 100 Tucker Park Road, University of New Brunswick, Saint John, New Brunswick, E2L 4L5, Canada

A R T I C LE I N FO

A B S T R A C T

Keywords: Glyphosate Optimal alpha Macrophytes Wetlands Recovery Critical effect size

Community-scale impacts of glyphosate-based herbicides on wetland plant communities and the magnitude of those impacts that should be considered biologically relevant are poorly understood. We contrast three different thresholds for setting biologically meaningful critical effect sizes for complex ANOVA study designs. We use each of the of the critical effect sizes to determine optimal α levels for assessment of how different concentrations of glyphosate-based herbicides affect wetland plant communities over two years of herbicide application (alone and in combination with agricultural fertilizers) and two subsequent years without herbicide (or fertilizer) application. The application of glyphosate-based herbicides was found to result in a decrease in macrophyte species richness, an increase in macrophyte species evenness, a decrease in macrophyte cover and a reduction in community similarity. There was little evidence that nutrient additions directly or indirectly affected plant community endpoints. The glyphosate effects were evident in the first year of herbicide application in 2009, and became more pronounced in the second year of herbicide application in 2010. However, when herbicides were not applied in 2011, recovery was observed in most endpoints, with the exception being species evenness, for which partial recovery was not observed until 2012. Optimal α levels differed among the three critical effect sizes for each ANOVA term and endpoint combination, however regardless of differences in α levels, conclusions were generally consistent across all critical effect sizes.

1. Introduction Glyphosate-based herbicides (N-phosphonomethyl glycine) are nonselective, systemic herbicides used for weed control in domestic, agricultural, and silvicultural contexts to reduce competing vegetation. Glyphosate inhibits plant growth by disrupting the enzyme enolpyruvylshikimate phosphate synthase (EPSPS) pathway, which is required to make several amino acids (Rubin et al., 1982; Williams et al., 2000). One reason glyphosate-based herbicides are the most widely used in the world (Woodburn, 2000; Grube et al., 2011) is because the amino acid biosynthesis pathway they target is exclusive to plants, apicomplexans, ascomycete fungi and prokaryotes, and is not shared by any members of the animal kingdom (Williams et al., 2000). Most studies concerning the effects of glyphosate-based herbicides on wetland macrophytes have been focused on controlling invasive wetland plant species such as the common reed (Phragmites australis) (Mozdzer et al., 2008; Crowe et al., 2011), water hyacinth (Eichhornia crassipes) (Jadhav et al., 2008) or cattail (Typha spp.) (Linz and Homan, 2011). Evaluations of the impacts of glyphosate-based herbicides on

*

non-target wetland macrophyte species are not common, but speciesspecific non-target negative impacts have been reported for duckweed (Lemna spp.) (Gardner and Grue, 1996; Kielak et al., 2011), Bolboschoenus maritimus (Mateos-Naranjo and Perez-Martin, 2013). Thus, unintended effects of glyphosate application could impact nutrient cycling because duckweed has been demonstrated to contribute to nutrient removal in wastewater (Zhou et al., 2018, 2019). Examinations of the effects of glyphosate-based herbicides on multiple non-target species are even less common. Dalton and Boutin (2010) have reported variable toxicity among a mixture of seven target and non-target wetland macrophyte species (Asclepias incarnata, Chelone glabra, Eupatorium maculatum, Eupatorium perfoliatum, Lycopus americanus, Phalaris arundinacea, and Verbena hastata) grown in a greenhouse and outdoor microcosm environments. They concluded that the impacts of glyphosate herbicides on wetland macrophyte species are variable among species, application times and other application procedures (Dalton and Boutin, 2010). However, the aforementioned multi-species research in greenhouse and outdoor microcosms has yet to be validated in natural systems.

Corresponding author. E-mail address: jeff[email protected] (J.E. Houlahan).

https://doi.org/10.1016/j.ecoenv.2019.109475 Received 16 May 2019; Received in revised form 18 July 2019; Accepted 24 July 2019 Available online 20 August 2019 0147-6513/ © 2019 Elsevier Inc. All rights reserved.

Ecotoxicology and Environmental Safety 183 (2019) 109475

J.F. Mudge and J.E. Houlahan

become naturalized to varying extent since their creation. Within this area we chose 24 wetlands intended to be representative of wetlands located in and around agricultural fields, that were relatively small (< 1 ha), had no permanent inflow or outflow, and had relatively homogenous macrophyte cover.

The critical effect size that would constitute a biological relevant difference in the wetland macrophyte community following glyphosate herbicide application can be conceptualized in different ways. Researchers using different critical effect size measures may dispute the biological relevance of results even though they agree on the statistical significance of the results using α = 0.05. This disconnect between statistical significance and biological relevance is a well-known flaw of null Hypothesis testing using α = 0.05 (Anderson et al., 2000). This disconnect can be resolved by incorporating critical effect sizes into the calculation of optimal α levels that both tie statistical significance to biological relevance while also minimizing the probabilities and/or costs of Type I and biologically relevant Type II errors of null hypothesis tests (Mudge et al., 2012). Our objectives are to contrast three potential approaches for setting a biologically meaningful critical effect size, using each potential critical effect size to determine optimal significance levels (Mudge et al., 2012) for assessment of how different concentrations of glyphosatebased herbicides affect wetland plant communities over two years of herbicide application (alone and in combination with agricultural fertilizers) and two subsequent years without herbicide (or fertilizer) application. This will allow us to: (1) apply the optimal α approach on real data, (2) illustrate how different approaches to setting a critical effect size can be used to make appropriate inferences, (3) evaluate the effects of glyphosate based herbicides on wetland plant communities in the years the herbicides were applied, (4) assess the interactive effect of herbicide and nutrient applications, (5) determine effects of varying glyphosate herbicide concentration, (6) test for lingering effects in the years after herbicides were applied (i.e. recovery).

2.2. Experimental design The study was designed as a replicated split-plot experiment. Each experimental wetland was bisected by an impermeable plastic barrier and one half of the system was randomly assigned as a control while the other half was assigned the treatment. Bisections were oriented along the shortest possible diameter through the wetland with similar topography, habitat features and plant community composition on each side. In two wetlands with suspected transient flow patterns (usually only in spring), the inflow side was designated control and the outflow side was treated. This prevented contamination of the control side of the wetland if the integrity of the barrier was compromised during the spring because the control side water flowed into the treatment side rather than the reverse. This split-plot design, where each wetland is partitioned into a control and an experimental area, minimizes the problem of confounding the natural variability among ponds with treatment effects. During the late summer of 2008 each of the 24 wetlands was split in half using an impermeable, stable, inert, UV resistant plastic barrier constructed from 0.76 mm opaque black high density polyethylene (HDPE) (Poly-Flex Inc. Geomembrane Lining Systems, Grand Prairie, Texas). The bottom of the barrier was weighted with gravel and buried 10 cm into the bottom substrate and extended approximately 2 m beyond the high water mark of each wetland. The barriers were approximately 1 m in height to ensure that they extended above the high water level for every wetland. The integrity of the barrier was inspected in the spring and fall and periodically throughout all field seasons. Barrier maintenance was required in all wetlands at some point during the study, but there was no evidence of contamination with glyphosate levels consistently below detection limits in water samples from control halves (Edge et al., 2012). The glyphosate-based herbicide Roundup WeatherMax™ (540 g acid equivalent (a.e.)/L, Monsanto, Winnipeg, MB, CAN) was applied at two different concentrations, alone and in combination with nutrients. This resulted in four treatment categories:

Hypothesis 1a/b/c. Glyphosate reduces plant cover with effect sizes of (a) η2 = 0.2, (b) max. μ – min. μ / within-group σ = 2, (c) maximum % difference from the control mean of 25%. (see Methods for detailed descriptions of effect sizes) Hypothesis 2a/b/c. Glyphosate reduces species richness with effect sizes of (a) η2 = 0.2, (b) max. μ – min. μ / within-group σ = 2, (c) maximum % difference from the control mean of 25%. Hypothesis 3 a/b/c. Glyphosate reduces and diversity with effect sizes of (a) η2 = 0.2, (b) max. μ – min. μ / within-group σ = 2, (c) maximum % difference from the control mean of 25%. Hypothesis 4. The effects of glyphosate application are mitigated by nutrient inputs.

1. A low target aqueous glyphosate concentration (210 μg a.e./L) (L). This is the upper 99th percentile concentration of glyphosate observed in wetlands located in agricultural areas in Southern Ontario (Byer et al., 2008; Struger et al., 2008), but still well below concentrations expected from direct overspray. 2. A high target aqueous glyphosate concentration (2, 880 μg a.e./L) (H). This is the direct overspray concentration expected in a 15 cm deep wetland with no intercepting vegetation that has been oversprayed at the maximum label rate (7.9 L/ha). 3. The low target aqueous glyphosate concentration together with nutrient additions (LN). 4. The high target aqueous glyphosate concentration together with nutrient additions (HN).

Hypothesis 5. The effects pf glyphosate application are persistent. 2. Methods 2.1. Study location The field site was located in the Long Term Experimental Wetlands Area (66°29′59.02″W, 45°40′48.62″N) on Canadian Forces Base Gagetown (CFB Gagetown) about 60 km from the city of Saint John in southern New Brunswick, Canada (Appendix A). CFB Gagetown is a 1200 km2 military installation primarily used for ground training. The CFB Gagetown site was selected for this large-scale, multi-year experiment because there is a large landbase, a single decision-making authority, a low probability that land ownership will change hands, and a desire to operate the area in a more environmentally sustainable way. It is a mixed-wood Acadian forest, dominated by red spruce (Picea rubens), white spruce (Picea glauca), red maple (Acer rubrum), white birch (Betula papyrifera) and yellow birch (Betula alleghaniensis). The six km2 study area was mechanically cleared of trees and the top soil pushed into berms in 1997 and 1998. Wetlands formed naturally in depressions and next to the windrows because of the disruption of drainage patterns of the area. This created a landscape containing hundreds of small wetlands ranging from < 1 ha to several hectares in size, which have

The two target concentrations represent the maximum environmentally observed concentrations in agricultural wetlands (L) and a reasonable worst case scenario assuming complete overspray of the entire wetland at the maximum label rate (H). The two treatments that combine herbicide and nutrients were done to simulate an agricultural multiple stressor scenario. The nutrient treatments were designed to increase the trophic status of each wetland one level (e.g. oligotrophic to mesotrophic) above background levels (calculated in 2008 – one year prior to chemical treatment). Each treatment (i.e. L, H, LN, HN) was replicated in the treatment half of 6 wetlands. 2

Ecotoxicology and Environmental Safety 183 (2019) 109475

J.F. Mudge and J.E. Houlahan S

The herbicide was applied to the wetlands twice in the spring of 2009 and twice in the spring of 2010. In 2009 the first herbicide application occurred on 15/16 May and the second application on 9/10 June. In 2010 the first application occurred on 26/27 April and the second application occurred on 24/25 May. We used two applications to simulate standard agricultural practices for glyphosate tolerant crops and timed to correspond with spring macrophyte emergence in each year. Post-application sampling of ponds showed that glyphosate concentrations were statistically significantly higher in the high concentration halves than the low concentration halves although we weren't able to precisely hit target concentrations. There was no evidence of contamination on the control sides of ponds (Edge et al., 2012). The 12 nutrient treatment wetlands were sprayed on the treatment half with nutrients commonly used in fertilizers (technical grade ammonium nitrate and phosphoric acid, purchased from Fisher Scientific). In 2008, wetlands were generally classified as oligotrophic/mesotrophic based on aqueous nutrient concentrations (Wetzel, 2001), and nutrients were added to increase the aqueous total Kjeldahl nitrogen (TKN) concentrations by one trophic level, and phosphorus was added to maintain the TKN:TP ratios observed in 2008. Nutrients were added to wetlands four times in each of 2009 and 2010. Timing of nutrient additions was determined by monitoring weekly wetland water chemistry data and adding nutrients to maintain TKN and TP at their desired levels. In 2009 nutrients were added to wetlands on 14 May, 29 May, 3 July and 19 Aug. In 2010 nutrients were added to wetlands on 22/23 April, 13/14 May, 24/25 June, and 30/31 July. Immediately before each application, total wetland volume was calculated from depth measurements taken at regularly spaced intervals along transects of each wetland. We then used the estimated wetland water volume to estimate the volume of herbicide/nutrient required to achieve a particular target concentration. The herbicide/nutrients were then mixed in ~3 L of wetland water and applied evenly to the surface of the wetland with a backpack sprayer. The same volume of water was applied to the control side using a separate backpack sprayer to ensure equivalent physical disturbance on both sides of each wetland. We estimated species composition of both emergent and submerged macrophyte species of each wetland half (Crow and Hellquist, 2000a, 2000b) and percent cover to the nearest 1% (Anderson, 1986) for each species maintaining the same search effort per unit area for each wetland half (approximately 3 s/m2 on average). Plant surveys were conducted in late May, between the first and second herbicide treatments, and twice more after the second herbicide treatment, once in late June and once in August in 2009 and 2010. To assess recovery from herbicide and nutrient applications we conducted plant surveys in late June, 2011, August 2011 and early August of 2012.

Evenness = ∑i = 1 pi *ln(pi ) ln(S ) ,

1

where S is the total number of species in the community and p is the proportion of the total cover that is comprised of the cover of species i. The Bray-Curtis similarity between wetland halves was calculated as S

2* ∑i = 1 minCoverTCi

(TotalCoverT + TotalCoverC ) ,

2

where S is all species in common between the treatment and control halves of a pond, minCoverTCi is the cover of species i on the side where it was lowest of the treatment and control sides, TotalCoverT is the total macrophyte cover on the treatment side and TotalCoverC is the total macrophyte cover on the control side. Species evenness and community similarity range between 0 and 1 so were arcsine transformed prior to statistical analysis. We used partially nested repeated measures linear mixed effect Type III ANOVA models to examine the influence of wetland half (treated vs. control), in combinations with herbicide concentration, nutrient concentration, and/or year, on species richness, Pielou's species evenness, and max. total cover. Herbicide concentration and nutrient treatments were each considered between subjects fixed factors, with pond as an incomplete random blocking factor, side (treated or control) as a within-subjects fixed factor and year as the within subjects repeated measure random factor. For the Bray-Curtis similarity endpoint, side was not possible to include as a factor because it was incorporated into the similarity calculation. Thus, we used an analogous ANOVA model to examine the influences of different combinations of herbicide concentration (between subjects fixed factor), nutrient treatment (between subjects fixed factor) and time (within subjects random repeated measure factor) on Bray-Curtis similarity. Optimal α: Statistical significance of ANOVA terms was evaluated using the optimal α approach for minimizing the combined probabilities or costs of Type I and II errors in null Hypothesis tests (Mudge et al., 2012). The calculation of optimal α levels for ANOVA terms requires estimates of the relative costs of Type I and II errors, the relative prior probabilities of null and alternate hypotheses, the numerator and denominator degrees of freedom for the ANOVA term, and an estimate of the critical (i.e. biologically relevant) effect size. As there was no strong rationale for considering one type of error to be more serious than the other or to assume different prior probabilities of the null and alternate hypotheses, we set the relative costs of Type I and II errors and prior probabilities to be equal. 2.4. Critical effect size Critical effect sizes can be difficult to conceptualize in ANOVA designs because researchers are often interested in pairwise differences among groups while effects in ANOVA designs usually incorporate variability across all groups, rather than pairwise differencesWe used (1) the η2 approach, (2) the max. – min./within-group σ approach and (3) the max. % difference from control mean approach. All 3 effects sizes must be converted to Cohen's f2 to be implemented (Cohen, 1988). The η2 critical effect size measure (Cohen, 1988) is analogous to R2 in that it represents the proportion of the total variance in the ANOVA model that is explained by group membership (either the factor levels or all the combinations of multiple factor levels for interaction terms). We chose η2 = 0.2 as the minimum biologically relevant proportion of total variance explained by group membership. f2 = η2/1- η2 therefore f2 = 0.25 for each endpoint and each ANOVA term. The maximum group mean – minimum group mean/within-group σ critical effect size is analogous to the Cohen's d critical effect size used as an effect size measure in t-tests (the difference between group means relative to the within-group standard deviation) in that it represents the smallest possible among-group variance that can produce a maximum difference between any two group means at least as large as some

2.3. Data analysis This study focuses on four plant community endpoints; i) species richness (Appendix B), ii) Pielou's species evenness index (Magurran, 2004) (Appendix C), iii) the sum of the peak cover for each species over the summer (Appendix D) and iv) the Bray-Curtis similarity index (Magurran, 2004) between the treated and control side of each wetland. We calculated community endpoints for each wetland half/year combination using two different techniques, maximum and mean cover, for pooling across the three samples. The sample containing the maximum cover of each species was used to estimatespecies richness and total cover for each wetland half in each year. The average cover of each species across all three samples was used to estimate evenness and community similarity. Thus, species richness for a wetland half was the count of all species that were detected in at least one of the surveys in any year. Total cover was calculated by summing the maxima for each species. When the largest observed abundance of a species across all three surveys was < 1%, the maximum abundance for that species was assumed to be 0.25% (the midpoint between 0% and 0.5%, which would round up to 1%). The Pielou's species evenness was calculated as 3

Ecotoxicology and Environmental Safety 183 (2019) 109475

J.F. Mudge and J.E. Houlahan

specified multiple of the within-group standard deviation (f2 = max. μ – min. μ/within-group σ*2 * k, where k is the number of groups). We chose the minimum biologically relevant difference between any two groups to be where max. μ – min. μ/within-group σ equals 2. This critical effect size value was chosen to be consistent with an effect size of 2 standard deviations proposed by Munkittrick et al. (2009) for a wide variety of biological and ecological endpoints. Because this critical effect size incorporates the number of groups the critical effect size will be smaller for factors with more groups, all other things being equal. The max. % difference from the control mean critical effect size represents the smallest possible among-group variance that can produce a maximum difference between any two group means at least as large as some specified percent difference from the control mean. The first two effect sizes were variance–standardized but max. % difference is not, and thus requires an estimate of both the expected control mean and the expected within-group variance. The estimates of the control mean and the within-group variance were taken from the observed data. For each ANOVA term, the max. % difference from the control mean critical effect size was converted to an f2 critical effect size by: (1) multiplying the control mean by 1 ± the critical maximum % difference from the control mean/100 to get the minimum and maximum of group means, (e.g. control mean = 5; critical maximum % difference = 25; 5*1 ± 25/100 = 3.75 and 6.25) (2) calculating the minimum amonggroup variance for the k group means by assuming all other group means are clustered at the midpoint between the minimum and maximum, then (3) dividing the minimum among-group variance (from Step 2) by the within-group variance from the data. We chose the minimum biologically relevant difference between any two groups to be a 25% difference from the control mean. This critical effect size value was also recommended by Munkittrick et al. (2009) for a variety of environmental monitoring endpoints. We calculated optimal α and associated β levels for each of the three potential critical effect sizes using an iterative process described in Mudge et al. (2012) (Table 1). The different critical effect sizes and degrees of freedom for each endpoint and ANOVA term combination produce different optimal α levels. Tests associated with larger optimal α levels imply weaker experimental design and thus, weaker evidence for effects (or lack thereof) than results associated with smaller optimal α levels. We classified the different optimal α levels into three categories corresponding to different levels of confidence in the statistical outcome. For optimal α ≤ 0.05, we considered any p-values ≤ α to constitute strong evidence for an effect and any p-values > α to be strong evidence that the effect was smaller than the critical effect size. For optimal α between 0.05 and 0.15, we considered any p-values ≤ α to be moderate evidence for an effect and p-values > α, to be moderate evidence that the effect was smaller than the critical effect size. Finally, for optimal alpha > 0.15 we considered any p-value ≤ α to be weak evidence for an effect, and a p-value > α to be weak evidence that the effect was smaller than the critical effect size.

f2 = 0.25 (Table 1, Fig. 1). The median f2 converted critical effect size value for the max. diff. ≥ 2 σ(within-group) measure was also 0.25, with larger converted f2 values for the ANOVA terms containing few groups (side, side*concentration and side*nutrients) and smaller converted f2 values for ANOVA terms containing many groups (side * year * concentration, side * year * nutrients, and side * year * concentration * nutrients). This is because the among group variance required to produce a constant maximum difference between any 2 groups decreases as the number of groups increases (additional group means added to the midpoint between the minimum and maximum would decrease the variance, holding the difference between the minimum and maximum groups constant). The max. diff. ≥ 25%μ(control) critical effect size measure translated to f2 values with a lower median value of 0.16. This effect size measure also had the largest range of f2 critical effect sizes (0.02–2.7) because the f2 critical effect size was dependent not only on the number of groups contained within each ANOVA term, but also on the size of a 25% difference from the control mean relative to the within-group variance. Optimal α levels ranged from 0.00010 to 0.35 among the ANOVA terms for the three potential critical effect sizes, with a mean α level of 0.13, a mean β level of 0.15 and a mean average of α + β of 0.14 (Table 1). Where we use α > 0.05 we increase the risk of committing a Type I error over the standard practice of using α = 0.05, but this is more than offset by the reduced risk of committing a Type II error. 3. Results 3.1. Species richness A total of 98 different taxa were found during plant surveys (Appendix E).The effect of wetland side on species richness differed among years for each of the three potential critical effect sizes (Table 2, Fig. 2). Richness was lower on treated sides than on control sides in 2009 and 2010, but in the 2011 and 2012 years when the treated sides did not receive further treatment, the sides had similar richness. For the η2 ≥ 0.2 and max. diff ≥ 2σ(within group) critical effect sizes, the α levels for the side*year interaction were each 0.018, and the observation of pvalues smaller than these low α levels indicated strong evidence for a real effect. By contrast, there was only weak evidence of a side*year interaction for the max. diff ≥25%μ(control) critical effect size, with p < α = 0.19. The consistently significant interaction between side and year for all potential critical effect sizes precluded interpretation of the effect of side independent of the effect of year. 3.2. Species evenness The effect of wetland side on species evenness also differed among years, regardless of which of the three potential critical effect sizes was used (Table 3, Fig. 3). The mean difference in species evenness between treated and control sides of wetlands increased from 0.038 to 0.049 from 2009 to 2011 (the year following the two years in which the treatments were applied), but in 2012, the mean difference dropped to 0.014. This observed effect constitutes strong evidence of a real side*year effect for species evenness using each of the three potential critical effect sizes, as the p-value was less than the already low optimal α levels of 0.018, 0.018 and 0.00010 for the η2 ≥ 0.2, max. diff ≥ 2σ(within group) and max. diff ≥25%μ(control) critical effect sizes, respectively. The significant interaction between side and year for all potential critical effect sizes precluded interpretation of the effect of side independent of the effect of year.

2.5. Critical effect sizes and optimal α levels After translating each of the three critical effect size measure into f2 values, each ANOVA term and endpoint combination had at least two different potential critical effect size values – occasionally two different critical effect size measures resulted in the same critical effect size value (22 of 31 ANOVA term and endpoint combinations had three different critical effect size values) (Table 1). There were eight ANOVA term and endpoint combinations for which the η2 ≥ 0.2 and the max. diff. ≥ 2σ(within-group) critical effect sizes corresponded to the same f2 value and one case where the η2 ≥ 0.2 and max. diff. ≥ 25%μ(control) critical effect sizes corresponded to nearly the same f2 value. The η2 ≥ 0.2 critical effect size (corresponding to at least 20% of the total variance in the data being explained by the factor level membership) translated to an among-group/within-group variance ratio of

3.3. Macrophyte cover The factors influencing macrophyte cover depend on what critical 4

Ecotoxicology and Environmental Safety 183 (2019) 109475

J.F. Mudge and J.E. Houlahan

Table 1 Critical effect sizes expressed as a ratio of among-group variance to within-group variance (f2), optimal α and associated β levels for each endpoint and ANOVA term combination under three different approaches to selecting a critical effect size. (values in bold, regular and italic font are strong, moderate and weak evidence respectively). endpoint

ANOVA term

richness

side side*concentration side*nutrients side*concentration*nutrientsside*year

evenness

side*year*concentration side*year*nutrients side*year*concentration*nutrients Side side*concentration side*nutrients side*concentration*nutrientsside*year

cover

side*year*concentration side*year*nutrients side*year*concentration*nutrients side side*concentration side*nutrients side*concentration*nutrientsside*year

similarity

side*year*concentration side*year*nutrients side*year*concentration*nutrients concentration nutrients concentration*nutrients year concentration*year nutrients*year concentration*nutrients*year

f2

0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25

η2 ≥ 0.2

max diff. ≥ 2σ(within

max diff. ≥ ± 25%μ(control)

group)

α

β

f

α

β

f2

α

β

0.18 0.22 0.22 0.220.018

0.22 0.25 0.25 0.25 0.019 0.019 0.019 0.019 0.22 0.25 0.25 0.25 0.019 0.019 0.019 0.019 0.22 0.25 0.25 0.25 0.019 0.019 0.019 0.019 0.19 0.19 0.19 0.068 0.068 0.068 0.068

1.0 0.50 0.50 0.25 0.25 0.13 0.13 0.063 1.0 0.50 0.50 0.25 0.25 0.13 0.13 0.063 1.0 0.50 0.50 0.25 0.25 0.13 0.13 0.063 1.0 1.0 0.50 0.50 0.25 0.25 0.13

0.033 0.13 0.13 0.22 0.018 0.088 0.088 0.20 0.033 0.13 0.13 0.22 0.018 0.088 0.088 0.20 0.033 0.13 0.13 0.22 0.018 0.088 0.088 0.20 0.022 0.022 0.077 0.011 0.063 0.063 0.16

0.028 0.13 0.13 0.25 0.019 0.10 0.10 0.24 0.028 0.13 0.13 0.25 0.019 0.10 0.10 0.24 0.028 0.13 0.13 0.25 0.019 0.10 0.10 0.24 0.019 0.019 0.078 0.011 0.068 0.068 0.19

0.26 0.13 0.13 0.064 0.064 0.032 0.032 0.016 2.7 1.3 1.3 0.67 0.67 0.34 0.34 0.17 0.31 0.16 0.16 0.078 0.078 0.039 0.039 0.019 0.64 0.64 0.32 0.32 0.16 0.16 0.079

0.18 0.29 0.29 0.34 0.19 0.29 0.29 0.35 0.0021 0.026 0.026 0.088 0.00010 0.0058 0.0058 0.050 0.15 0.27 0.27 0.33 0.16 0.26 0.26 0.34 0.054 0.051 0.13 0.039 0.12 0.12 0.22

0.21 0.37 0.37 0.47 0.23 0.36 0.36 0.45 0.0013 0.020 0.020 0.083 0.000085 0.0061 0.0061 0.057 0.18 0.34 0.34 0.45 0.19 0.33 0.33 0.43 0.051 0.051 0.15 0.041 0.14 0.14 0.28

0.018 0.018 0.018 0.18 0.22 0.22 0.220.018 0.018 0.018 0.018 0.18 0.22 0.22 0.220.018 0.018 0.018 0.018 0.16 0.16 0.16 0.063 0.063 0.063 0.063

2

Fig. 1. Critical effect sizes (standardized as f2 values, the ratio of the size of the among-group variance to the within-group variance), optimal α and associated β levels summarized over each ANOVA term and endpoint (richness, evenness, cover and similarity) combination, for each of the three critical effect size measures.

effect size is used. For the max. diff ≥ 2σ(within group) and max. diff ≥25%μ(control) critical effect sizes, there was weak evidence of a side*year*concentration*nutrients effect (Table 4). The significant side * year * concentration * nutrient interaction represented the larger mean differences between treated and control sides of H wetlands relative to other wetland treatments (HN, L, and LN) in 2010 (Fig. 4). For the larger, η2 ≥ 0.2 critical effect size, there was strong evidence suggesting that the side*year*concentration*nutrient interaction effect

does not explain at least 20% of the variability in macrophyte cover, p > α and β = 0.019. There was, however, strong evidence for a side*year interaction effect under the η2 ≥ 0.2 critical effect size (p < α = 0.018) whereby the macrophyte cover was lower on treated sides than control sides during years receiving treatments (2009 and 2010) but not in the following years (2011 and 2012), and weak evidence for a side*concentration*nutrient interaction effect under the η2 ≥ 0.2 critical 5

Ecotoxicology and Environmental Safety 183 (2019) 109475

J.F. Mudge and J.E. Houlahan

Table 2 Degrees of freedom, observed effect sizes and p-values of the ANOVA for the difference in species richness between wetland sides. P-values not in bold were not significant using optimal α with any of the three potential critical effect sizes.

Table 3 Degrees of freedom, observed effect sizes and p-values of the ANOVA for the difference in Pielou's species evenness between wetland sides. P-values not in bold were not significant using optimal α with any of the three potential critical effect sizes.

df

df

intercept side side*concentration side*nutrients side*concentration*nutrients side*year side*year*concentration side*year*nutrients side*year*concentration*nutrients

1, 1, 2, 2, 2, 6, 6, 6, 6,

120 17 17 17 17 120 120 120 120

F-value

p-value

691.001 4.912 1.019 0.053 0.317 55.069 0.909 0.887 0.584

< 0.0001 0.0406*‡ 0.382 0.9485 0.7325 < 0.0001*†‡ 0.4912 0.507 0.7425

intercept side side*concentration side*nutrients side*concentration*nutrients side*year side*year*concentration side*year*nutrients side*year*concentration*nutrients

1, 1, 2, 2, 2, 6, 6, 6, 6,

120 17 17 17 17 120 120 120 120

F-value

p-value

3416.889 18.05 0.535 1.046 0.464 5.862 0.349 0.648 1.096

< 0.0001 0.0005*†‡ 0.5952 0.3728 0.6366 < 0.0001*†‡ 0.9093 0.692 0.369

*p-value ≤ optimal α for the η2 ≥ 0.2 effect size. †p-value ≤ optimal α for the max diff. ≥ 2σ(within group) effect size. ‡p-value ≤ optimal α for the max diff. ≥ ± 25%μ(control) effect size.

*p-value ≤ optimal α for the η2 ≥ 0.2 effect size. †p-value ≤ optimal α for the max diff. ≥ 2σ(within group) effect size. ‡p-value ≤ optimal α for the max diff. ≥ ± 25%μ(control) effect size.

effect size (p < α = 0.22), with H wetlands having slightly less macrophyte cover on treated sides relative to their controls compared to all other treatment categories (HN, L, and LN). The significant side*year*concentration*nutrient interaction for the max. diff ≥ 2σ(within group) and max. diff ≥25%μ(control) critical effect sizes precluded interpretations of main effects, 2-way and 3-way interactions, while the significant side*year and side*concentration*nutrient interactions for the η2 ≥ 0.2 critical effect size precluded interpretation of the main effects.

similarity between wetland sides in 2010, relative to 2009, 2011 and 2012 (Table 5, Fig. 5). 4. Discussion 4.1. Effects of glyphosate on wetland plant communities The application of glyphosate-based herbicides was found to result in a decrease in macrophyte species richness, an increase in macrophyte species evenness, a decrease in macrophyte cover and a reduction in community similarity. These effects were evident in the first year of herbicide application in 2009, and became more pronounced in the second year of herbicide application in 2010. However, when herbicides were not applied in 2011, recovery was observed in most endpoints, with the exception being species evenness, for which partial recovery was not observed until 2012.

3.4. Similarity in community composition Averaged across years, there was weak evidence for a difference in Bray-Curtis similarity between the H and L treatment concentrations under the η2 ≥ 0.2 critical effect size (p < α = 0.16) (Table 5). Wetlands receiving the high glyphosate treatment had a mean similarity between control and treated sides of 62.0%, while wetlands receiving the low glyphosate concentration had a mean similarity between control and treated sides of 69.0%. There was also moderate evidence (critical effect size η2 ≥ 0.2, p < α = 0.063) to strong evidence (critical effect sizes max. diff ≥ 2σ(within group) and max. diff ≥25%μ(control), p < α = 0.011 and p < α = 0.039 respectively) for a significant difference in Bray-Curtis similarity among years. This difference in similarity among years is primarily due to a lower observed Bray-Curtis

4.2. Richness In contrast to the lower species richness we observed on treated sides of wetlands relative to their control sides during years where the herbicide was applied, other studies have reported no difference in species richness, or sometimes higher species richness following glyphosate herbicide application. Wetland restoration experiments that

Fig. 2. Mean differences in species richness ( ± 1 SE) between treated and control sides of wetlands from 2009 to 2012. 6

Ecotoxicology and Environmental Safety 183 (2019) 109475

J.F. Mudge and J.E. Houlahan

Fig. 3. Mean differences in Pielou's species evenness ( ± 1 SE) between treated and control sides of wetlands from 2009 to 2012. Table 4 Degrees of freedom, observed effect sizes and p-values of the ANOVA for the difference in the sum of maximum macrophyte cover for each species between wetland sides. P-values not in bold were not significant using optimal α with any of the three potential critical effect sizes. df intercept side side*concentration side*nutrients side*concentration*nutrients side*year side*year*concentration side*year*nutrients side*year*concentration*nutrients

1, 1, 2, 2, 2, 6, 6, 6, 6,

120 17 17 17 17 120 120 120 120

F-value

p-value

1132.551 46.3185 0.7718 0.3446 2.0251 14.2379 0.6355 1.8215 2.2881

< 0.0001 < 0.0001*†‡ 0.4777 0.7133 0.1626*†‡ < 0.0001*†‡ 0.7016 0.1004‡ 0.0398†‡

Table 5 Degrees of freedom, observed effect sizes and p-values of the ANOVA for BrayCurtis similarity in community composition between wetland sides. P-values not in bold were not significant using optimal α with any of the three potential critical effect sizes. df intercept concentration nutrients concentration*nutrients Year concentration*year nutrients*year concentration*nutrients*year

1, 1, 1, 1, 3, 3, 3, 3,

60 20 20 20 60 60 60 60

F-value

p-value

2658.3748 4.0712 0.2241 0.5823 6.501 0.1958 1.2276 0.2778

< 0.0001 0.0572* 0.6411 0.4543 0.0005*†‡ 0.8988 0.3076 0.8412

*p-value ≤ optimal α for the η2 ≥ 0.2 effect size. †p-value ≤ optimal α for the max diff. ≥ 2σ(within group) effect size. ‡p-value ≤ optimal α for the max diff. ≥ ± 25%μ(control) effect size.

*p-value ≤ optimal α for the η2 ≥ 0.2 effect size. †p-value ≤ optimal α for the max diff. ≥ 2σ(within group) effect size. ‡p-value ≤ optimal α for the max diff. ≥ ± 25%μ(control) effect size.

Burge et al., 2017). In terrestrial agricultural systems, glyphosate herbicide use for weed control in glyphosate-tolerant soybean fields can, under some circumstances, decrease weed species richness early in the growth season but tended to increase weed species richness in the same

have applied glyphosate herbicides to control invasive plant populations have reported increases in non-target plant species richness following herbicide use (Paveglio and Kilbride, 2000; Ailstock et al., 2001;

Fig. 4. Mean differences in the sum of the maximum percent cover for each species ( ± 1 SE) between treated and control sides of wetlands from 2009 to 2012. The double line with diamond markers represents wetlands treated with the high concentration of the glyphosate herbicide on one side, the solid line with triangle markers represents wetlands treated with the high concentration of the glyphosate herbicide plus nutrients on one side, the dotted line with diamond markers represents wetlands treated with the low concentration of the glyphosate herbicide and the dashed line with square markers represents wetlands treated with the low concentration of the glyphosate herbicide plus nutrients on one side.

7

Ecotoxicology and Environmental Safety 183 (2019) 109475

J.F. Mudge and J.E. Houlahan

Fig. 5. Mean Bray-Curtis similarity ( ± 1 SE) between treated and control sides of wetlands in 2009–2012 (solid line). For comparison, the upper dashed line represents the mean control side similarity between two consecutive years, and the lower dashed line represents the mean control side similarity between wetlands for each year.

fields later in the growth season (Vitta et al., 2004). A review of the effects of glyphosate use on terrestrial plant communities by Sullivan and Sullivan (2003) shows overlapping 95% CI of species richness between control and glyphosate treated plots (primarily from forestry scenarios). This study may have had a greater capability of observing reductions in species richness following glyphosate herbicide application due to the comparatively large level of replication relative to other studies (24 replicate wetlands), the elimination of among-wetland variability as a confounding factor, through use of a split-plot design, and the use of the optimal α approach which ensures optimal statistical power to detect biologically relevant effects. The reductions in species richness were small (i.e. 3–5 species) and only present in the years in which herbicide was applied, making them potentially easy to miss.

application provided residual habitat structure which may reduce potential indirect effects of glyphosate herbicide application on invertebrates and/or amphibians. The lack of continuing differences in macrophyte cover between treated and control sides of wetlands in nonapplication years (2011 and 2012) may also decrease the likelihood of indirect impacts of glyphosate herbicides in wetlands at higher trophic levels, if runoff or accidental overspray events do not occur every year. 4.5. Community composition Few studies have examined the influence of glyphosate herbicide application on wetland plant community composition, aside from studies focused on the effectiveness of glyphosate herbicides in reducing dominance by an invasive species (Jadhav et al., 2008; Mozdzer et al., 2008; Linz and Homan, 2011). Dalton and Boutin (2010) showed variable susceptibility to a glyphosate herbicide among seven nontarget wetland macrophyte species in a laboratory setting and community composition changes when these species were grown together in greenhouse microcosms. In this study, only minor changes in community composition were observed, and only after two years of herbicide application, with similarity rebounding in non-application years. However, we did not examine changes at the individual species level and we will be assessing species-level effects and recovery in a future paper.

4.3. Evenness Species evenness has not commonly been examined as a plant community endpoint in response to glyphosate herbicide application, however several studies have reported increases in plant species diversity following herbicide application (Paveglio and Kilbride, 2000; Ailstock et al., 2001), with some attributing the diversity increases to reductions in the abundances of previously dominant species (Freedman et al., 1993) or increases in the abundances of previously less common species (Horsley, 1994; Mozdzer et al., 2008). Although the size of the difference in species evenness between treated and control halves of wetlands had begun to decrease by 2012, the persistence of higher species evenness on treated sides relative to control sides in non-spray years (2011 and 2012) suggests that glyphosate herbicide application may be a useful component of wetland restoration management and control of invasive wetland plant species.

4.6. Concentration and nutrient dependent effects There was very little evidence of different effects of herbicide application on macrophyte communities between wetlands treated with the environmentally observed concentration of glyphosate and wetlands treated with the predicted maximum environmental concentration of glyphosate, except for the weak evidence that there may be lower community similarity between treated and control sides of wetlands treated with the higher glyphosate concentration. This lack of different response strengths at different glyphosate concentrations is in contrast to most studies of effects of glyphosate herbicides on wetland plants, which generally describe a concentration dependent response (Jadhav et al., 2008; Dalton and Boutin, 2010; Linz and Homan, 2011; Kielak et al., 2011). The lack of a concentration dependent effect in this study likely results from the additional spray at the maximum recommended label rate targeted specifically at emergent macrophytes that all wetlands received immediately following the application of the target concentration to the water. The purpose of the additional glyphosate application directly targeting the macrophyte community was

4.4. Cover The observation of a significant reduction in plant cover on treated sides of wetlands relative to their control sides following glyphosate herbicide application is not surprising, as reducing plant cover is the intended purpose of glyphosate herbicides. Studies are frequently designed to determine the most efficient herbicide concentration for reducing the cover of invasive wetland plants (Jadhav et al., 2008; Mozdzer et al., 2008; Linz and Homan, 2011). More surprising was that cover differences between treated and control sides of wetlands, although present, were not more extreme. The persistence of > 50% macrophyte cover on treated sides of wetlands post herbicide 8

Ecotoxicology and Environmental Safety 183 (2019) 109475

J.F. Mudge and J.E. Houlahan

as a % difference from the control mean. Also, while this effect size measure will remain consistent among different endpoints using the same ANOVA design, implementation in statistical analysis is complicated by the need to translate the effect size into a unique f2 critical effect size for each ANOVA term (for minimum among-group variance, f2 = d2/2 k, where k is the number of groups). Seven of the 8 ANOVA term and endpoint combinations that would have been significant using α = 0.05 were also significant using a maximum difference between any two groups as large as 2 standard deviations of the data as the critical effect size, with one additional significant ANOVA term and endpoint combination that would not have been significant using α = 0.05. The η2 critical effect size measure is relatively easy to implement into statistical analysis because there is a simple formula to translate this effect size into an f2 critical effect size (f2 = η2/(η2-1)), and the critical effect size remains the same among endpoints and among ANOVA terms. This critical effect size measure can be relatively easy to conceptualize because it is analogous to the R2 estimate of regression models in that the effect is expressed as the ratio of the variability explained by group membership relative to the total variability in the dataset. Specification of a threshold η2 value that represents biological relevance can be problematic, however, because in many study designs the pairwise differences between groups are of greater importance than the variability among all groups. Seven of the 8 ANOVA term and endpoint combinations that would have been significant using α = 0.05 were also significant using η2 = 0.2 as the critical effect size, with two additional significant ANOVA term and endpoint combinations that would not have been significant using α = 0.05.

to maximize the possibility of indirect impacts of glyphosate herbicides on the invertebrate or amphibian communities through direct effects to the plant community. However, this consistent amount of herbicide applied directly to the plant community on the treated sides of all wetlands was also much higher than the dose received through the different treatment concentrations applied directly to the water's surface. There was also very little evidence of different effects of herbicide application on macrophyte communities between nutrient and no-nutrient treatments. The nutrients added directly to the water column may have been rapidly sequestered by algae and may not have been available to macrophytes. 4.7. Critical effect sizes, optimal α and Hypothesis testing with ANOVA designs 4.7.1. Critical effect size measures in ANOVA The results for 27 of the 31 ANOVA term and endpoint combinations were consistent among the three different critical effect sizes. This suggests that the optimal α approach can be robust to moderate variation in choice of critical effect size among researchers. When obvious and large effects exist, they will be detected using the optimal α approach even when there is disagreement about the most appropriate minimum threshold for biological relevance. This is because, if there are large effects then, by definition, there are medium and small effects. Thus, if you set the threshold for a large effect size and obtain a statistically significant result you will also get a statistically significant result for thresholds targeting small or medium effects. The four instances where the statistical outcome was not consistent among the three critical effect sizes can be interpreted as having weak evidence for an effect if smaller effect sizes are considered biologically relevant but strong evidence that there is not a large effect if only large effects are biologically relevant. There is no clearly superior critical effect size measure for ANOVA designs. Among the three critical effect size measures we examined there is a tradeoff between ease of conceptualization for determination of what constitutes biological relevance and ease of implementation in statistical analysis. The max. % difference from the control mean critical effect size measure is easy to conceptualize and relate to biological relevance, but in order to incorporate this critical effect size measure into the calculation of an optimal α level, it requires estimates of the control mean and the within-group standard deviation. The biologically relevant percent difference from the control mean must be translated to an absolute difference between the minimum and maximum group means and then divided by the within-group standard deviation to produce a Cohen's d value which must then be converted to an f2 value (for minimum among-group variance, f2 = d2/2 k, where k is the number of groups). The difficulty in implementing this critical effect size in statistical analysis is further compounded by the need to translate it into a unique f2 critical effect size for each endpoint and ANOVA term combination (because the Cohen's d value is calculated using the control mean, which differs among endpoints, and because the number of groups differs among ANOVA terms). All of the 8 ANOVA term and endpoint combinations that would have been significant using α = 0.05 were also significant using a maximum difference between any two groups as large as a 25% difference from the control mean as the critical effect size, as well as two additional significant ANOVA term and endpoint combinations that would not have been significant using α = 0.05. The max. – min./within-group σ critical effect size (Cohen's d) is still somewhat easy to relate to biological relevance in that it focuses on maximum differences between any two groups, however the specification of a biologically relevant difference between any two groups in units relative to the within-group standard deviation may be more conceptually difficult than expressing a biologically relevant difference

4.7.2. Hypothesis testing using ANOVA designs Critical effect sizes and statistical power are frequently ignored in ANOVA designs. This is because ANOVA study designs are frequently complex and often used in situations where the pairwise differences between groups are of greater relevance than the overall variability among groups. More explicit consideration of critical effect sizes can, however, improve statistical results in ANOVA designs. Incorporating critical effect sizes into the calculation of optimal α levels can both decrease the chances of statistical errors when sample sizes are large and when only large effects are biologically relevant, while at the same time increasing the likelihood of detecting biologically meaningful effects in situations that would typically have low statistical power to detect effects using α = 0.05. While choosing appropriate measures and magnitudes of critical effect sizes and incorporating them into the statistical analyses may not always be easy, especially for complex ANOVA designs, it can improve research quality through increases in strength of statistical inference and ease of interpretation by associating statistical significance with biological significance. 4.8. Conclusions (1) Glyphosate application reduced species richness and % cover by at least 25% relative to the control side. (2) Glyphosate application increased evenness by, at least, 25% relative to the control sides. (3) All impacts of glyphosate applications in 2009 and 2010 were dramatically reduced or absent by 2012. (4) There was no evidence of concentration-dependent effects or interactive effects of nutrient applications. (5) In most tests, the targeted effect size did not affect statistical significance but occasionally evidence for small to medium effects was found but there was not evidence for large effects. Acknowledgements Thanks to the Natural Sciences and Engineering Research Council of Canada for funding of this project through both their Discovery and 9

Ecotoxicology and Environmental Safety 183 (2019) 109475

J.F. Mudge and J.E. Houlahan

Strategic Grant programs.

Appendix A. Map of the field sites on CFB Gagetown in New Brunswick, Canada

Appendix B. Species richness of macrophytes in treatment and control halves in each year

Treatment

Wetland

2009

H

2 13 22 24 37 43 7 12 19 21 36 41

42

HN

46

2010 44 51 28 46 40 45 65 37 44 47 52 28

41

41

41 44 36 45 38 43 49 36 47 45 42 28

48

50

2011 45 58 37 50 41 58 70 39 45 53 50 41

43

43

46 44 37 45 37 48 48 32 45 53 49 32

37

42

10

2012 35 40 30 41 33 41 53 39 35 46 45 31

36

39

42 33 38 38 23 40 40 38 33 45 45 30

32

31

All years 32 35 29 30 30 33 48 21 27 34 31 26

29

31

31 32 33 30 19 28 42 23 27 37 31 23

40

42

39 46 31 42 36 44 59 34 38 45 45 32

37

38

40 38 36 40 29 40 45 32 38 45 42 28

Ecotoxicology and Environmental Safety 183 (2019) 109475

J.F. Mudge and J.E. Houlahan L

LN

3 5 10 23 39 42 11 15 20 26 32 40

All wetlands

43

45

41 49 38 38 50 42 48 59 44 53 38 26

44

38

45

40 42 36 35 43 31 53 58 38 56 45 19

41

45

47

48 61 35 41 45 42 52 56 51 46 47 29

48

41

42

48 54 29 34 42 36 45 53 37 46 44 26

42

35

35

43 47 30 22 36 32 31 40 36 40 36 24

37

36

39

39 46 35 25 39 31 37 44 41 43 42 28

37

25

24 42 17 21 26 22 28 38 37 26 33 26

31

30

28

31

34 37 28 16 30 25 33 35 37 27 35 19

30

37

39 50 30 31 39 35 40 48 42 41 39 26

39

40

36

39

40 45 32 28 39 31 42 48 38 43 42 23

38

Appendix C. Pielou's species evenness of macrophytes in treatment and control halves in each year

Treatment

Wetland

2009

H

2 13 22 24 37 43 7 12 19 21 36 41 3 5 10 23 39 42 11 15 20 26 32 40

0.85

HN

L

LN

All wetlands

0.77

0.82

0.80

2010 0.87 0.84 0.78 0.93 0.81 0.89 0.81 0.78 0.72 0.84 0.77 0.70 0.75 0.88 0.84 0.85 0.85 0.78 0.75 0.88 0.85 0.84 0.89 0.60

0.81

0.87

0.86

0.81

0.85

0.88 0.86 0.91 0.93 0.89 0.77 0.82 0.89 0.82 0.88 0.88 0.89 0.73 0.83 0.83 0.87 0.84 0.78 0.91 0.90 0.88 0.84 0.93 0.66

0.85

0.80

0.74

0.74

0.74

2011 0.80 0.81 0.80 0.86 0.73 0.80 0.80 0.81 0.60 0.85 0.77 0.63 0.74 0.77 0.70 0.71 0.84 0.67 0.72 0.80 0.81 0.81 0.87 0.43

0.76

0.84

0.79

0.77

0.80

0.89 0.85 0.88 0.84 0.71 0.89 0.82 0.69 0.69 0.81 0.86 0.88 0.69 0.79 0.80 0.77 0.81 0.77 0.84 0.81 0.84 0.86 0.85 0.59

0.80

0.80

0.79

0.78

0.80

2012 0.83 0.86 0.76 0.85 0.68 0.80 0.86 0.81 0.72 0.85 0.82 0.69 0.84 0.86 0.80 0.66 0.86 0.66 0.81 0.88 0.79 0.84 0.93 0.51

0.79

0.86

0.84

0.81

0.85

0.90 0.86 0.91 0.89 0.70 0.92 0.81 0.85 0.74 0.89 0.89 0.82 0.81 0.89 0.87 0.70 0.86 0.73 0.89 0.90 0.86 0.91 0.92 0.64

0.84

0.82

0.78

0.78

0.81

All years 0.82 0.89 0.85 0.79 0.78 0.78 0.92 0.79 0.63 0.81 0.91 0.62 0.76 0.88 0.75 0.76 0.88 0.64 0.83 0.88 0.85 0.87 0.91 0.50

0.80

0.85

0.80

0.76

0.83

0.86 0.88 0.84 0.90 0.73 0.90 0.90 0.84 0.64 0.84 0.87 0.73 0.76 0.89 0.78 0.73 0.76 0.61 0.87 0.92 0.81 0.81 0.90 0.68

0.81

0.82

0.77

0.78

0.79

0.83 0.85 0.80 0.86 0.75 0.82 0.85 0.80 0.67 0.84 0.82 0.66 0.77 0.85 0.77 0.74 0.86 0.68 0.78 0.86 0.83 0.84 0.90 0.51

0.79

0.86

0.82

0.79

0.83

0.88 0.86 0.88 0.89 0.76 0.87 0.84 0.82 0.72 0.86 0.88 0.85 0.75 0.85 0.82 0.77 0.82 0.72 0.88 0.88 0.85 0.86 0.90 0.64

0.83

Appendix D. Sum of maximum % cover of macrophytes in treatment and control halves in each year

Treatment

Wetland

2009

H

2 13 22 24 37 43 7 12 19 21 36 41 3 5 10 23 39 42 11 15 20 26 32 40

103

HN

L

LN

All wetlands

110

97

120

107

2010 133 101 97 88 82 117 151 70 112 109 127 89 108 115 60 88 124 92 114 150 75 194 83 103

79

82

82

96

84

97 86 62 66 47 116 134 58 82 73 76 64 96 91 54 80 111 59 80 127 57 161 90 59

102

102

101

101

101

2011 98 91 94 105 87 137 153 66 100 99 95 98 108 117 85 100 87 108 107 102 99 132 73 93

589

71

75

75

70

68 50 37 65 54 76 98 51 75 104 58 39 92 84 56 58 77 82 66 86 62 115 63 59

103

106

101

102

103

11

2012 111 100 95 106 99 108 123 103 91 107 108 105 101 111 82 101 106 105 92 112 111 125 88 92

101

108

102

106

104

121 103 101 90 92 98 101 127 97 110 108 107 99 117 89 97 113 94 100 124 106 117 96 95

88

84

92

89

88

All years 91 74 74 111 79 98 90 88 96 90 46 92 91 93 86 94 86 101 94 93 87 88 81 91

78

83

90

90

85

88 75 92 67 65 81 83 75 95 89 66 90 97 85 85 89 94 92 93 102 88 94 81 82

99

100

98

103

100

108 91 90 102 87 115 129 82 99 101 94 96 102 109 78 96 100 101 102 114 93 135 81 94

79

86

87

92

86

93 78 73 72 64 93 104 78 87 94 77 75 96 94 71 81 99 82 85 110 78 122 82 74

Ecotoxicology and Environmental Safety 183 (2019) 109475

J.F. Mudge and J.E. Houlahan

Appendix E. Species detected at one or more of the wetland sites

Latin Name Acer rubrum Agrostis scabra Alisma plantago-aquatica Alnus incana Anaphalis margaritacea Aster lateriflorus Aster umbellatus Aulacomnium pallustre Betula papyrifera Betula populifolia Bidens cernua Bidens frondosa Brasenia schreberi Calamagrostis canadensis Callitriche heterophylla Calopogon tuberosus Carex canescens Carex crinite Carex echinata Carex lurida Carex magellanica Carex scoparia Carex stricta Carex tuckermanii Chara sp. Cornus canadensis Danthonia spicata Drosera rotundifolia Eleocharis acicularis Eleocharis obtusa Eleocharis ovata Eleocharis palustre Epilobium angustifolium Epilobium leptophyllum Equisetum fluviatile Equisetum pratense Eriophorum sp. Eupatorium perforiatum Fragraria vesca Fragraria virginiana Galium obtusum Gylceria borealis Glyceria striata Hieracium sp. Hypericum canadense Hypericum majus Impatiens capensis Iris versicolor Juncus brevicaudatus Juncus canadensis

Abbreviation acerub agrsca alipla alninc anamar astlat astlat aulpal betpap betpop bidcer bidfro brasch calcan calhet caltub carcan carcri carech carlur carmag carsco carstr cartuc chara corcan danspi drorot eleaci eleobt eleova elepal epiang epilep equflu equpra eri sp. eupper fraves fravir galobt glybor glystr hie sp. v hypmaj impcap iriver junbre juncan

Latin Name Juncus effusus Juncus filiformis Juncus nodosus Juncus tenuis Kalmia angustifolia Larix laricina Leersia oryzoides Lycopodium inundatum Lycopus uniflorus Lysimachia terrestris Mnium sp. Onoclea sensibilis Persicaria hydropiper Persicaria sagittata Phleum pratense Picea mariana Pinus resinosa Pinus strobus Polytrichum sp. Populus grandidentata Populus tremuloides Potamogeton epihydrus Potamogeton pusillus Potamogeton spirillus Rhododendron canadense Rhododendron groenlandicum Rosa nitida Rubus idaeus Salix bebbiana Salix eriocephala Salix pedicellaris Salix petiolaris Salix pyrifolia Scirpus cyperinus Solidago graminifolia Solidago rugosa Solidago uglinosa Sparganium americanum Sphagnum sp. Spiraea alba Spiranthes cernua Spiraea tomentosa Torreychloa pallida Toxicodendron rydbergii Triadenum fraseri Typha latifolia Utricularia intermedia Utricularia macrorhiza Veronica scutellata Viola sp.

Abbreviation juneff junfil junnod junten kalang larlar leeory lycind lycuni lyster mni sp. onosen perhyd persag phlpra picmar pinres pinstr pol sp. popgra poptre potepi potpus potspi rhocan rhogro rosnit rubida salbeb saleri salped salpet salpyr scicyp solgra solrug solugl spaame sph sp. spialb spicer spitom torpal toxryd trifra typlat utrint utrmac verscu vio sp.

Appendix F. Supplementary data Supplementary data to this article can be found online at https://doi.org/10.1016/j.ecoenv.2019.109475.

a beach on Georgian Bay. J. Great Lakes Res. 37, 616–624. Crow, G.E., Hellquist, C.B., 2000a. Aquatic and Wetland Plants of Northeastern North America: Vol. 1. Pteridophytes, Gymnosperms and Angiosperms, Dicotyledons. University of Wisconsin Press, Madison Wisconsin, pp. 480. Crow, G.E., Hellquist, C.B., 2000b. Aquatic and Wetland Plants of Northeastern North America, vol. 2. University of Wisconsin Press, Madison Wisconsin, pp. 400 Angiosperms, monocotyledons. Dalton, R.L., Boutin, C., 2010. Comparison of the effects of glyphosate and atrazine herbicides on nontarget plants grown singly and in microcosms. Environ. Toxicol. Chem. 29, 2304–2315. Edge, C.B., Thompson, D.G., Hao, C., Houlahan, J.E., 2012. A silviculture application of the glyphosate-based herbicide Visionmax to wetlands has limited direct effects on amphibian larvae. Environ. Toxicol. Chem. 31, 2375–2383. Freedman, B., Morash, R., MacKinnon D, 1993. Short-term changes in vegetation after the silvicultural spraying of glyphosate herbicide onto regenerating clearcuts in Nova Scotia, Canada. Can. J. For. Res. 25, 2300–2311. Gardner, S.C., Grue, C.E., 1996. Effects of Rodeo and Garlon 3A on nontarget wetland species in central Washington. Environ. Toxicol. Chem. 15, 441–451. Grube, A., Donaldson, D., Kiely, T., Wu, L., 2011. Pesticide Industry Sales and Usage:

References Ailstock, M.S., Norman, C.M., Bushmann, P.J., 2001. Common reed Phragmites australis: control and effects upon biodiversity in freshwater nontidal wetlands. Restor. Ecol. 9, 49–59. Anderson, D.R., Burnham, K.P., Thompson, W.L., 2000. Null Hypothesis testing: problems, prevalence, and an alternative. J. Wildl. Manag. 64, 912–923. Anderson, E.W., 1986. A guide for estimating cover. Rangelands 8, 236–238. Burge, O.R., Bodmin, K.A., Clarkson, B.R., Bartlam, S., Watts, C.H., Tanner, C.C., 2017. Glyphosate redirects wetland vegetation trajectory following willow invasion. Appl. Veg. Sci. 20, 620–630. Byer, J.D., Struger, J., Klawunn, P., Todd, A., Sverko, E., 2008. Low cost monitoring of glyphosate in surface waters using the ELISA method: an evaluation. Environ. Sci. Technol. 42, 6052–6057. Cohen, J., 1988. Statistical Power Analysis for the Behavioral Sciences, second ed. Lawrence Earlbaum Associates Inc., Hillsdale New Jersey, pp. 569. Crowe, A.S., Leclerc, N., Struger, J., Brown, S., 2011. Application of a glyphosate-based herbicide to Phragmites australis: impact on groundwater and near-shore lake water at

12

Ecotoxicology and Environmental Safety 183 (2019) 109475

J.F. Mudge and J.E. Houlahan

Paveglio, F.L., Kilbride, K.M., 2000. Response of vegetation to control of reed canarygrass in seasonally managed wetlands of southwestern Washington. Wildl. Soc. Bull. 28, 730–740. Rubin, J.L., Gaines, C.G., Jensen, R.A., 1982. Enzymological basis for herbicidal action of glyphosate. Plant Physiol. (Wash. D C) 70, 833–839. Struger, J., Thompson, D., Staznik, B., Martin, P., McDaniel, T., Marvin, C., 2008. Occurrence of glyphosate in surface waters of Southern Ontario. Bull. Environ. Contam. Toxicol. 80, 378–384. Sullivan, T.P., Sullivan, D.S., 2003. Vegetation management and ecosystem disturbance: Impact of glyphosate herbicide on plant and animal diversity in terrestrial systems. Environ. Rev. 11, 37–59. Vitta, J.I., Tuesca, D., Puricelli, E., 2004. Widespread use of glyphosate tolerant soybean and weed community richness in Argentina. Agricult. Ecosyst. Environ. Times 103, 621–624. Wetzel, R., 2001. Limnology: Lake and River Ecosystems. Academic Press, Orlando Florida. Williams, G.M., Kroes, R., Munro, I.C., 2000. Safety evaluation and risk assessment of the herbicide Roundup and its active ingredient, glyphosate, for humans. Regul. Toxicol. Pharmacol. 31, 117–165. Woodburn, A.T., 2000. Glyphosate: production, pricing and use worldwide. Pest Manag. Sci. 56, 309–312. Zhou, Q., Lin, Y., Li, X., Han, Z., Zeng, G., Lu, L., He, S., 2018. Effect of zinc ions on nutrient removal and growth of Lemna aequinoctialis from anaerobically digested swine wastewater. Bioresour. Technol. 249, 457–463. Zhou, Q., Li, X., Lin, Y., Yang, C., Tang, W., Wu, S., Li, D., Lou, W., 2019. Effect of copper ions on removal of nutrients from swine wastewater and on release of dissolved organic matter in duckweed systems. Water Res. 158, 171–181.

2006 and 2007 Market Estimates. Biological and Economic Analysis Division, Office of Pesticide Programs, Office of Chemical Safety and Pollution Prevention. U.S. Environmental Protection Agency, Washington DC, USA. Horsley, S.B., 1994. Regeneration success and plant diversity of Allegheny hardwood stands after Roundup application and shelterwood cutting. N. J. Appl. For. 11, 109–116. Jadhav, A., Hill, M., Byrne, M., 2008. Identification of a retardant dose of glyphosate with potential for integrated control of water hyacinth, Eichhornia crassipes (Mart.) Solms-Laubach. Biol. Control 47, 154–158. Kielak, E., Sempruch, C., Miodusyewska, H., Klocek, J., Leszcyznski, B., 2011. Phytotoxicity of Roundup Ultra 360 SL in aquatic ecosystems: biochemical evaluation with duckweed (Lemna minor L.) as a model plant. Pestic. Biochem. Physiol. 99, 237–243. Linz, G.M., Homan, H.J., 2011. Use of glyphosate for managing invasive cattail (Typha spp.) to disperse blackbird (Icteridae) roosts. Crop Protect. 30, 98–104. Magurran, A.E., 2004. Measuring Biological Diversity. Blackwell Publishing Company, MA, USA. Mateos-Naranjo, E., Perez-Martin, A., 2013. Effect of sub-lethal glyphosate concentrations on growth and photosynthetic performance of non-target species Bolboschoenus matitimus. Chemosphere 93, 2631–2638. Mozdzer, T.J., Hutto, C.J., Clarke, P.A., Field, D.P., 2008. Efficacy of imazapyr and glyphosate in the control of non-native Phragmites australis. Restor. Ecol. 16, 221–224. Mudge, J.F., Baker, L.F., Edge, C.B., Houlahan, J.E., 2012. Setting an optimal α that minimizes errors in null Hypothesis significance tests. PLoS One 7, e32734. Munkittrick, K.R., Arens, C.J., Lowell, R.B., Kaminski, G.P., 2009. A review of potential methods of determining critical effect size for designing environmental monitoring programs. Environ. Toxicol. Chem. 28, 1361–1371.

13