The Great Lakes Basin is characterized by a connected set of five lakes, with generally comparable climatic conditions (Scott and Huff, 1997; McDermid et al., 2015), and varied urban development conditions, creating the ideal region for applying and testing the idea that species richness accrual is faster in downstream watershed positions – even in urban watersheds. The Canadian side of the GLB, includes streams that drain into Lake Superior, Lake Huron, Lake Erie, and Lake Ontario (Figure 1). This study does not include streams or land areas from the United States side of the GLB because sampling efforts and methods differed across these states, and it was important to have consistent and comparable samples. The Great Lakes Basin hydrology layers were derived from the Great Lakes – St. Lawrence River Basin Watersheds (GLBW) package (Ontario Ministry of Natural Resources and Forestry – Provincial Mapping Unit, 2009). We further divide our study area into the Ecoregions that capture different portions of the GLB (Table 1). Human activities and development are influenced by the climatic conditions defined within Ecoregions. The northwest of Ontario is largely forested and has an historic forestry industry presence (Crins et al., 2009). The population density of this region is also relatively low resulting in lower incidences of urbanized landcover (Crins et al., 2009). To the southeast, where the climate is generally milder and the growing season is substantially longer, there is a high percentage of agriculture and urbanized landcover (Crins et al., 2009).

The Ontario Landcover Compilation Version 2 (OLCCv2) was used to define land cover classes in the GLB. The OLCCv2 is a combination of the following landcover packages: the Far North Landcover Version 1.4 (FNLCv1.4), the Southern Ontario Land Resource Information System, and The Provincial Landcover 2000 Edition (PLC2000) (Ontario Ministry of Natural Resources and Forestry, 2014). The OLCCv2 contained Landsat imagery for the entire Province of Ontario taken between the years 1999–2011: 1999–2002 (PLC2000), 1999–2002 (SOLRIS) and 2005–2011 (FNLCv1.4, Ontario Ministry of Natural Resources and Forestry, 2014). The standardized resolution of the OLCCv2 is fifteen (15) meters and the compilation uses the projection, NAD83 Lambert Conformal Conic (Ontario Ministry of Natural Resources and Forestry, 2014). The OLCCv2 contains twenty-nine (29) landcover classes, which are listed fully in Supplementary Table S1. The OLCCv2 was downloaded from Ontario GeoHub.¹

Using the “reclassify” tool in ArcMap, the initial 29 landcover classes were condensed into six classes for the purpose of developing separate species accrual relationships by cover class. The six classes were: Urban, Greenspace, Agriculture, Forest, Water, and an additional class labeled Other (see Supplementary Table S1). “Urban” landcover defines built-up areas and their associated anthropogenic structures (Ministry of Natural Resources and Forestry Science and Research Branch, 2019). “Greenspace” typically includes all public green areas, including urban parks, wetlands, and recreational facilities (Hahs and McDonnell, 2006; Taylor and Hochuli, 2017; Suarez-Rubio and Krenn, 2018). We use the term “Greenspace” to represent natural landcover that did not fall into the definition of Forest (e.g., Hedgerows, Disturbance, Tallgrass Prairie, etc.): this landcover class included the most varied of the original OLCCv2 landcover classes. The “Agriculture” land class consisted of all agricultural areas including semi-natural areas (Ministry of Natural Resources and Forestry Science and Research Branch, 2019). “Forest” included all areas with tree cover greater than 60% (Lee et al., 1998). “Water” refers to water bodies such as lakes (Ministry of Natural Resources and Forestry Science and Research Branch, 2019), but for the purpose of this study streams were included implicitly by looking directly at the potential impact of watershed position. The “Other” category was used to classify any features that did not fall under any of the five categories (i.e., Bedrock, Sand, Gravel, and Cloud or Shadow, Supplementary Table S1).

In this study, “stream” refers to any size of river or flowing body of freshwater, and these are characterized using varying Strahler orders (Figure 2; Strahler, 1957). Ontario hydrology datasets and layers (within quaternary watersheds) were based on the Ontario Integrated Hydrology data (OIH; Provincial Mapping Unit, 2016). Stream orders were grouped to ensure adequate bird point count sample sizes in downstream areas as there are more headwater areas than mid-order and downstream ones. Watershed positions (WSP) were attributed using the Strahler order classification of stream reaches as follows:

First and second-order streams that drained directly into the Great Lakes were removed. Their exclusion was necessary to ensure that these smaller stream areas do not skew the results due to lake effects, as they would be classified as headwaters by the study criteria, but they are located, and possibly behaving like, downstream WSPs given their position relative to the Great Lakes.

