Etobicoke Creek, Toronto, Ontario, Canada

Introduction

As per previous investigations (Bazinet et al 2010), there is evidence to suggest that the increasing urbanization in the Greater Toronto Area correlates with a decrease in biotic community population and diversity across affected lotic systems across time. Therefore, as part of this investigation, a thorough analysis of watersheds extending from Etobicoke to Ajax was conducted and subsequently, assessed in conjunction with the previous results and surveys (Wallace et al 2013).

library(tidyverse)
library(ggplot2)
library(readr)
library(SciViews)
library(lubridate) 
library(vegan)
library(ggpubr)
library(gghighlight)
#loading datasets and necessary packages
#Data from the TRCA’s Regional Watershed Monitoring Program
Counts_2013 <- read_csv("rwmp-benthic-macroinvertebrate-data-2013-2020.csv") #Loading data and tools, Using data from Toronto Regional and Conservation Authority
## Warning: One or more parsing issues, see `problems()` for details
#Note that the stats all use Kick and Sweep Method and 100+ fixed count method and thus, are removed from stats for simplicity
#also removing site name 
Counts <- Counts_2013[,-c(1,3:4, 19:20)]

First, we establish the datasets we will be using, of which we have many. The Toronto and Region Council Authority (TRCA) website has open access data as far back as 2013 for our intended purposes. With this, we have completed an detailed analysis of macroinvertebrate diversity with a focus on sensitive taxa as well as potential environmental stressors including temperature and dissolved oxygen using the appropriate reference guides and methodologies where possible. The analysis techniques used in this report were derived from previously utilized methods in order to maintain relative comparability. Following, a brief discussion regarding the results of current and past trends is provided. Finally, a conclusion to round out how to improve the watersheds, if necessary will be brought to attention.

As per the data that was taken from the Toronto and Region Conservation Authority (TRCA), figure 1 shows the testing sites across the Toronto region for the data taken and was created via QGIS.

Figure 1. Map of the Greater Toronto Region and recorded macroinvertebrate survey sites between 2013 and 2020.

Macroinvertebrate Counts

Using the data over the period 2013-2020, Shannon-Wiener index will be used to determine the diversity of macroinvertebrate communities in the GTA Region.

\[ H = -\sum_{i = 1}^{S}{P_iln(P_i)} \] Figure 2. Shannon-Weiner Index for the surveyed creeks in the Greater Toronto Area, 2013 - 2020

H = Shannon diversity index

S = Number of Species Present

P = proportion (n/N) of individuals of a species (n) found divided by the total number of individuals found (N)

ln = Natural Log

E = Summation of the following equation

Counts <- Counts[,-c(2:5,7:9,11,14,15)]  #cleaning unnecessary columns 
Counts <- Counts %>% group_by(Watershed, `Collection Date`) %>% mutate(Prop = `Total Count`/sum(`Total Count`)) #Calculating P 
Counts <- Counts %>% group_by(Watershed, `Collection Date`) %>% mutate(lnProp = (ln(Prop))) #Calculating natural log
Counts <- Counts %>% group_by(Watershed, `Collection Date`) %>% mutate(PLP = Prop*lnProp) #Multiplying the proportions by natural log of the proportions
Counts <- Counts %>% group_by(Watershed, `Collection Date`) %>% mutate(Shannon_final = abs(sum(PLP))) #summarizing the shannon values by  watershed and then absolute valuing it 
## Clean Display of Data
Counts %>% group_by(`Watershed`, `Collection Date`) %>% summarise(Shannon = diversity(`Total Count`)) -> Test

Figure 3. Scatterplots depicting Shannon-Wiener Index across tested watersheds in the GTA during the period 2013-2020. NB:The higher the number, the more diverse the community.

