Temporal change in urban fish biodiversity—Gains, losses, and drivers of change

Abstract Our aim was to examine temporal change in alpha and beta diversity of freshwater fish communities in rivers that have urbanized over the same period to understand the influence of changes in land use and river connectivity on community change. We used biological (2001–2018), land use (2000–2015), and connectivity data (1987–2017) from Toronto, Ontario, Canada. We used linear mixed effects models to determine the strength of upstream land use, connectivity, and their changes over time to explain temporal change in alpha and beta diversity indices. We examined beta diversity using the temporal beta diversity index (TBI) to assess site‐specific community change. The TBI was partitioned into gains and losses, and species‐specific changes in abundance were assessed using paired t‐tests. There were more gains than losses across the study sites as measured by TBI. We found little to no significant differences in species‐specific abundances at aggregated spatial scales (study region, watershed, stream order). We found different relationships between landscape and connectivity variables with the biodiversity indices tested; however, almost all estimated confidence intervals overlapped with zero and had low goodness‐of‐fit. More fish biodiversity gains than losses were found across the study region, as measured by TBI. We found TBI to be a useful indicator of change as it identifies key sites to further investigate. We found two high value TBI sites gained non‐native species, and one site shifted from a cool‐water to warm‐water species dominated community, both of which have management implications. Upstream catchment land use and connectivity had poor explanatory power for change in the measured biodiversity indices. Ultimately, such spatial–temporal datasets are invaluable and can reveal trends in biodiversity useful for environmental management when considering competing interests involved with urban sprawl in the ongoing “Decade on Restoration.”