All GIS processing was completed using ArcMap 10.7.1. Using the “Clip” tool in ArcMap, the enhanced stream flow layer from the OIH was cropped to the extent of the GLB study area. Streamflow lines were then “Dissolved” so the only identifying information was Strahler order. This simplification allowed for the compilation of numerous enhanced flow datasets to be compared. The “Multiplepart to Singlepart” feature was then applied to separate streams of similar orders into independent features allowing for the selection of these streams individually. All separate enhanced flow layers were then combined using the “Merge” feature. The merged enhanced flow of the GLB was then used to select the boundaries located on the shore of the Great Lakes. Headwater streams that flow directly into the lakes were then spatially selected and removed from analysis.

Point count data were obtained from the Ontario Breeding Bird Atlas (OBBA) 2001–2005. This was the second such Atlas compiled for the Province of Ontario. The first OBBA covered 1981–1985 and the upcoming 3rd OBBA surveys began in 2021. The OBBA dataset contains species observations expected to represent most breeding birds within Ontario. Sampling counts occurred between May and July (dates varied depending on location) of the years 2001 through 2005 (Bird Studies Canada, Environment Canada’s Canadian Wildlife Service, Ontario Nature, Ontario Field Ornithologists and Ontario Ministry of Natural Resources, 2008). Surveys for the OBBA were conducted within 10 × 10 km squares covering the province (Bird Studies Canada, Environment Canada’s Canadian Wildlife Service, Ontario Nature, Ontario Field Ornithologists and Ontario Ministry of Natural Resources, 2008). This covered the regions within the Great Lakes Basin Watershed. Each point count was observed for 5 min and all bird species within a 100 m radius were recorded by birders able to identify birds by their song (Bird Studies Canada, Environment Canada’s Canadian Wildlife Service, Ontario Nature, Ontario Field Ornithologists and Ontario Ministry of Natural Resources, 2008). A minimum of twenty-five (25) of these observations were completed per square by volunteer birders (for the majority of the OBBA) to ensure that the abundance of each species was well represented, but a minimum of 10-point counts per square was accepted in the final year (Cadman et al., 2007). Finally, all observations were validated by professional ornithologists and Regional Coordinators of the OBBA sampling team (Bird Studies Canada, Environment Canada’s Canadian Wildlife Service, Ontario Nature, Ontario Field Ornithologists and Ontario Ministry of Natural Resources, 2008).

The OBBA dataset including rare species was retrieved via a data request to the Nature Counts portal. The dataset was then modified and cleaned in R version 4.0.3 (R Core Team, 2020) before being assessed further in ArcMap. Point counts were initially located by latitude and longitude, which were then transformed and projected to NAD 1983 Ontario Lambert Conformal Conic. For this study, we selected point counts that were located within the limits of the GLB. This was completed using the “Clip” tool in ArcMap. In some instances, a particular point was sampled on multiple occasions; for these instances, the sample with the highest total species richness was selected to represent that point count.

We used the observation distance radius, i.e., 100 m to characterize the underlying landcover around point count sampling locations. The majority landcover (LCMAJ) of each of these buffers was then calculated using the “zonal statistics” tool in ArcMap. This provided a more accurate landcover description of the surveyed area than that of the underlying landcover values at the exact point count location. The distance of each buffer to the nearest stream was also included using the “near” tool. This feature was used to identify the Strahler order associated with that sample point.

Upon completion of our spatial analysis, each point count contained information on the following:

Point count data were divided into subsets based on Ecoregion-Landcover Classification-Watershed position (ER-LCC-WSP) using the “subset” function in R. From here on we will refer to these subsets as subgroups. The number of sites in each subgroup was evaluated using the “BiodiversityR” package and the function “diversitycomp,” which provided a count of sites at each subgrouping (Kindt and Coe, 2005). The number of sites located in the OUT WSP in each subgroup had the smallest number of sites per subgroup category. The smallest number was then used as the number of samples to create our Species Accumulation Curves (SACs). Through this method, we assessed each WSP in a specific subgrouping as evenly as possible to glean a better understanding of the impact of watershed position. Additionally, species diversity metrics including richness, abundance, Shannon, and Simpson were calculated for each subgroup.

