The study was conducted in the Monjolinho River Basin, located in the State of São Paulo, south-eastern Brazil (Fig. 1). The regional climate is subtropical humid mesothermal or CWA according to the Köppen classification. The Basin has an area of approximately 275 km² and includes the entire urban zone of São Carlos Municipality. The main river, Monjolinho, has a length of approximately 43.25 km, originates in the east of the city at an elevation of about 900 m and flows in the east-west direction (Espíndola 2000). The administrative area of São Carlos covers 1137.303 km² with a population of 221,950 inhabitants, 96% of which is urban (IBGE 2010). The land cover is composed of approximately 60% agricultural areas, where the most common activity is sugarcane cultivation; the remainder is almost equally divided between natural vegetation and urban areas (Ferreira and Cunha-Santinho 2014).

Figure 1.

Characterisation of the study area in the Monjolinho River Basin, central region of the State of São Paulo, south-eastern Brazil.

The population of São Carlos has grown significantly in the last decades (Costa et al. 2013) together with intensive agricultural activities which have been increasingly impacting the local native vegetation and natural resources (Costa et al. 2013, Ferreira and Cunha-Santinho 2014, Neves et al. 2018). The urbanisation has also created a heat island effect which is worsened by air pollution due to the use of fossil fuels. Besides that, the annual dry period (about 6 months), a characteristic of the climatic region, contributes to this effect (Costa et al. 2013, CETESB 2017, Neves et al. 2018).

The original vegetation in the Monjolinho River Basin was characterised by the presence of Cerradão (forest), open and dense savannah phytophysiognomies, mesophytic forest (semi-deciduous seasonal forests) and other vegetation types, such as riparian forests and open wetlands (Soares et al. 2003).

In this study, we carried out a comprehensive assessment of biodiversity and ecosystem services, focusing on key indicators: nature-based recreation, cooling effect (microclimate regulation), air pollutant removal (measured using PM₁₀) and bird communities (measured using bird richness).

We selected 43 sample points in the watershed that are representative of all ecosystems and landscape types, including transition areas. The points have a minimum distance of 200 m from each other and were located along watercourses (Fig. 1). The sampling points were randomly selected, taking into account factors such as minimum distance, accessibility and the effort to encompass all types of riparian environments within the watershed, ensuring a comprehensive assessment of biodiversity. In a subsequent step, we evaluated bird richness (BR), land use and land cover (LULC) and three Ecosystem Services (ES) employing the methodologies described below.