As seen here, a variety of trends can be depicted across different systems in the GTA. The most definitive conclusion we can draw from this is that there are upward trends towards higher diversity in macroinvertebrates within such systems. For example, systems such as the Humber River and Petticoat Creek show an increasing trend in diversity over the last decade. However, we have not looked into less tolerant taxa such as tricoptera, plecoptera and ephemeroptera who serve as better overall ecological indicators than species such as Oligochaeta and Culicidae for instance (Thorne et al 2000).

Sensitive Taxa

#Counting the number of sensitive taxa (Mayfly, Stonefly, Caddisfly)
subset(Counts, grepl("Mayfly", `Common Name`)) -> Mayfly #Found em i think
subset(Counts, grepl("Stonefly",`Common Name`)) -> Stonefly
subset(Counts, grepl("Caddisfly", `Common Name`)) -> Caddisfly
rbind(Mayfly, Caddisfly, Stonefly) -> End
End %>% group_by(Watershed, `Collection Date`) %>% summarise("Total Sum" = sum(`Total Count`)) %>% ggplot(aes(`Collection Date`, `Total Sum`)) + facet_wrap(~Watershed) + geom_jitter()

Figure 4. Scatterplot for distribution and abundance of sensitive taxa, tricoptera, plecoptera and ephemeroptera

As alluded to by various sources, whilst some river systems display varying abundances of the EPT taxa total population, systems such as Carruthers Creek, Frenchman’s Bay and Petticoat Creek show significantly lower population or near absent populations. In addition, these underperforming systems are located primarily in Pickering and Ajax in the Durham Region. There may be some correlation between the locations and their current land management.

Whilst overt trends are not so present in the data, potential trends may be more present in environmental abiotic data. Therefore, a thorough data analysis of such is also conducted at the following sites in figure 2. In accordance with Qu et al 2005 and Jiang et al 2010, Elevation and levels of silicon dioxide are primary factors determining macroinvertebrate populations. Hence, these factors warrant investigation. Other studies have shown negative relationships between urbanization effects and macroinvertebrate populations such as increased chemical input into such systems and a decrease in diversity, which is more likely the case for watersheds in our area of study (Wallace et al 2013).

For the sake of simplicity, zinc, temperature and conductivity will be calculated and compared across the streams in this study.

Dissolved Oxygen

With benthic macroinvertebrates, dissolved oxygen is of the utmost importance for such organisms. It serves as an indicator for the presence of benthic macrophytes and the lack thereof can be a warning that can give rise to algal blooms and other wanton effects (Ogbeibu and Oribhabor 2002). Thus, such a parameter is very important to look at in terms of our it may affect communities across the GTA.

#Loading data from csv files
read.csv("2015-wq-data-from-database-no-project-sites.csv") -> Data2015
read.csv("2016-wq-data-from-database-no-project-sites.csv") -> Data2016
read.csv("trca-2017-wq-data-rwmp-sites.csv") -> Data2017
read.csv("trca-2018-wq-data-rwmp-sites.csv") -> Data2018
read.csv("trca-2019-wq-data-rwmp-sites.csv") -> Data2019
read.csv("2020-wq-data-from-database-no-project-sites-1.csv") -> Data2020
#Extracting Dissolved Oxygen Data and removing uneccessary columns
Data2015[,c(5,10:12,13, 15:16)] %>% filter(ParameterDefinitionName == "Oxygen, Dissolved") -> DO15
Data2016[,c(5,10:12,13, 15:16)] %>% filter(ParameterDefinitionName == "Oxygen, Dissolved") -> DO16
Data2017[,c(3,8:11, 13:14)] %>% filter(ParameterDefinitionName == "Oxygen, Dissolved") -> DO17
Data2018[,c(5,10:12,13, 15:16)] %>% filter(ParameterDefinitionName == "Oxygen, Dissolved") -> DO18
Data2019[,c(5,10:12,13, 15:16)] %>% filter(ParameterDefinitionName == "Oxygen, Dissolved") -> DO19
Data2020[,c(5,10:12,13, 15:16)] %>% filter(ParameterDefinitionName == "Oxygen, Dissolved") -> DO20
#Coagulating all the data together
rbind(DO15,DO16, DO17, DO18, DO19, DO20) -> DO
#Graphing it
ggplot(DO) + geom_jitter(aes(Year, FinalValue))  + gghighlight::gghighlight(FinalValue < 6.5) + facet_wrap(~Watershed) + ylim(0,15)
## Warning: Using `across()` in `filter()` is deprecated, use `if_any()` or
## `if_all()`.
## Warning: Removed 559 rows containing missing values (geom_point).
## Warning: Removed 32 rows containing missing values (geom_point).