The first steps in creating the SACs involved using the “vegan” package function: specaccum (Oksanen et al., 2020). The “exact” method, which provides the expected mean richness for a specified number of sampling sites, was used (Kindt and Coe, 2005). Using both AIC and BIC to test model fit, the Lomolino model outperformed other SAC models (i.e., Gleason, Gitay, and Arrhenius) in most of the curves created. The resultant parameters of Lomolino fitted SACs are in the form of a slope, an asymptote, and an xmid. Slopes of the Lomolino models represents the slope at the inflection point, or the point on the slope where the concavity of the curve changes (Lomolino, 2001). In terms of the Lomolino model, this inflection point also represents the maximum rate of increase in species richness (Lomolino, 2001). Fitted slope is therefore the parameter that allows for estimating the rate of species accrual at its greatest point. The asymptote represents the maximum richness of the sites sampled (Lomolino, 2001). The xmid is the number of points at which half the total species richness is achieved, this parameter however was not included for further analysis.

Upper and lower confidence intervals around these Lomolino model parameters (slopes and asymptotes) were then used to create a range of possible curves and sampled 1,000 times. From these repeated samples, we extracted the curve coefficients, i.e., slopes and asymptotes as the response variable in our statistical models, but accrual asymptotes greater than the observed 216 species (total species count) were removed. The 1,000 SAC iterations increased the likelihood of covering the full range of sampling possibilities from the existing community.

Using mixed effect models (MEM), we evaluated the relationship between our response variables, i.e., slope and asymptote, respectively, and the predictor variables, i.e., landcover class (LCC), watershed position (WSP), and Ecoregion in which the points making up the SACs were located. Species accrual slopes or asymptotes were in the form of positive nonzero numbers, therefore the Gamma distribution with a log link was used (Bolker et al., 2009; Harrison et al., 2018). The use of a mixed model approach allowed us to estimate the fixed effects of land cover class (LCC) and watershed position (WSP) on our response variables while controlling for Ecoregion as a random effect influencing both slopes and intercepts (see GLMM section below) allowing us to characterize random variation caused by the geographical Ecoregion of sample sites. Several models were tested and compared based on their AIC using the “anova” function. The most parsimonious models that did not include confounded variables were included.

Generalized linear mixed models were fit using the Gamma family with a log link and estimated using maximum likelihood (ML) and the BOBYQA optimizer. Use of a log link produced log transformed estimates. Two separate models were included, (1) slope as the response variable and (2) asymptote as the response variable. To predict our response variables, LCC and WSP were included as fixed effects with interactions, and Ecoregion was included as a random effect impacting parameter slopes and intercepts.

Additionally, we employed the use of Generalized Additive Models to further assess non-linearities in the relationship between our response variables and predictors as they vary by Ecoregion. Predictor variables were converted to a logit binary classification with Forest = 1, Greenspace = 2, Agricultural = 3, and Urban = 4. This was also done with WSP where HW = 1, MW = 2 and OUT = 3. All analyses were completed in R version 4.0.3 (R Core Team, 2020) using the mgcv package (Wood, 2011).

Across the sites included in this study, a total of 216 species of birds were observed. This total richness represents approximately 80% of the roughly 270 confirmed breeding species of Ontario (Cadman et al., 2007). A total of 34,936-point count stations were included in our analysis. Ecoregion sampling varied with Lake Simcoe-Rideau being the most sampled region while Pigeon River was under sampled resulting in its exclusion from further analysis. The highest richness was observed in the Lake Simcoe-Rideau Agricultural landcover, the highest Shannon diversity was observed in Georgian Bay Greenspace landcover class, and the highest inverse Simpson was calculated in Georgian Bay Agricultural landcover (Supplementary Table S2). The values of these diversity indices are likely impacted by the levels of sampling, further supporting our decision to use SACs to equalize efforts among watershed position.

Across the Great Lakes Basin Watershed, slopes ranged from 1.47 to 3.14 whereas asymptotes ranged from 66.30 to 215.50 (Supplementary Tables S3, S4). The highest mean slope (2.30) was found in the Lake Abitibi – Forest – MW subgroup with the highest mean asymptote (162.47) found in the Georgian Bay – Forest – MW subgroup (Supplementary Tables S3, S4). The slope and asymptote parameters from the SACs are presented in boxplots (Figures 3, 4). Boxplots show variation in the rates of accrual due to a combination of landcover type, watershed position, and Ecoregion in which the points were found. Outliers were found at both the upper and lower ranges of the slope parameter with more outliers located at the upper range. In the asymptote data, outliers were only located at the upper range of the boxplots. This issue was addressed by iterating these curves 1,000 times and removing all asymptotes greater than 216 which was the total species count in this study from further analysis. The likelihood of overestimation rather than underestimation is a function of the species accrual methodology whereby it is not possible for the number of species to decrease with increasing sample size. With slopes, this is not the case. That is, it was possible to have a high slope in one iteration of the curve followed by a lower slope in another. A summary of slopes and asymptotes from the complete set of SACs can be found in supplementary tables and figures (Supplementary Tables S3, S4).