| INTRODUC TI ON
Land-use change is one of the leading causes of global biodiversity change (Jung et al., 2019).Urbanization is one of the most common and dramatic forms of land-use change, and is characterized by increased human population density, increased artificial impermeable surface (e.g., roads, parking lots, and sidewalks), increased pollution, elevated temperature, and habitat fragmentation (Kondratyeva et al., 2020;Walsh et al., 2005).These changes can have profound impacts on all ecosystems experiencing urbanization.Land-use change can be used as a proxy for numerous finerscale aspects of urban change as several environmental variables that drive ecological change are often correlated with land-use change.For example, transitioning from forest to urban land cover is expected to impact hydrological dynamics and lead to more runoff, thus increasing the concentration of and types of contaminants flowing into water bodies (Arnold & Gibbons, 1996;Walsh et al., 2005).Potential ecological stressors include the concentration of the contaminant, and knowledge of both requires intensive sampling which is difficult to conduct at large spatial scales.
Calculating land-use change can be done consistently over large spatial extents through remote sensing and can reveal broader patterns of change across a watershed (Turner et al., 2007).The downside of using land-use change as a proxy for ecological change is the potential loss of fine-scale variation; however, the benefit of increased power to detect change and draw inferences over large spatial extents is appealing.
Currently, 56% of the global human population (83% in North America) lives in urban regions, and this is expected to grow to 68% by 2050 (89% in North America), making further urbanization imminent (Khor, 2022).Understanding how organisms respond to urban pressure is a fundamental ecological goal in the ongoing United Nations "Decade on Restoration" (2021)(2022)(2023)(2024)(2025)(2026)(2027)(2028)(2029)(2030).This is particularly crucial for freshwater organisms, which are estimated to have declined by 84% since the 1970s (WWF, 2020).The interaction between further urbanization and restoration goals is important to consider, especially due to the increased valuation of urban green and blue spaces, which are naturalized vegetated spaces such as parks and forests (green), and aquatic spaces like ponds, wetlands, and rivers (blue).The value of urban blue space was recently recognized within the 2023 Kunming-Montreal Global Biodiversity Framework (GBF), demonstrating the high level of global concern for the ecological integrity of urban water bodies (CBD, 2022).Target 12 of the GBF calls for more biodiversity-inclusive planning and greater access to higherquality blue and green spaces in urban regions (CBD, 2022).Target 19 calls for more effective resource mobilization, including domestic resource mobilization for conservation, and Target 21 reiterates the importance of biodiversity monitoring (CBD, 2022).Regional specific biodiversity monitoring data provide opportunities to study how urbanization causes ecological changes, which in turn informs how resources can be effectively mobilized for conservation research and better biodiversity-inclusive urban planning into the future.We leveraged regional monitoring data to further the understanding of how freshwater ecological communities have changed in urban regions.We analyzed trends in alpha and beta diversity in a large urban region using long-term fish community monitoring data from a local watershed management agency.Fish species distributions are closely tied to water quality and stream condition and thus change in fish communities can be indicative of broader environmental changes (Trautwein et al., 2012).Further, long-term datasets collected through standardized sampling methodology provide invaluable opportunities to study ecological change (Hughes et al., 2017).The standardized methodology of long-term monitoring datasets provides critical data which can be compared across time facilitating confident conclusions on biodiversity change (Turner et al., 2015).Long-term ecological datasets allow for comparisons of within-site change, instead of relying on space-for-time substitutions which are common in urban ecological studies when long-term datasets are unavailable (Damgaard, 2019).Temporal datasets better control for aspects of site uniqueness, as many urban sites are impacted differently through time due to historic patterns of landuse (i.e., legacy effects), which may lead to artifacts in space-for-time substitution approaches (Foster et al., 2003).
Using several measures of biodiversity change better encapsulates ecological change as each provides a different interpretation or estimate of the same underlying data patterns (Lyashevska & Farnsworth, 2012).Alpha diversity measures community diversity at a particular site, with different metrics placing emphasis on the number of species or the abundance of individual species (Jost, 2007).
To capture community change, measures of alpha diversity can be compared among sites or at the same site through time.However, when alpha diversity is considered alone, information on species identity can be lost.Beta diversity measures community variation among sites (Koleff et al., 2003), or among time periods at the same site, with comparisons across space being more common, in part due to a lack of long-term within-site data.Importantly, measures of beta diversity can retain information on species identity, providing information on the loss and gain of individual species at sites (Baselga, 2008(Baselga, , 2010)).Long-term datasets with repeated measures at multiple sites allow for comparisons across both space and time (Hughes et al., 2017).Spatial-temporal datasets like the one we use are becoming more common as long-term local monitoring datasets mature.Used in tandem, alpha and beta diversity indices are useful to understand broad patterns of biodiversity change through both number of species and composition of communities (Magurran, 2021).
Throughout history, urban centers have been tied to water (Hein, 2016).Living near freshwater provides access to drinking water, irrigation sources, transportation corridors, and cultural practices (Hein, 2016).Expansion of cities to accommodate increased human population growth, and rapid increases in migration from rural to urban regions in recent decades, has led to increased density of humans on the shores and in the watersheds of freshwater systems (Haidvogl, 2018).As a result, freshwater rivers and lakes have become degraded over time through multiple mechanisms including fragmentation due to small dams and mills, anthropogenic pollution (e.g., sewage, industrial waste, and household waste), and increased prevalence of introduced nonnative species (Haidvogl, 2018).The removal of soils and natural vegetation and increases in impermeable surfaces have fundamentally reduced, if not eliminated, the hydrological buffering capacity of urban watersheds leading to flashy freshwater systems prone to flooding (Arnold & Gibbons, 1996).Changes in hydrological dynamics, species dispersal ability due to in-stream barriers, and increased pollutant loading can inevitably alter ecological dynamics in freshwater bodies, often causing declines in species sensitive to environmental change (Walsh et al., 2005).
Hydrological connectivity is important for riverine dynamics as changes in connectivity impact not only species dispersal patterns but also nutrient cycling and productivity (Ward & Stanford, 1983).
Reduced ecological connectivity due to historic mills and dams, and contemporary roadway culverts, buried stream sections, and dams impacts a variety of ecological processes such as access to spawning grounds, access to prey, and ability to avoid predators (Choy et al., 2018;Edge et al., 2017).At the metacommunity level, decreased connectivity can reduce species colonization of sites from larger source populations which may supplement smaller populations across a watershed (Brown & Swan, 2010).Therefore, isolated sites may see higher rates of local extinction due to the lack of immigration from other sites (i.e., less dispersal).
To understand how freshwater biodiversity is changing in a highly urbanized region, we examined changes in riverine fish communities and potential drivers of change across 15 years in Toronto, Ontario, Canada's largest urban center.Toronto has experienced increased urban sprawl since the beginning of European settlements in the region in the mid-18th century (Careless, 1984).
Urbanization progressed steadily from the northern shore of Lake Ontario northward throughout the region's river valleys (Careless, 1984).Through time, the City buried rivers in underground pipes and dredged or filled wetlands to reduce exposure to polluted water and increase developable lands (Bonnell, 2014).
Three large river systems and several smaller systems are found in Toronto, of which each has varying degrees of changes in their land-use history due to uneven patterns of urban sprawl through time.In our study, we compared alpha and beta diversity of fish communities through time and examined drivers of change through analyzing land-use change due to wetland loss, woodland loss, and anthropogenic land intensification over the study time period.Additionally, we examined the impact of connectivity through time by considering dendritic connectivity and change in connectivity as a potential driver of fish community change.The first objective of our study was to examine whether fish communities are experiencing more gains or losses in species or abundance-per-species through time (Question 1).We predicted fish communities would show more losses than gains due to increased urbanization in the study region over time and changing environmental quality ultimately contributing to species decline and local extirpation.Our second objective was to determine whether changes in upstream catchment land use and dendritic connectivity explained variation in fish community indices through time (Question 2).We predicted land-use change, specifically loss of more natural land-use cover (woodland, and wetland), would best explain variation in withinsite biodiversity through time.