Figure 5. Scatterplot for dissolved oxygen (DO) values across the GTA over the period 2015-2020. Data highlighted indicates below expected DO values (< 6.5mg/L).

6.5mg/L was used as the baseline per recommendations from Environment Canada (Unknown). Hence, they are highlighted in Figure 5. Here we can see that overall, dissolved oxygen is not necessarily below the norm outside of a few points in time. Notably, Etobicoke Creek and Rouge River have multiple points in which the DO levels were lower than the threshold and may prove to be causes for concern.

Electrical Conductivity

As electrical conductivity serves as an indicator for both dissolved metals as well as salinity in watersheds, it is important to look at such a metric in addition to other parameters in order to develop a more complete picture of the ecological health of our waterways.

Data2015[,c(5,10:12,13, 15:16)] %>% filter(ParameterDefinitionName == "Conductivity, Specific (Field)") -> CON15
Data2016[,c(5,10:12,13, 15:16)] %>% filter(ParameterDefinitionName == "Conductivity, Specific (Field)") -> CON16
Data2017[,c(3,8:11,13:14)] %>% filter(ParameterDefinitionName == "Conductivity, Specific (Field)") -> CON17
Data2018[,c(5,10:12,13, 15:16)] %>% filter(ParameterDefinitionName == "Conductivity, Specific (Field)") -> CON18
Data2019[,c(5,10:12,13, 15:16)] %>% filter(ParameterDefinitionName == "Conductivity, Specific (Field)") -> CON19
Data2020[,c(5,10:12,13, 15:16)] %>% filter(ParameterDefinitionName == "Conductivity, Specific (Field)") -> CON20
rbind(CON15, CON16, CON17) -> CON #combining it all together
rbind(CON18, CON, CON19, CON20) -> CON
CON %>% arrange(Year) -> CON #sorting out by year
CON %>% ggplot(aes(Year, FinalValue)) + geom_jitter() + facet_wrap(~Watershed) + xlab("Year") + ylab("Microsiemens (uS/cm)") + ggtitle(label = "Figure 6A" )  #jitter plot for Conductivity 
## Warning: Removed 35 rows containing missing values (geom_point).

CON %>% ggplot(aes(Year, FinalValue)) + geom_jitter() + facet_wrap(~Watershed) + ylim(0, 5000) + xlab("Year") + ylab("Microsiemens (uS/cm)") + gghighlight::gghighlight(FinalValue > 1000) + ggtitle(label = "Figure 6B")  #jitter plot again but with thresholds put in place and highlighted
## Warning: Using `across()` in `filter()` is deprecated, use `if_any()` or
## `if_all()`.
## Warning: Removed 206 rows containing missing values (geom_point).
## Warning: Removed 169 rows containing missing values (geom_point).

Figure 6a and b. a) displays a electrical conductivity surveys across watersheds in the Greater Toronto Area. b) shows all watersheds with results that exceed the reference levels for Ontario (1000 uS/cm).

In accordance with the Environmental Protection Agency (EPA, 2012), ideal electrical conductivity reference levels are between 0-500 uS/cm. In that regard, half of all the tests show values greater than that. This is particularly apparent for Don River, Etobicoke Creek and Humber where these values are regularly exceeded. However, values from Environment Canada state that anything below 1000uS/cm is considered acceptable an even expected of most rivers.