To gain a better understanding of the relationship between the asymptotes and slopes of our species accrual curves, we created scatter plots for each subgroup to examine these variables in concert. There was an inverse curvilinear relationship between slopes and their corresponding asymptotes such that steeper slopes corresponded with lower asymptotes (Figure 5). This inverse relationship was not entirely surprising given that with less heterogeneity in an LCC, species accrual is expected to be faster and corresponding asymptotes lower; whereas increased heterogeneity in an LCC should result in more unique species observed as sampled points increase, leading to a slower rate of accrual and higher asymptotes (Thompson and Withers, 2003). When accounting for watershed position, this relationship varied further with the Ecoregion – landcover combinations. Looking at Figure 5, we can see that even at the same slope, many HW, MW, and OUT accrual curves differed in their corresponding asymptote. This provides evidence that watershed position plays a role in the rate of species accrual.

The rate of species accrual tells the observer how rapidly a system increases in species with each additional point sampled (Colwell and Coddington, 1994; Chiarucci et al., 2008). Here, the rate of accrual provides an idea of the relative favourability of underlying habitat conditions, measured as the combination of Ecoregion, landcover type, and watershed position. One perspective about the information derived from a steep slope is that it is a system level metric where a smaller sample is required to reveal most of that system’s species richness (Thompson and Withers, 2003; Dove and Cribb, 2006). Alternatively, curves with shallower slopes, i.e., a slower rate of accrual, require a larger number of samples to achieve a similar level of richness (Thompson and Withers, 2003; Dove and Cribb, 2006). Although, not necessarily equal to a high level of richness, high rates of accrual may be favorable to a highly functional ecosystem particularly in remnant landscapes (Flather, 1996). Additional work is needed to further explore this idea.

Ecoregions are defined based on their climatic conditions, historical landcover, as well as surface water presence, all having the potential to affect the types of species present in that area (Crins et al., 2009; Yang et al., 2011). Based on our results, estimates for Lake Abitibi, Lake Nipigon, and Lake Temagami, which are located to the Northwest of the GLB showed faster rates of species accrual than Georgian Bay, Lake Erie – Lake Ontario and Lake Simcoe – Rideau, located in the more southern parts of the GLB (Figure 1 and Table 3). These patterns however were modified when looking at our subgroups in their entirety (Supplementary Tables S3 and Figure 3). Independently, our GAM results show that Lake Temagami had the highest parameter estimate for slope, and Georgian Bay exhibited the highest asymptote estimate (Table 3). Climatic variability in these Ecoregions (Table 1) plays a role in affecting rates of species accrual (Figure 3, Kallimanis et al., 2008). The Ecoregions of Lake Nipigon, Lake Abitibi, and Lake Temagami are all described as having cooler mean temperatures and shorter summers than Georgian Bay, Lake Erie - Lake Ontario and Lake Simcoe – Rideau. Historical land use patterns also varied among Ecoregions with higher levels of agriculture to the south and maintained forests to the northwest of the GLB (Crins et al., 2009; Figure 1). As expected, regional context and climatic conditions were captured by the Ecoregional (random) effects in our models. Surprisingly, Agriculture and Urban landcovers in Georgian Bay had the steepest slopes detected (after mostly northern forests) indicating relative uniformity in these landcover types because a smaller sample would be required to reveal most of the system’s richness. This result further highlights the effect of Ecoregion on species accrual rates. However, the same was not true for Agricultural and Urban landcover in southern regions (Simcoe, Erie-Ontario) where slopes were shallow indicating the need for more samples to determine maximum richness. Shallower slopes in more human managed landcover types are in keeping with the results of other studies (Flather, 1996). Several of these combinations showed that these factors further compound to affect the rates of accrual seen at a specific watershed position (Figure 3).