| Fish community data sampling
Fish community data were obtained from the Toronto and Region Conservation Authority through their Regional Watershed Monitoring Program (TRCA, 2022).Sites in the same watershed were sampled once every 3 years, and watersheds were sampled in different years.Fish sampling was performed using the Ontario Stream Assessment Protocol and single-pass electrofishing (Jarvie & Jackson, 2022;Stanfield, 2010).Fifty-four sampling locations from five watersheds were included in our analysis: Etobicoke Creek (n = 11), East Don River (n = 13), Highland Creek (n = 5), Mimico Creek (n = 3), and the Rouge River (n = 22) (Figure 1).Individuals were identified to species, and abundances were tallied by the TRCA.Where samples were not identified to species level, species were assigned based on site history and/or expert knowledge or were removed from the dataset (see Table S1 for details).Fish data were summarized at the site level for two different time periods, T1: 2001-2003 (T1) and T2: 2016-2018 (T2).Each time period included one sampling event per site.Site Strahler order ranged from 1 to 6.

| Ecological diversity indices
Species occurrence and abundance data were used to compute three alpha diversity metrics at T1 and T2: species richness, Pielou's evenness, and Shannon diversity.Richness measures the number of species occurring within a site.Pielou's evenness measures the relative abundances of species and is bounded between 0 and 1, with 1 indicating a completely even community (Jost, 2010).Pielou's evenness was calculated through the equation J = H′/log(S) where J is Pielou's evenness, S is richness, and H′ is the Shannon diversity index.The Shannon diversity index accounts for both species richness and relative abundance, and is considered to be more influenced by richness than Pielou's evenness in relatively small communities (Dejong, 1975).Each metric is commonly used to summarize biodiversity, and there are benefits and downsides to each (Magurran, 2012).Using all three metrics simultaneously better encapsulates the structure and change of biological communities.
Alpha diversity metrics were calculated with the vegan package (Oksanen et al., 2020) using the specnumber() function for Pielou's evenness and diversity() for Shannon's diversity.
Change in fish community composition over time was determined using the temporal beta diversity index (TBI) (Legendre, 2019).
To calculate the TBI, a relevant dissimilarity index is chosen, and a dissimilarity matrix is created.We use percent difference (i.e., Bray-Curtis dissimilarity), which has been shown to be well suited to analyze community data (Legendre, 2019).When using percent difference dissimilarity, the direction of community change within a site can be calculated through D%diff = (B + C)/(2A + B + C), where D%diff is percentage difference, B is where abundance is higher in T1 than T2 (losses), C is where abundance is higher in T2 than T1 (gains), and A is common abundance between T1 and T2 (similarity) (Legendre, 2019).Larger values of the TBI (D%diff) indicate more change between T1 and T2.Significance of change within a site between time periods in relation to other sites can be tested using a permutation approach as outlined by Legendre (2019).This approach uses a one-tailed test of significance toward the upper tail; thus, only sites with exceptionally large change will be found to be "significantly" different (Legendre, 2019).Permutation is then followed by multiple test correction, wherein we used a Holm correction for multiple testing (Legendre, 2019).Exceptional sites can then be analyzed further to understand what underlying processes may be driving potential change, we consider both pre and post multiple testing correction results (see Legendre (2019) for details).The TBI was calculated using the TBI() function within the adespatial package of R using the method parameter set to percent difference (%difference) and 9999 permutations (Dray et al., 2021;Legendre, 2019).