With regards to figure 6B, each assessed waterway displays high electrical conductivity across a number of points each year. That being said, there does not seem to be a temporal trend so to speak of but alas, such a commonplace trend across all the rivers is cause for concern and should be remedied where possible.

Zinc Levels

Heavy metals such as zinc have the potential to bioaccumulate in macroinvertebrates. In some cases, levels of trace metals found in surveyed organisms are equivalent to those found in nearby substrate (Goodyear and McNeil, 1999). Hence, Zinc levels in the water have been assessed accordingly. Going off surface water concentration guidelines set out by the Canadian Water Quality Guidelines for the Protection of Aquatic Life (2018), acceptable concentrations seem to sit between <2 to 537ug/L. Hence, trace amounts are expected to occur in any given waterbodies.

#loading zinc data
Data2015 %>% filter(ParameterDefinitionName == "Zinc, Total") -> ZN2015 #Filtering data
Data2016 %>% filter(ParameterDefinitionName == "Zinc, Total") -> ZN2016
Data2017 %>% filter(ParameterDefinitionName == "Zinc, Total") -> ZN2017
Data2018 %>% filter(ParameterDefinitionName == "Zinc, Total") -> ZN2018
Data2019 %>% filter(ParameterDefinitionName == "Zinc, Total") -> ZN2019
Data2020 %>% filter(ParameterDefinitionName == "Zinc, Total") -> ZN2020
ZN2015[,-c(1:4,6:9,13:14,17:22)] -> ZN2015 #eliminating unnecessary data
ZN2016[,-c(1:4,6:9,13:14,17:22)] -> ZN2016
ZN2017[,-c(1:2,4:7,11:12,15:19)] -> ZN2017
ZN2018[,-c(1:4,6:9, 13:14,17:22)] -> ZN2018
ZN2019[,-c(1:4,6:9, 13:14,17:22)] -> ZN2019
ZN2020[,-c(1:4,6:9,13:14,17:25)] -> ZN2020
rbind(ZN2015, ZN2016, ZN2017) -> ZN #combining it all together
rbind(ZN2018, ZN, ZN2019, ZN2020) -> ZN
ZN %>% ggplot(aes(Year, FinalValue)) + geom_jitter() + facet_wrap(~Watershed) + xlab("Year") + ylab("µg/L")  #jitter plot for zinc
## Warning: Removed 34 rows containing missing values (geom_point).

Figure 7. Zinc concentrations for watersheds in the Greater Toronto Area.

Figure 7 shows us that zinc levels across the GTA are within acceptable levels. Of note, results in Duffins Creek and Etobicoke Creek display single anomalies that are close to the 537ug/L but at this point in time, there seems to be no significant trend so to speak of.

Temperature

Whilst water temperature may not necessarily correlate directly with macroinvertebrate health and population diversity within a stream (Chinnayakanhalli et al 2011), it is worth considering due to its property to directly correlate with other factors such as dissolved oxygen levels and changes in local climate (EPA, 2016).