We assessed the land use and land cover and the canopy cover of urban forest inside a buffer of 500 m from each survey point. We used Google images from Open Layers plugin in the Quantum GIS (QGIS) software and data from OpenStreetMap (OSM) and then classified according to the CORINE legend with some adaptation for local features (Suppl. materials 1, 2 - https://land.copernicus.eu/user-corner/technical-library/corine-land-cover-nomenclature-guidelines/html). The land use was classified by using photo-interpretation of remotely-sensed data and field visits for ground-truthing for the year 2018.

For each model, we generated a final raster image, from which we extracted the mean values of ES, using the zonal statistics of QGIS algorithms on the processing toolbox, for each 500 m buffer of the 43 surveyed points. Based on the results, we assigned a score to each sample area. First, we calculated the mean value of each ES in the 500-m buffer around each sample point and bird richness and then normalised the results between 0-1, where 1 was the highest value amongst the sample areas. After that, we classified the results into three classes, dividing the range 0-1 into three equal parts (low – 0-0.33, medium – 0.33-0.66, high – 0.66-1.00) to facilitate the visualisation and comparison amongst the results.

To assess synergies and trade-offs between all the ES and biodiversity indicators, we ran the correlation bivariate model using the Reduced Major Axis (RMA) algorithm with log transformation to data aiming to avoid errors related to different measures used to calculate each indicator (Hammer et al. 2001). To investigate the variation across the urban-rural-natural gradient, we classified and grouped points according to five major landscape types defined by the land-cover proportion and dominance of the following landscape units: native vegetation and superficial water, crops and pastures and urbanised areas. Thus, we categorised the areas surrounding the sampling points in five landscape types:

1. urban (URB);
2. peri-urban (PERI), i.e. the transition between urban areas and non-urban areas;
3. farmland (FARM);
4. transition (TRAN), defined as the transitional landscapes between farmland and native vegetation; and
5. nature (NAT), comprising native vegetation and superficial water (streams, river, lakes) (Suppl. material 1).

We tested the influence of these five landscape types on the mean values of ES performing the analysis of variance (ANOVA) and the post-hoc Tukey's Q test to identify significant differences between pairs of landscape types, using PAST Program (Hammer et al. 2001), v.4.03 released in 2020. For ES datasets with at least one out of five matrix types showing non-normal distribution confirmed with the Shapiro-Wilk normality test, we proceeded with a log transformation. We considered the significance level alpha = 0.05 to accept an alternative hypothesis of each test.

We identified available public areas which could be targeted by a reforestation programme (Fig. 2) to simulate the future scenarios with an improvement in the canopy cover of public green areas. As available areas, we used the difference between areas covered by trees and not covered on the selected public green areas. Inside the areas available for future reforestation, we simulated possible scenarios considering the use of 25%, 50%, 75% and 100% for the improvement of canopy cover. Then we recalculated the cooling effect and air pollutant removal models for each of these four percentages of possible scenarios. To run the air pollutant removal, we adopted a LAI index of 0.9 m²/m² on the newly-planted areas (this value was selected through a good and tangible value, based on the values of existing canopies in good conditions).

Figure 2.

Division of micro-basins of the Monjolinho River (main basin) in the urban area, with indication of the buffers where the future possible scenarios of reforestation were estimated and the planting sites (yellow) identified for the scenarios.

To better address the urban management and facilitate the decision-making, we assessed the results of the future scenarios for each urban watershed area (Fig. 2). The micro-basins were delimited in QGIS using GRASS 7 (r.watershed model) in the processing toolbox with a 30 m resolution Digital Elevation Model (DEM) (images 21S48 and 22S48) elaborated from SRTM data from USGS, available in the geomorphological database of the Brazilian National Institute for Space Research (INPE) (http://www.webmapit.com.br/inpe/topodata/).

With the increase of forested areas, we assumed that this scenario could influence the bird richness and the available services offered by these animals. Thus, we tried to predict the consequences of a management focused on reforestation by identifying the relationship between the vegetation cover area and the number of bird species. We used the polynomial regression analysis, based on a least-squares criterion and singular value decomposition, to select the best-fit model using the Akaike Information Criterion – AIC value. For the selected model, we tested the significance of the fit using the F test and obtained the coefficient of determination (R²), i.e. the proportion of the variance explained by the model (Hammer 2020). This procedure was implemented in PAST Program (Hammer et al. 2001) using data collected in the selected micro-basins inside the urban zone, which represent the current land use and added data collected in the peri-urban zone, representing a scenario with more vegetation and ecological corridors, but with a similar species pool of the bird community composed by those that can inhabit completely urbanised habitats.

The indicators of Ecosystem Services (ES) assessed along watercourses in the Monjolinho Basin were categorised into three intensity classes (low, medium, high) for each sampling point, based on the normalised values (detailed data in Suppl. material 2). These data were plotted on the maps showing the spatial patterns of ES provisioning intensity for pollutant removal and cooling effect, the availability of recreation opportunities for people and the bird richness (number of bird species sampled in situ) occurring in each sampling point (Fig. 3).

Figure 3.

Supply of Ecosystem Services in sampling points along waterbodies of Monjolinho Basin, south-eastern Brazil.

The different zones of land use at landscape scale influenced the mean values of ES provided by nature: pollutant removal (F = 23.64, p < 0.0001), cooling effect (F = 25.11, p < 0.001) and bird richness (F = 11.29, p < 0.001) as a surrogate for services offered by animals; and the nature-based recreation opportunities for people (F = 20.19, p < 0.001). The variation of ES values of each zone is detailed in Fig. 4. The post-hoc test indicated that the points completely immersed in the urban matrix were significantly different from other landscape zones, considering paired comparisons for each of four ES (p < 0.05, see Suppl. material 3).

Figure 4.

Variation in Ecological Services values across different landscape zones along waterbodies in the main basin of Monjolinho River, south-eastern Brazil and respective ANOVA results.

While the urban landscapes showed the lowest means of pollutant removal, cooling effect and bird richness, those areas have the highest means of recreation opportunities (Fig. 4). Additionally, we found significant differences between points located in landscapes with a higher proportion of native vegetation and the points in peri-urban areas, regarding cooling effect and recreation (p < 0.01). Additionally, areas with native vegetation are significantly different from all other zones in terms of pollutant removal (p < 0.012).

The assessment of relationships between each ES (Fig. 5) indicated a significant positive correlation between PM₁₀ removal and cooling effect (r = 0.827, R² = 0.684, p < 0.001). Towards the same direction, the bird richness-ES correlations pointed to a positive and significant relationship with PM₁₀ removal (r = 0.539, R² = 0.29, p < 0.001) and cooling effect (r = 0.649, R² = 0.42, p < 0.001). Beside these positive relationships, we found trade-offs with the recreation opportunities, that exhibited a significant negative correlation with cooling effect (r = -0.578, R² = 0.335, p = 0.048), with PM₁₀ removal (r = -0.419, R² = 0.176, p = 0.005) and with bird richness (r = -0.59, R² = 0.348, p < 0.001).

Figure 5.

Linear models indicating the correlation between Ecosystem Services. The diagonal above shows the surveyed points representing the five landscape matrices plotted in linear correlation models. The diagonal below shows the correlation results. Dark and light blue colours: positive correlations. Red and orange colours: negative correlation.

Four possible scenarios were simulated to understand how the increase in forest cover area could enhance the two ES which are primarily determined by the presence of vegetation, the cooling effect and pollutant removal, as shown in Table 1. The available sites to allocate restoration efforts correspond to only 1.46% of the total area. However, in the best scenario with the restoration actions applied in 100% of the available area, the pollutant removal indicator can be increased from 2.61% up to 7.73% and cooling effect from 1.35% up to 3.34%, according to the sections (micro-basins) of the main river basin which varied in available areas to be restored. There is one micro-basin without available area to increase these two ES, but it showed the highest current values. In another case, the available area of one micro-basin is so small (0.973 ha) that no changes in the cooling effect indicator were observed in the four possible future scenarios (Table 1).

Bird species richness, which represents a set of ES being offered in urban areas, can be increased with the addition of forested areas, based on the model that considers data collected in points of urban and peri-urban zones (Fig. 5). Sampled points immersed in the urban matrix are the focus for implementation of nature-based solutions and showed the lowest mean number of bird species (Fig. 3). Since the peri-urban zone is the only type of landscape directly connected to urban zones, we included data from peri-urban points in the model as the most plausible optimum state of bird richness that some of the urban environments could reach if they have similar features to the peri urban areas, i.e. larger vegetated area. The regression model that best fits the data (AIC = 4.3) is a linear or one order polynomial regression, which indicates a significant positive influence of vegetation cover area on the richness of bird species (Fig. 6). About 39.2% of data on bird richness in urban and peri-urban landscape matrices can be explained by forested areas (R² = 0.392; F = 18.045; p < 0.001).

Figure 6.

Linear regression model fitted to data of bird richness and vegetation cover area, considering restored forests and similar vegetation types in n = 30 surveyed points located in urban (filled circles) and peri-urban (empty circles) zones. This model can be represented by the equation 0.223x+0.542 (log-transformed data).

Through the ecosystem services (ES) assessment, it was possible to identify clear patterns of intensity spatially distributed according to land cover and land use, mainly derived from a gradient from natural to urbanised areas. While air pollutant removal, cooling effect and bird richness were strongly associated with vegetation cover (well-structured vegetation, for example, forests), the nature-based recreation opportunities, which also depend on human infrastructure, presented the highest values in urbanised areas, indicating a trade-off amongst the availability of these ES.

The assessment data map (Fig. 3) favours the comprehension about the differences amongst points and areas in terms of the distribution of bird richness and the intensity of ES provision in a land-use gradient composed by urban, peri-urban, farmland, transition and natural landscape matrices. Each point encompasses the current land use at a medium scale (500 m buffer), that was used to calculate ES values. Hence, it can be asserted that the sense of well-being an individual may experience varies across the river basin due to the distinct combinations of ES provided in different locations, thereby promoting diverse benefits (Rendón et al. 2019).

Looking at the current state of the ecosystem services, for the cooling effect, it is possible to see a positive scenario considering that over half of the sampled areas have a high level of cooling, as shown in Fig. 3. On the other hand, there was a clear drop in this effect in the urban area, the place with the highest demand for this service by citizens, since it has a higher human population density of the municipality (IBGE 2010, Costa et al. 2013) and its physical characteristics, including low vegetation cover, enhance the possibility of heat island phenomena (Sandifer et al. 2015). There are some vegetation patches in assessed urban areas, but only small fragments or green areas with low canopy cover (e.g. anthropogenic fields, parks, gardens).

The analysis of the PM₁₀ removal throughout the Monjolinho Basin indicated better removal where there is more vegetation cover. The air pollutant removal has a direct relationship with the leaf area of plants, which depends on vegetation type and its conservation state (Manes et al. 2016, Fusaro et al. 2017). These results pointed to a bad current situation, in which the places where the greatest amounts of pollutants are produced are precisely those where there is the greatest demand and the lowest availability of the pollutant removal services promoted by vegetation.

The current state of the nature-based recreation model indicated more opportunities in urban areas, despite these areas presenting the lowest levels of cooling effect and air pollutant removal and comparatively low bird richness. This scenario highlights important issues with the green infrastructure in the urban areas. Considering that the locales intended for recreation are the places where there is a concentration of people using them, the low levels of regulation ES is an indicator of the need to improve green infrastructure (i.e. tree cover) in those areas in order to increase and restore ES levels (Elmqvist et al. 2015, Fusaro et al. 2017, Cortinovis et al. 2018). The assessment of nature-based recreation opportunities shows better opportunities in the most urbanised areas, followed by the peri-urban zone, with a lack of opportunities in areas beyond city limits towards rural and natural areas (Fig. 3).

Even with the landscape types influencing the bird richness and ecosystem services levels, it is possible to notice in the maps one region on the north, where we found medium and high levels of bird richness and the other ES. That area can be considered a biodiversity and ES hotspot in the studied watershed. It is an important area with a large continuous vegetation fragment, belonging to the Federal University of São Carlos and to a neighbouring Zoopark (Parque Ecológico de São Carlos). The natural fragment has some trails that are used by people for exercise and by University students to do their research. The Zoopark attracts local people and visitors of other municipalities for recreation and education about wildlife. At the same time, the presence of conserved vegetation not only regulates the climate and reduces the air pollutants at the local level, but also can play a role as an important peri-urban barrier for the air pollutants (Fusaro et al. 2017). That vegetation also harbours at least two water springs of the Monjolinho River, the main river of the Basin, which is one of the main water sources for local agriculture and human consumption. It is a good example of a multifunctional landscape offering many ecosystem services for people's demands (Zasada 2011, Santos-Martín et al. 2019).

The current state of bird richness and ecosystem services shown in Fig. 3 indicated an apparent synergy amongst the bird richness, cooling effect and air pollutant removal, but there is a trade-off with nature-based recreation opportunities. These trends were confirmed by the correlation statistics (Fig. 4). This synergy between cooling effect and air pollutant removal was expected due the relationship of both models with canopy cover and the consequent synergy with bird richness once this vegetation contributed with more habitats for birds (MacGregor-Fors and Schondube 2011, Fusaro et al. 2017, Zardo et al. 2017). However, the trade-offs between nature-based recreation and other ES were unexpected, once the recreation opportunities considered here are according to the nature-based spectrum, so a better green infrastructure (i.e. the presence of vegetation, tree and canopy cover, environmentally-friendly playground) was expected in these places destined for recreation (Wen et al. 2018). This trade-off scenario could be reversible with nature improvement in the urban areas and/or recreation improvement in natural areas, which could be reached by planning a multifunctional landscape supported by Nature-based solutions (Cortinovis et al. 2018, Escobedo et al. 2019). Those solutions must consider offering more ecosystem services promoted by vegetation according to the population demand and create or enhance bird habitats throughout the Monjolinho Basin.

The variations of ES indicators throughout the Basin were driven by the landscape matrix, i.e. the predominant land use inside the assessed areas. The mean values of ES in urbanised areas significantly differed from all other landscape matrices, indicating that anthropogenic environments are the most distinct and can cause intense impact on services offered by nature.

The influence of the urban-rural-nature gradient on bird richness and ecosystem services was clearly shown by this assessment (Figs 3, 4). It is possible to notice two main patterns. The first was that the gradient flows in a positive direction from landscapes more impacted to less impacted or from urban (worst scenario) to nature (better scenarios with more suitable habitats) for bird species richness and the availability of regulation services (cooling effect and PM₁₀ removal). Additionally, in the second, the gradient drives the results in the opposite direction, being positive from less impacted to more impacted landscapes or from (worse) nature to urban (better) for recreation services, even though considering the nature-based recreation opportunities.

The lower of bird richness in urban areas may further increase the negative picture for ES in these areas, as birds also have their ecological roles that can be converted to ES and greater bird species richness levels may provide a greater richness of ecological functions and a consequent supply of ES (Şekercioḡlu 2006, Lessi et al. 2020). Furthermore, many Neotropical birds in the assessed community are seed dispersers, i.e. they contribute to forest maintenance, as well as the enhancement of current vegetation and the recovery of degraded areas (Wijedasa et al. 2020), which could lead to an enhancement of other ES, including air pollutant removal and cooling effect.

The different responses to the gradient and the different correlations show that the difference on the land use of point surroundings have a major impact on biodiversity and the supply of ES. Thus, the ES and bird richness assessed in this study follow the gradient of the landscape, decreasing from natural to urban (except recreation services), as found in other surveys (Haase et al. 2012, Hou et al. 2015), while population density and demand for ecosystem services increase in the same gradient. The trade-off relationship between recreational opportunities and climate-regulated ES in this type of gradient was also evidenced by Haase et al. (2012), as well as the trade-off relationship between recreation and biodiversity in different landscapes (Turkelboom et al. 2018). In this way, the loss of ES in the natural-rural-urban gradient and, at the same time, the increase in the demand for ES due to the increase in the population with the increase of urban landscape in this gradient are clear. The substitution of areas with vegetation for anthropogenic uses and the lack of policies that encourage the preservation or restoration of areas with native vegetation in urban areas, such as the surrounding streams, are identified as the main causes of this lack of ES supply. On the other hand, once the distribution of ES and biodiversity are diagnosed, we may assess how different areas could interconnect in a compensatory system for ES supply and biodiversity aiming to benefit areas with less ES offerings and more human demand.

It is known that urban areas do not support the same biodiversity as natural landscapes, but with the improvement of green areas and the creation of heterogeneous landscapes using native vegetation and wetlands, it is possible to increase the support capacity of the environment, increasing biodiversity rates and ES supply. On the other hand, the large areas of native vegetation close to urban zone (peri-urban and rural landscapes) can play a role as local hotspots for both biodiversity and ES, a situation that can be used for planning nature-based solutions, since those hotspots can serve as source areas to support improvements in the most degraded environments. Likewise, the cooling effect and air pollutant removal services within the urban zones are important for climate regulation and air quality in the landscape scale (Larondelle and Haase 2013, Fusaro et al. 2017, Cortinovis and Geneletti 2018).

The analyses of the landscape gradient helped to understand how to plan a compensation system, based on policies for the preservation of natural areas and the enrichment of urban areas using green infrastructure. The identification of synergies, trade-offs and the land-use types that contribute to the ES provision allows policy-makers to better understand the hidden consequences of preferring one ES to another (Haase et al. 2012). In addition, the results indicated a necessity to develop conservation planning (Chan et al. 2011), especially in areas subject to intensive human activities (Wu et al. 2013). Planning a multifunctional landscape can be the solution for implementing actions and creating infrastructures to contribute to ES supply and biodiversity conservation (Peña et al. 2018).

The results made it evident that simple public policies, such as tree planting as a nature-based solution to urban problems, can achieve good results in short and medium-term periods (Fig. 5 and Table 1). However, this conclusion was only possible after mapping the environment and proceeding with an assessment of current ES, which allows us to understand the interactions amongst landscape elements and ecosystem services.

The selected ES indicators for modelling the future scenarios, i.e. air pollutant removal and cooling effects, depend on green areas, mainly the vegetation type and structure, which are also precursors of bird assemblage structure and the services they can offer. In this context, the most valuable nature-based solution to enhance these three ES in areas with low values of indicators and to provide contact with nature during recreation activities inside the city, is the restoration of non-vegetated areas. If public sites without arboreal vegetation that are abandoned and many are extremely degraded, were used to create green areas with well-structured vegetation in urbanised landscapes (urban and peri-urban zones), our results showed that the three ES provided exclusively by nature could be significatively enhanced (Table 1) and consequently people's well-being (Sandifer et al. 2015, Rendón et al. 2019).

Nature-based recreation opportunities are not determined by forested areas as the other ES are. Instead, urban green infrastructures can be used to increase the opportunities, but only if there is planning focused on recreation, generally associated with infrastructure which allows such activity (Cortinovis et al. 2018). In this sense, it is not viable to try to estimate this ES indicator in scenarios with land restoration efforts. However, it is possible to state that the enlargement of vegetation cover area can serve for both, to promote a higher quality of recreation at sites already being used for this purpose (e.g. outdoor exercises, birdwatching) (Whelan et al. 2008, Rendón et al. 2019, Teixeira et al. 2019) or to create potential new areas for future use after a well-designed plan, focused on providing opportunities for the population (Sandifer et al. 2015, Cortinovis et al. 2018, Frantzeskaki 2019).

There are already several techniques for planning, decision-making and project execution, i.e. nature-based solutions, to improve the urban and agricultural scenarios of ecosystem services (Kabisch et al. 2016, Burkhard and Maes 2017, Zulian et al. 2017) and bring to these areas a greater supply of ES and a better quality of life for people as well. Nature-based recreations can even be one of the planning guidelines for bringing more biodiversity to the cities and achieving higher ES rates in those centres that are the places with the highest ES demand because of the higher concentration of people (Cortinovis and Geneletti 2018).

The synergy found between cooling effect and air pollutant removal services (Fig. 5) allows an evident improvement of these services in the future scenarios evaluated. Despite this, we see that the improvement is not achieved in the same proportion for all services (Fig. 5) mainly for recreation opportunities that depend on other actions and gain in biodiversity that depends on several factors (Belaire et al. 2015, Elmqvist et al. 2015, Frantzeskaki 2019).

The synergy present between ecosystem services allows a win-win relationship in a future scenario on well-planned landscapes, as discussed above. Additionally, it is important to consider a scenario based on native vegetation with high diversity to support biodiversity (Freitas et al. 2020, Berthon et al. 2021). Contrarily, it could cause a loss of synergy and create a trade-off relationship between services and biodiversity (Haase et al. 2012), as occurs in monoculture sites (Barros 2017, Berthon et al. 2021).

The assessment, based on micro-basins, may facilitate the decision-makers to access results of each section of the main basin, to rank the priorities to better address policies at small scales (due to a common lack of investment) and to meet local demands in each area. A focused design which considers the results of an ES assessment can play an important role as a facilitator instrument for natural resources management (Campanelli 2012, Tundisi and Scheuenstuhl 2012). As an example, this study shows Tijuco Preto and Gregório micro-basins with low and medium levels of ES as shown in the maps in Fig. 2 and indicate a lack of vegetation in the middle of their urbanised area (Table 1 and Suppl. material 1) even after the improvement in the future scenarios. Therefore, the assessment indicates that these regions should be prioritised in planning of public urbanisation policies. We found that only one of the possible tasks that ES mapping and assessment can perform to support urban planning and how we can activate that using a nature-based solution. From that, we can perform a range of elements, including the preferences, perspectives and demands from stakeholders and then apply other possibilities and multi-criteria analysis for even more accurate planning (Adem Esmail and Geneletti 2018, Geneletti et al. 2020).

In the end, some limitations on the ES assessment and analysis must be acknowledged. Part of the methodology was adapted to deal with the study area characteristics to run all models. There is a lack of studies on the same topic that could enhance and/or support the model development. The adaptation and mapping of land use, soil and canopy cover classification and the lack of an official map in high resolution could generate inaccurate data and assessments. Additionally, the mapping of the recreation opportunities was made in the laboratory by OpenStreetMap (OMS) infrastructure and confirmed in loco, but it would require the active involvement of local stakeholders able to score the elements, based on local knowledge and experience.