| Land cover
Land-cover data were derived from the Southern Ontario Land Resource System, created by the Ontario Ministry of Natural Resources and Forestry.Land-cover categories for T1 were derived from SOLRIS v1.2 by condensing definitions into three categories: anthropogenic, woodland, and wetland (Table S2).Land-cover change between T1 and T2 was derived from SOLRIS v3.0 (OMNRF, 2019) by condensing the change categories into three categories: anthropogenic intensification, woodland loss, and wetland loss (Table S3).
Both datasets are to a resolution of 0.5 hectares (OMNRF, 2019).data were summarized at the entire upstream catchment level for each site through ArcGIS hydrotools (Figures S1 and S2).ArcGIS Pro Version 3.0.3was used for all land cover analyses.

| Connectivity
Dendritic connectivity was estimated using the Dendritic Connectivity Index (DCI) (Cote et al., 2009).For each stream segment, we calculated DCIs, a measure of how well an individual stream segment is connected to all other stream segments in 1987 (T1) (Choy et al., 2018) and 2017 (T2) (Edge et al., 2017).The DCI at T2 was calculated by recording the location and measuring down-

| Explanatory modeling
Linear mixed effects models were used to determine which variables best explained the change in alpha and beta diversity indices between T1 and T2.A Gaussian distribution was used for all models as the response variables, which are measures of change, were normally distributed.Candidate models included predictor variables of land-cover categories at T1, land-cover change categories between T1 and T2, dendritic connectivity at T1, and dendritic connectivity change between T1 and T2 (Table 1).The TBI compares two time periods; therefore the TBI itself is used as a temporal response variable (Legendre, 2019).Watershed was retained as a random effect for each response variable.Differences in diversity patterns between watersheds are possible due to varying spatial patterns of urbanization and restoration through time creating legacy effects (Foster et al., 2003).
The candidate models were compared using an AIC modelselection framework with a correction for small sample size (AICc) (Harrison et al., 2018;Symonds & Moussalli, 2011).Models were considered top models in our model-selection framework if their ΔAIC < 2 (Burnham et al., 2011;Harrison et al., 2018).We tested models containing interactions between each variable and its change component; however, models with interactions had high variance inflation factors so the interactions were not included in the final candidate models (Table 1).All predictor variables were standardized to z-scores prior to analyses.Model averaging was performed for all models with ΔAIC < 2 for each response variable, and subset estimates are reported.Model selection through AICc was performed using the function aictab() through the package AICcmodavg (Mazerolle, 2020) in R with models with ΔAIC < 2 retained for model averaging.Model averaging was performed using the function model.avg()from the package MuMIn (Barton, 2022;Harrison et al., 2018).R 2 for mixed models were calculated using the function r.squaredGLMM from the MuMIn R package (Barton, 2022;Harrison et al., 2018).R version 4.2.2 was used for all statistical analyses and visualizations.

| Individual species abundance
Changes in individual species abundance across the study region, each watershed, and each stream order were tested individually using paired t-tests with Holm corrections for multiple testing (Holm, 1979;Legendre & Condit, 2019).This approach allows for the testing of trends in significant change in abundance of individual species through aggregating sites to a higher level of organization (e.g., study, watershed, and stream order).Paired t-tests to analyze individual species change were computed with 9999 random permutations followed by a Holm correction for multiple testing through the adespatial package and functions paired.

| Data summary
A total of 42 fish species were observed within the compilation of Changes in site species richness between the two time periods ranged between −7 and 6, (i.e., sites lost between 0 and 7 species, and gained between 0 and 6 species; Figure S3).Changes in Pielou's evenness between T1 and T2 ranged between −0.76 and 0.55 (Figure S3).Changes in Shannon diversity ranged between −0.76 and 1.57 (Figure S3).Across all sites, the TBI ranged between 0.26 and 0.96 (Figure S3).