#Loading temperature data
read.csv("rwmp_watertemperature_data_2014.csv") -> Temp_2014
Temp_2014[,-c(1:7,12:19,24:26)] -> Temp_2014
Temp_2014 %>% group_by(Watershed, Month_Number, Day_Number) %>% summarise(n_temp = mean(Obs_Value)) -> Temp_2014 
#Need to do for the rest of the years
Year <- c(rep(2014))
#2015
read.csv("rwmp_watertemperature_data_2015.csv") -> Temp_2015
Temp_2015[,-c(1:7,11:17,23:25)] -> Temp_2015
Temp_2015[,-c(2:3)] -> Temp_2015
Temp_2015 %>% group_by(Watershed, Month_Number, Day_Number) %>% summarise(n_temp = mean(Obs_Value)) -> Temp_2015
#2016
read.csv("rwmp_watertemperature_data_2016.csv") -> Temp_2016
Temp_2016[,-c(1:7,9:17,23:25)] -> Temp_2016
Temp_2016 %>% group_by(Watershed, Month_Number, Day_Number) %>% summarise(n_temp = mean(Obs_Value)) -> Temp_2016
#2017
read.csv("rwmp_watertemperature_data_2017.csv") -> Temp_2017
Temp_2017[,-c(1:7,9:17,23:25)] -> Temp_2017
Temp_2017 %>% group_by(Watershed, Month_Number, Day_Number) %>% summarise(n_temp = mean(Obs_Value)) -> Temp_2017
#2018
read.csv("rwmp_watertemperature_data_2018_a.csv") -> Temp_2018a
Temp_2018a[,-c(1:7,9:17,23:25)] -> Temp_2018a
Temp_2018a %>%  group_by(Watershed, Month_Number, Day_Number) %>% summarise(n_temp = mean(Obs_Value)) -> Temp_2018a
read.csv("rwmp_watertemperature_data_2018_b.csv") -> Temp_2018b
Temp_2018b[,-c(1:7,9:17,23:25)] -> Temp_2018b
Temp_2018b %>% rename(Watershed = Humber, Timetag = X2018.05.31.23.00.00.000, Year_Number = X2018, Month_Number = X5, Day_Number = X31, Obs_Value = X19.912) -> Temp_2018b
Temp_2018b %>%  group_by(Watershed, Month_Number, Day_Number) %>% summarise(n_temp = mean(Obs_Value)) -> Temp_2018b
read.csv("rwmp_watertemperature_data_2018_c.csv") -> Temp_2018c
Temp_2018c %>% rename(Watershed = Rouge, Timetag = X2018.08.25.00.15.00.000, Year_Number = X2018, Month_Number = X8, Day_Number = X25, Obs_Value = X21.378) -> Temp_2018c
Temp_2018c  %>%  group_by(Watershed, Month_Number, Day_Number) %>% summarise(n_temp = mean(Obs_Value)) -> Temp_2018c
rbind(Temp_2018a, Temp_2018b, Temp_2018c) -> Temp_2018
#Altogether
Temp_2014 %>% mutate(Year = 2014) -> Temp_2014
Temp_2015 %>% mutate(Year = 2015) -> Temp_2015
Temp_2016 %>% mutate(Year = 2016) -> Temp_2016
Temp_2017 %>% mutate(Year = 2017) -> Temp_2017
Temp_2018 %>% mutate(Year = 2018) -> Temp_2018
rbind(Temp_2014, Temp_2015, Temp_2016, Temp_2017, Temp_2018) -> Temp_Full
#Trying something new
Temp_2014 %>% group_by(Watershed, Month_Number) %>% summarise(N = mean(n_temp)) -> Simp2014
Temp_2015 %>% group_by(Watershed, Month_Number) %>% summarise(N = mean(n_temp)) -> Simp2015
Temp_2016 %>% group_by(Watershed, Month_Number) %>% summarise(N = mean(n_temp)) -> Simp2016
Temp_2017 %>% group_by(Watershed, Month_Number) %>% summarise(N = mean(n_temp)) -> Simp2017
Temp_2018 %>% group_by(Watershed, Month_Number) %>% summarise(N = mean(n_temp)) -> Simp2018
#Cleaner graph with lines and no dodgy stuff
ggplot(Simp2014, aes(Month_Number,N)) + geom_line(aes(colour = "2014")) + geom_line(data = Simp2018, aes(color = "2018")) + facet_wrap(~Watershed)

Figure 8. Temperature changes over the period 2014-2018 in watersheds across the Greater Toronto Area.

Across the watersheds, a majority of them show an increase between the years 2014-2018 in terms of water temperature changes. On average, the summer temperatures display a higher increase than the winter temperatures. Albeit, some of the data is missing from certain sites. This is certainly a case that should be investigated further as higher temperatures generally indicates a reduction in dissolved oxygen for example. This in turn will lead to a decrease in macroinvertebrate populations.