Our first research question focused on observing the ability of urban landscapes to retain a relative amount of species richness and rates of accrual. The ability of urban environments to maintain and attract ecological richness is of the utmost importance to the future of our cities. Our results show that urban environments have the potential to maintain a relatively high species richness that varies with watershed position (Figure 4). Several studies speak to the variability of species diversity patterns observed within urban settings (Ferenc et al., 2014; Hensley et al., 2019). This pattern is further enhanced when looking at urban areas in the context of the Ecoregions in which they are embedded. In our study, Ecoregions that were historically more urbanized, i.e., Lake Erie – Lake Ontario and Lake Simcoe – Rideau, have the potential to maintain species diversity (Figures 1, 4). Particularly for the Lake Simcoe – Rideau Ecoregion, Urban landcover classes at mid watershed positions exhibited a mean asymptote (159.32 – MW) and maximum species richness (138) that were among the highest overall (Figure 1; Supplementary Tables S2, S4). Much of the surrounding landcover associated with the urban environment in these regions is classified as agricultural, which in previous studies is noted to modify the species that occur within these landscapes (Sirami et al., 2019; Lee et al., 2022). In fact, Lake Erie – Lake Ontario and Lake Simcoe – Rideau are both described as having cropland (Agricultural) as their primary landcover type (Crins et al., 2009). This heterogeneous combination of landcover under the more southern climatic conditions of the Mixed Wood Plain Ecozone may provide ideal conditions for maintaining relatively high species richness (Figure 1 and Supplementary Table S2).

Greenspace estimates outperformed Forest in terms of rates of accrual and asymptotes (Table 2). Though this was unexpected, it may be explained by several factors. Our classification of Greenspace is comprised of different (12) landcover types that were characterized by some unique open spaces (e.g., Tallgrass Prairie, Swamp, and Marsh). Landcover heterogeneity is understood to influence species diversity metrics, and several of the landcover types included in our Greenspace index were found within our study area only in small amounts. All landcover types included in the Greenspace class are characterized by possessing <60% tree cover and our grouping was done to ensure sufficient bird sampling points were captured. Whereas the Forest cover class was less heterogeneous (Supplementary Table S1). Our Forest class is comprised mainly of coniferous, deciduous, and mixed forest. The variability of the cover classes included in Greenspace likely provides habitats that are favorable to a greater number of avian species than the habitat provided by the Forest landcover class due to the greater heterogeneity of the types included in the Greenspace class (Turner and Tjørve, 2005; Kallimanis et al., 2008; Stein et al., 2014).

We hypothesized that species accrual rates would be higher as you move toward downstream positions. From our results, when looking solely at the effect of watershed position, we see that MW and OUT exhibit estimates that are higher than those of HW (Table 2). However, in the context of our subgroups, the impact of watershed position was highly variable in relation to Ecoregion and the landcover types that these streams flowed through (Figure 3 and Figure 4). Indeed, like riparian areas, the ecological function and impact of all watershed positions is not the same (Baker et al., 2001; Thorp et al., 2006). The variability of SAC slopes and asymptote estimates with WSP are in-line with other papers that emphasize the effects of the variability of stream reaches and the additional impact of their surrounding landscape (Baker et al., 2001; Wiens, 2002; Yang et al., 2011; Soininen et al., 2015). Our results show that landcover type is further modified by watershed position in its effect on species diversity metrics. As our approach is based on observed patterns and statistical relationships, we cannot directly deduce causality; nor can we extract the exact mechanisms that cause these differences. However, there is an extensive literature that can be used to explore the mechanisms of WSP on species. The impact of watershed position can be explained best by understanding the potential effect of stream size on aquatic species accrual rates (Minshall et al., 1985). It can be hypothesized that the patterns that increase discharge and in-turn species richness, may also translate to impacts on neighboring terrestrial areas. In addition to the size of the streams, the structure of riparian areas and the impact that they have on the surrounding landscape also changes with a change in watershed position (Baker et al., 2001; Wiens, 2002; Allan, 2004). This relationship changes when looking at mean asymptotes with MW exhibiting the greatest impact on maximum richness (Table 2). Streams may also flow through several Ecoregions from headwater to outlet (longitudinal) resulting in further differences in hydrology, geomorphic features such as soil and climatic conditions including precipitation (Yang et al., 2011). There may also be differences in subsidies from freshwaters systems to the surrounding landscape at different locations along the stream (Lamberti et al., 2020). The variability of watershed position and its impact on species accrual rates remained throughout our study.

We face limitations however as we did not investigate SACs within specific watersheds, so we cannot speak to the impact of sampling the same watershed in different positions. This limitation was overcome partially by our use of randomizations within the 95% CIs of each slope or asymptote. We are confident that we captured the full range of variability in SACs across LCC and WSPs without explicitly accounting for within watershed effects.