| Temporal beta diversity index
Across all sites sampled, there were more gains than losses for species or abundances-per-species average.Averaged gains (0.59) were larger than losses (0.41) (i.e., species gains averaged 59% and losses averaged 41%; Figure 2).We found seven sites with significant change in species or abundances-per-species (p < .05) prior to Holm correction (Figure 2).These included two sites on the Don River (DN011WM, DN010WMb), one site on the Rouge River (RG019WM), one site on Highland Creek (HL005WM), one site on Mimico Creek (MM002WM), and two sites on Etobicoke Creek (EC001WM, EC010WM) (Figure 3; Table S4).When the data were re-analyzed using occurrence rather than abundance data, two sites were found to have significant change with one gain site (EC010WM) and one loss site (EC001WM) prior to Holm correction (Table S5).We present p-values both pre-and post-Holm correction for multiple testing, as both values may be of interest to ecologists.The TBI values themselves can be used in modeling without significance testing (see Kuczynski et al., 2018;Robichaud & Rooney, 2022); however, we chose to include significance testing for consideration.The 7 exceptional sites pre-correction are subsequently referred to as high TBI sites, as the pre-correction significant sites represent the seven highest TBI sites, all reporting a TBI value >0.89.These seven sites represent the sites with the largest compositional differences through time (Legendre & Condit, 2019).

| Individual species change
Paired t-tests were performed across all sites, within each watershed, and within each stream order with a correction for multiple testing.When corrected for multiple testing, there were no significant differences in species abundances between the two time periods when sites were examined at the study level, watershed level, or stream order level (Tables S6-S8) other than a slight mean decrease in white sucker (Catostomus commersonii) in the Don River watershed (Table S7).
F I G U R E 2 B-C plot for the temporal beta diversity index (TBI).Squares represent sites where gains > losses.
Circles represent sites where gains < losses.Colors denote watersheds.The location of the centroid line is above the 1-1 line, suggesting the system reports more gains overall than losses.Asterisks represent high TBI sites (prior to holm correction p < .05sites).Size of points represents relative TBI values.
For the high TBI sites, four had more gains than losses, and three had more losses than gains (Figures 2 and 3).There were no significant differences in species abundances between the two time periods when sites were compared across all high TBI sites (Table S9), to all high TBI sites with overall gains (Table S10), or all high TBI sites with overall losses (Table S11).A closer examination of species-specific abundance changes across all the sampling sites revealed an increase in the invasive round goby (Neogobius melanostomus) from 0 to 78 with 74 sampled at one site (DN001WM), and an increase from 0 to 4 sampled at one high TBI site (EC001WM).Additionally, introduced green sunfish (Lepomis cyanellus) increased in two high TBI sites (EC001WM and EC010WM) to abundances of 1 and 255, respectively.

| Model results
The two averaged top models (Table 2) for richness indicated a positive relationship between change in richness and wetland at T1, wetland loss, and DCI at T1, although estimate 95% confidence intervals overlapped with zero except for wetland loss (Figure 4).There was a negative relationship between change in richness and DCI change, although the estimate 95% confidence intervals overlapped with zero (Figure 4).
The four averaged top models (Table 2) for Shannon diversity indicated positive relationship between change in Shannon diversity and woodland T1, wetland T1, wetland loss, woodland loss, anthropogenic T1, anthropogenic intensification, and DCI T1, although estimate 95% confidence intervals overlapped with zero (Figure 4).
There was a slight negative relationship between change in Shannon diversity and DCI change, although estimate 95% confidence intervals overlapped with 0 (Figure 4).
The three averaged top models for Pielou's evenness indicate a positive relationship between change in Pielou's evenness and woodland T1, wetland T1, anthropogenic T1, although estimate confidence intervals overlapped with zero (Figure 4).There was a negative relationship with woodland loss, wetland loss, and anthropogenic intensification, although all confidence intervals overlapped with zero (Figure 4).
The three averaged top models (Table 2) for TBI indicate a negative relationship between TBI and woodland T1, woodland loss, wetland T1, and anthropogenic T1, although estimate confidence intervals overlapped with zero (Figure 4).There was a positive relationship between TBI and wetland loss and anthropogenic intensification, although estimate confidence intervals overlapped with zero (Figure 4).