Conclusion

Whilst macroinvertebrate communities show a wide variety of diverse taxa in the GTA, the effects of urbanized streams become more apparent when certain watersheds show a lack of sensitive taxa as well as increasing temperatures and high levels of electrical conductivity. In order to tackle these problems, a number of solutions can be implemented which when compounded can reduce the effects of urban development on these systems. Riparian plantings can provide much needed shade as well as filtering out heavy metals travelling in soil. Additionally, for the most degraded environments, more active volunteering efforts to clean the perimeter of their respective watersheds such as Don River, Etobicoke and Humber will greatly benefit the ecology of these areas.

References:

Bazinet, N. L., Gilbert, B. M., & Wallace, A. M. (2010). A comparison of urbanization effects on stream benthic macroinvertebrates and water chemistry in an urban and an urbanizing basin in Southern Ontario, Canada. Water Quality Research Journal, 45(3), 327-341.

Canadian Water Quality Guidelines for the Protection of Aquatic Life. (2018). Accessed April 8 2022 from https://ccme.ca/en/res/zinc-en-canadian-water-quality-guidelines-for-the-protection-of-aquatic-life.pdf.

Chinnayakanahalli, K. J., Hawkins, C. P., Tarboton, D. G., & Hill, R. A. (2011). Natural flow regime, temperature and the composition and richness of invertebrate assemblages in streams of the western United States. Freshwater Biology, 56(7), 1248-1265.

Environment and Natural Resources Canada. (Date Unknown). Accessed July 28 2022 from https://www.enr.gov.nt.ca/sites/enr/files/conductivity.pdf.

Environment and Natural Resources Canada. (Date Unknown). Accessed April 10 2022 from https://www.enr.gov.nt.ca/sites/enr/files/dissolved_oxygen.pdf

Environmental Protection Agency. EPA (2012). Accessed April 9 2022 from https://archive.epa.gov/water/archive/web/html/vms59.html.

Environmental Protection Agency. EPA (2016). Accessed July 28 2022 from https://www.epa.gov/sites/default/files/2016-08/documents/print_stream-temperature-2016.pdf

Goodyear, K. L., & McNeill, S. (1999). Bioaccumulation of heavy metals by aquatic macro-invertebrates of different feeding guilds: a review. Science of the Total Environment, 229(1-2), 1-19.

Jiang, X. M., Xiong, J., Qiu, J. W., Wu, J. M., Wang, J. W., & Xie, Z. C. (2010). Structure of macroinvertebrate communities in relation to environmental variables in a subtropical Asian river system. International Review of Hydrobiology, 95(1), 42-57.

Ogbeibu, A. E., & Oribhabor, B. J. (2002). Ecological impact of river impoundment using benthic macro-invertebrates as indicators. Water research, 36(10), 2427-2436.

Ontario Ministry for the Environment. (2012). Accessed March 20 2022 from https://www.ontario.ca/document/water-quality-15-streams-agricultural-watersheds-southwestern-ontario-2004-2009#section-5.

Qu, X., Tang, T., Xie, Z., Ye, L., Li, D., & Cai, Q. (2005). Distribution of the macroinvertebrate communities in the Xiangxi River system and relationships with environmental factors. Journal of Freshwater Ecology, 20(2), 233-238.

Thorne, R. S. J., Williams, W. P., & Gordon, C. (2000). The macroinvertebrates of a polluted stream in Ghana. Journal of Freshwater Ecology, 15(2), 209-217.

Toronto Regional Conservation Authority. (2021).Retrieved February 27 2022 from https://data.trca.ca/dataset?q=TRCA%E2%80%99s+Regional+Watershed+Monitoring+Program.

Wallace, A. M., Croft-White, M. V., & Moryk, J. (2013). Are Toronto’s streams sick? A look at the fish and benthic invertebrate communities in the Toronto region in relation to the urban stream syndrome. Environmental monitoring and assessment, 185(9), 7857-7875.