| DISCUSS ION
Our TBI analysis revealed a set of seven sites with substantial change in community composition over time.Further exploratory analysis using occurrence data rather than abundance data revealed a drop in the number of sites signaling significant change pre-Holm correction for multiple testing.This result may indicate species have already been lost prior to the fish surveys in 2001, making fluctuations in abundance of remaining species more common and influential in our TBI analysis.The different results in using abundance over occurrence data also indicate some sites may show signs of local extinction debt, where poor conditions reduce abundances, but the species continues to persist at the site for the time being (Figueiredo et al., 2019).The TBI is a useful metric for identifying sites with large community change and can subsequently be used in further modeling to understand differences in drivers between gains and losses of species of abundance-per-species (Legendre, 2019;Legendre & Condit, 2019).Further analysis of abundance trends related specifically to gains and losses may reveal underlying community dynamics related to extinction and colonization patterns (Carvalho et al., 2012).
The type of change differed between the high TBI sites, with three sites showing losses and four sites showing gains overall, one high TBI site gaining the invasive round goby, and two high TBI sites gaining green sunfish.The observed 4:3 ratio of gains to losses is the opposite direction of our prediction of more losses than gains.
The TBI is determined by change in species presence and/or species abundance; we used abundance-per-species due to availability of count data.species (Padayachee et al., 2017); thus, we may expect TBI to increase through time if invasive species presence and abundance increase, which can have an immediate impact on the interpretation of biodiversity metrics, and a longer-term impact on local biodiversity through changing community dynamics (Gaertner et al., 2017).For example, round goby invasion of tributaries in the Great Lakes region is of concern, as round goby have been found to impact native species through competition for resources and predation on early life stages of native species (Poos et al., 2010).individual exceptional change sites within a broader urban matrix which may become of management interest.
The TBI can be scaled up to a higher level of biological organization than species (Legendre, 2019).While the focus of the current paper was species-level biodiversity change, higher levels of biological organization such as family, or other classifications based on thermal guild or functional traits can be used in future studies (Austen et al., 1994;Benoit et al., 2021).Guild level analyses may reveal broader patterns of change which can be obscured by species-level investigation.Our results suggest thermal guilds may be an important categorization to consider, as we found one high TBI site shifted from a cooler-water fish dominated community to a warm-water fish dominated community (EC010WM).It is possible land-use change upstream of this site contributed to changes in water temperature, as the site saw anthropogenic land intensification directly upstream, along with small pockets of woodland and wetland loss (Figure 3).
While our temporal beta diversity analysis revealed gains and losses in fish biodiversity across the study sites, our chosen potential drivers of fish community gains and losses had low explanatory power.Model-averaged predictor-estimate confidence intervals overlapped with zero for almost all estimates and reported relatively low R 2 values (Figure 4; Table 2).to analyze pre-existing long-term datasets must be a regional and national priority to meet international biodiversity goals.
Percentage land-cover change per land-use category was measured as the change in land cover per category divided by the watershed area per site.Anthropogenic intensification represents anthropogenic lands which were further urbanized, for example, if an area changed from permeable urban to impermeable urban.Land-cover F I G U R E 1 Map of study sites across Etobicoke Creek, Mimico Creek, Highland Creek, the Don River, and the Rouge River.Black points represent study sites.Blue lines represent river or creek watercourse.
stream pool depth and barrier height of all existing barriers.The DCIs at T1 were calculated by compiling a list of mitigated or removed barriers in the streams and adding those barriers to the surveyed list at T2.Mitigated or removed barriers were identified by contacting groups and agencies (e.g., Trout Unlimited, Ontario Streams, Ontario Ministry of Natural Resources, and Department of Fisheries and Oceans) that have been involved with barrier mitigation projects in the Toronto Region, reviewing permit applications under Ontario Regulation 166/06 in the Fisheries and Oceans Canada database, and reviewing the corporate records database of the Toronto and Region Conservation Authority.Additional details of the methodology to identify removed or mitigated barriers can be found in Choy et al. (2018).The permeability of each barrier was estimated using outlet drop and baseflow of the stream(Choy et al., 2018;Edge et al., 2017).
sites and time periods; the range in species richness observed differed among the sites within each watershed: Rouge River (2

F
Map of sites with the temporal beta diversity index (TBI) results.Size represents relative TBI value.Plus sign sites have more gains than losses; circle points represent more losses than gains.Red sites represent high TBI sites (pre-Holm correction significant sites).Gray polygons represent urban intensification.Green polygons represent woodland loss.Purple polygons represent wetland loss.Site identification numbers increase in the upstream direction.
Such competition and predation by non-native round goby may further impact already declining species by further reducing the abundance of native species, a theory supported by the results ofKindree et al. (2023).Comparing change in individual species abundances across the study region, each individual watershed, and each stream order revealed little significant change in individual species abundances when corrected for multiple comparison testing.This result may be due to diverse spatial patterns of urbanization across the region causing site-specific legacy effects(Perring et al.,   Wohl, 2019), and the inherent dendritic nature of rivers creating heterogeneous habitat conditions between sites(Benda et al., 2004).Spatial patterns of urbanization may differ among sites due to time since land-use change, proximity to direct land-use change, or type of land-use change occurring within a watershed.For example, some sites further downstream see little change in upstream catchment land use as they have been urbanized for over 100 years, whereas some sites closer to headwaters see higher relative levels of land-use change as headwaters are relatively less urbanized areas than downstream sites in the study region (Figure S2).Such site specificity may obscure site-level change when aggregated to a broader spatial extent.While inconvenient for making broaderscale conclusions on drivers of change in species abundances, this difference suggests metrics like the TBI are useful for identifying F I G U R E 4 Model-averaged estimates for richness, Shannon diversity, Pielou's evenness, and temporal beta diversity index (TBI) with 95% confidence intervals.The TBI intercept value of 0.62 with 95% CI [0.57, 0.67] not shown due to data visualization limitations.All estimates represent subset estimates.
While calculating landuse change can be done consistently over large spatial extents through remote sensing, fine-scale variation can be obscured since land use is a proxy encapsulating many urban changes, and this may have contributed to our results.Finer geographic scale like local upstream catchment at smaller distance buffers or riparian zone land use may reveal stronger relationships to biodiversity indices; however, we were interested in cumulative upstream land-use impacts at the entire upstream catchment scale.Degree of land-use change may have fallen below thresholds previously found to explain ecological change.For example, while thresholds for land-use change tend to be geographically specific, they were found to range from 1% to 12% for urbanization and 2% to 37% for agriculture in a recent global study(Chen & Olden, 2020).The effects of land use may also be statistically stronger when considering factors such as habitat composition.While data to this scale are often not available due to collection limitations through time, it is possible finer-scale changes in habitat may be stronger drivers of biodiversity change(Chase et al., 2018).Additionally, it is possible physical variables like change in water temperature, may explain some variation(Jarvie & Jackson, 2022), as one site of interest appears to show a change in thermal community composition.As suggested byKuczynski et al. (2018), including predictors such as specific pollutants or flow may also improve explanation of variation within sites through time.As temporal datasets mature and ecological monitoring continues, analyzing trends over time periods longer than the 15 years analyzed in the current study may reveal stronger relationships between drivers of biodiversity change if time lags are present(Watts et al., 2020).Ultimately, while curating temporal biotic datasets is of the utmost importance, collecting physical and chemical datasets alongside biological data will aid in future explanatory modeling efforts.The value of open-source data gathered through standardized sampling methods in long-term monitoring programs should not be underestimated.Such long-term datasets provide the basis to answer numerous ecological and conservation-oriented questions.Long-term datasets are particularly valuable in regions experiencing ongoing land-use change and because we cannot logistically (nor ethically) plan and implement the corresponding large-scale experimental studies that would allow testing of these ecological responses (i.e., experiments to test the effects of environmental stressors common to urbanized regions at the scale of urban centers).The evolution of the "sustainable city" from the "industrial" and "sanitary" city makes for an interesting comparison of ecological change over time(Kaushal et al., 2015).The approach in this paper can be used to examine ecological change, whether due to further urbanization or restoration.Future studies can include more land-use categories where land is restored to more natural land cover, unlike the current study where only land use change in the direction of further anthropogenic intensification was included.However, such studies will only be possible with continued monitoring of ecosystems through time using standardized approaches such as those used in long-term monitoring programs.Considering the ongoing "Decade on Restoration" and the recent 2023 Kunming-Montreal GBF, the collection of data to inform biodiversity-inclusive planning and the funding of research While sites indicate overall gains or losses, investigating individual species change is important to identify whether gains are gains in species or abundance, and whether gains are native or non-native species.Gains in native species may be reflective of recovery due to restoration efforts but gains in non-native species be indicative of broader invasion dynamics ongoing in urban regions.Model selection results for richness, Shannon diversity, Pielou's evenness, and the temporal beta diversity index (TBI).Anthro.represents anthropogenic model.Models with <2 delta AIC include R 2 values.Models <2 delta AIC were included in model averaging and data visualization.