Ecosystem functioning in urban grasslands: The role of biodiversity, plant invasions and urbanization

Summary

Urbanization is driving the transformation of natural and rural ecosystems worldwide by affecting both the abiotic environment and the biota. This raises the question whether urban ecosystems are able to provide services in a comparable way to their non-urban counterparts. In urban grasslands, the effects of urbanization-driven ecological novelty and the role of plant diversity in modulating ecosystem functioning have received little attention. In this study, we assessed the influence of biodiversity, abiotic and biotic novelty on ecosystem functioning based on in situ measurements in non-manipulated grasslands along an urbanization gradient in Berlin (Germany). We focused on plant aboveground biomass (AGB), intrinsic water-use efficiency (iWUE) and ¹⁵N enrichment factor (Δδ¹⁵N) as proxies for biomass production, water and N cycling, respectively, within grassland communities, and tested how they change with plant biogeographic status (native vs alien), functional group and species identity. Approximately one third of the forb species were alien to Berlin and they were responsible for 13.1% of community AGB. Community AGB was positively correlated with plant-species richness. In contrast, iWUE and Δδ¹⁵N were mostly determined by light availability (depicted by sky view factor) and urban parameters like the percentage of impervious surface or human population density. We found that abiotic novelty potentially favors aliens in Berlin, mainly by enhancing their dispersal and fitness under drought. Mainly urban parameters indicating abiotic novelty were significantly correlated to both alien and native Δδ¹⁵N, but to AGB and iWUE of alien plants only, pointing to a stronger impact of abiotic novelty on N cycling compared to C and water cycling. At the species level, sky view factor appeared to be the prevailing driver of photosynthetic performance and resource-use efficiency. Although we identified a significant impact of abiotic novelty on AGB, iWUE and Δδ¹⁵N at different levels, the relationship between species richness and community AGB found in the urban grasslands studied in Berlin was comparable to that described in non-urban experimental grasslands in Europe. Hence, our results indicate that conserving and enhancing biodiversity in urban ecosystems is essential to preserve ecosystem services related to AGB production. For ensuring the provision of ecosystem services associated to water and N use, however, changes in urban abiotic parameters seem necessary.

Keywords:

Species diversity – Biodiversity – Urban ecosystems – Urban ecology – Grasslands – Invasive species – Urban environments – Ecosystem functioning

Introduction

Increasing human population and urban sprawl are driving the fast transformation of natural and rural environments into urban ecosystems worldwide [1]. Over the past six decades, 25% of the global population has moved from rural into urban settlements and this trend is expected to continue [2]. Urbanization encompasses significant changes in all components of urban ecosystems and their functioning [3], such as the so-called urban heat island effect [4], higher concentrations of heavy metals in the soil and soil compaction [5–7] increased atmospheric nitrogen deposition [8], habitat fragmentation and isolation [9] as well as the development of transportation networks that favor the dispersal of alien plants [10, 11]. Urbanization is also a powerful driver of change to biodiversity patterns, associated with filtering species according to their pre-adaptation to urban environments [12–14]. While alien plant species are frequently found in different urban settings, particularly in novel urban ecosystems [15], cities can also harbor a considerable richness of native plant species [12, 16]. However, urbanization has been also linked to a loss of biodiversity, ecological homogenization and changes in community composition, including a higher proportion of alien species pre-adapted to novel habitats [17–20] and species compositions that cannot be found in near-natural environments.

Consequently, species in urbanized settings face situations to which they have not been exposed during their evolutionary history. Urbanization leads to the transformation of (near-)natural into novel urban ecosystems that prevail in many cities and can be characterized by both novel habitat configurations and novel species assemblages without analogues to natural landscapes [18]. These conditions thus can be characterized as ‘ecologically novel’ [21], and ‘ecological novelty’ in this sense can be seen as an overarching conceptual framework to address the complex eco-evolutionary consequences of anthropogenically induced changes for organisms and ecosystems. In the following, we use the term ‘abiotic novelty’ to characterize environmental conditions that differ from near-natural grassland conditions typical for a particular area. ‘Biotic novelty’ refers to the presence of species, and thus interaction partners, that have not been present in the area in historic times (i.e. alien species).

Grasslands in rural settings provide a wide variety of relevant ecosystem services such as microclimate regulation, water flow regulation and nutrient cycling [22]. Grassland vegetation supports the sequestration of carbon dioxide and soil formation through primary production and modulates hydrological flows by fluctuations in its water use. Moreover, through the acquisition and storage of nitrogen (N), grassland vegetation plays a major role in terrestrial nitrogen cycling. Worldwide, grasslands make up a considerable area in cities. In Berlin, for instance, they cover as much as 5% of the city area [23]. Urban grasslands encompass meadows and lawns in domestic gardens, parks, vacant land, remnants of rural landscapes, and areas along transportation corridors [23–28]. The role of urban grasslands for biodiversity conservation has gained increasing attention in the last decade [23, 24, 29, 30]. Yet in urban grasslands, the role of plant diversity in modulating ecosystem functioning, including proxies such as aboveground biomass (AGB) production or resource use has received little attention [24, 31].

A positive relationship between biodiversity and ecosystem functioning (often exemplified by AGB and less frequently by plant water or N use) has been widely described in non-urban experimental grasslands resembling natural or agricultural grassland communities [22, 32–39]. This relationship has been also found in experimental urban grassland assemblages resembling North America turfgrass communities [40] and experimental grassland assemblages resembling the soil and vegetation on abandoned urban sites in Berlin [41]. However, it remains unknown whether the positive link between biodiversity and ecosystem functioning prevails in urban grasslands where biodiversity and other environmental factors are not manipulated (i.e. non-experimental grasslands)[42]. Ecosystem functioning in urban grasslands is critically understudied, with only a few studies focusing on urban lawns as a subgroup of urban grasslands (e.g. [20, 40, 43, 44]). Furthermore, studies comparing ecosystem functions of urban vs non-urban grasslands are rare (but see [43, 44]). One of these few studies has shown that in the Front Range of Colorado (USA), lawns were more productive than native grasslands, likely as a result of increased irrigation or fertilization associated to management [44].

Along with biodiversity, both biotic novelty (the presence of alien species) and abiotic novelty (altered environmental conditions) can be expected to ultimately affect ecosystem services by modifying ecosystem functioning. For example, strongly altered urban soils may be less productive than rural soils; and a high abundance of alien plants may modulate biomass production or nutrient cycling–with negative effects on ecosystem services [45]. Plant invasions have been often associated with negative effects on native plant and animal species, ecological functions and ecosystem services [45–47]. For instance, a global meta-analysis of field studies assessing the impact of alien plants revealed that overall species diversity declined but primary production increased, probably due to a sampling effect [47]. Yet, invasion effects on ecological functions are clearly understudied in urban environments (e.g., [31, 46]). Given the high abundance of alien plants in urban environments, and particularly in novel urban ecosystems [15], enhanced insights into biodiversity-ecosystem functioning relationships in cities are important for urban planning and biodiversity conservation, and ultimately for human well-being.

Nonetheless, the available information on the role of ecological novelty in influencing ecosystem functioning proxies such as AGB production, intrinsic water-use efficiency (iWUE) and N cycling in urban grasslands is limited. Furthermore, the relative role of biotic versus abiotic factors, and their significance in comparison to biodiversity as a driver of ecosystem functioning is unclear. Evidence on the key drivers of various ecosystem properties compiled for a wide range of aquatic and terrestrial ecosystems is diverse ([36] and references therein): biotic changes were generally found to have an equally important effect on ecosystem properties than abiotic changes [48]. For instance, in boreal forests, changes in microclimate (i.e. decreased heat flux into the soil) induced by moss cover may modulate permafrost stability in a comparable way to environmental change. Conversely, abiotic conditions may sometimes offset the effects of species richness on ecosystem processes, as shown for instance by a theoretical modelling study revealing that abiotic conditions mask the effect of plant species richness on plant productivity [49]. Such ambiguity may be due to different mechanisms acting upon different groups of plant species, such as functional groups or native versus alien species, or different mechanisms acting during different periods of the growing season.

Despite the fact that the performance of individual species, or functional groups of species, might ultimately determine ecosystem level functioning (e.g. [50, 51]), current knowledge on species-specific responses to biodiversity and ecological novelty [52] and on seasonal variations in such response is scarce, especially in the context of urban grassland research. In this regard, the combined estimation of multiple photosynthetic parameters as well as leaf stable isotope composition can provide insights into the physiological mechanisms controlling plant individual´s performance and AGB production, especially when measured over the course of the growing season [53]. Chlorophyll fluorescence, for instance, reflects the plant ability to put up with abiotic stress factors and indicates the damage level of the photosynthetic apparatus [54], while leaf C isotope composition reflects the average intrinsic water use efficiency throughout the time period during which the organic matter assessed was formed [55].

Based on the positive relationship between diversity and ecosystem functioning found in numerous experimental grasslands and semi-natural grasslands, namely between diversity and AGB production [22, 33, 37, 40], diversity and water use [50, 56] as well as diversity and N cycling [57], we hypothesize that species richness is the main controlling factor for plant community AGB, iWUE and N use in urban grasslands, overriding both abiotic and biotic novelty effects. In this study we aim to improve the understanding of the relative importance of biodiversity, abiotic and biotic novelty for the functioning of urban grasslands. To our knowledge, this is the first study to do so based on in situ measurements in non-manipulated urban grasslands. Indeed, previous results based on experimental urban grasslands pointed to the importance of such field assessments [40]. In particular, we aim to:

• identify the main factors (i.e. parameters related to biodiversity, abiotic and biotic novelty) driving the variability in AGB, iWUE and N cycling of whole plant communities, plant functional groups and plants with different biogeographic status (native vs alien) across 20 urban grasslands located in Berlin;

• assess the biodiversity-ecosystem functioning relationship in the studied urban grasslands and compare it with that found in natural and agricultural grasslands;

• test if parameters related to AGB production, iWUE and N use at the species level significantly differ between the mid- and peak growing season and assess the relationship between those parameters and biodiversity, abiotic and biotic novelty in both time periods.

Material and methods

Study area and study system

The study was carried out in Berlin, Germany. Berlin has an area of 891.1 km² and a population of ≈ 3.6 mio. inhabitants [58]. Mean annual precipitation amounted to 576 mm and mean annual air temperature was 9.9°C for the period from 1981 to 2010 at the central city location of Tempelhof [59]. Our study system is dry grasslands according to the classification of the “Catalogue of Natural Habitats and Species of Appendices I and II of the Habitats Directive in Brandenburg” [60]. Dry grasslands are a vegetation type that spans a range of near-natural to strongly human-influenced sites throughout the city. For this reason, urban dry grasslands had been selected as a model ecosystem within the CityScapeLabs, an experimental platform established by the Berlin-Brandenburg Institute of Advanced Biodiversity Research (BBIB) with a network of permanent plots for the investigation of biodiversity and ecosystem functioning in urban environments [61]. For this study, we selected a subset of 20 dry grassland plots of 16 m² that were relatively evenly distributed across the city (Fig 1) and whose surroundings were subject to different levels of urbanization, indicated, e.g. by human population density or percentage of impervious (sealed) surface. These plots comprised several land use types (e.g. parks, cemeteries, forest clearings), but were selected to minimize human management inputs such as fertilization, irrigation and mowing. The parameters related to biodiversity and abiotic novelty at the study sites (with the exception of air temperature and air relative humidity) were provided by the CityScapeLabs, and the degree of biotic novelty of the study sites was calculated based on the CityScapeLabs biodiversity data (see Table 1 for more information).

Biodiversity

Vascular plant diversity was characterized both by taxonomic and functional diversity, estimated as species richness and Rao’s quadratic entropy (Rao´s Q, [62]), respectively. To assess species richness, vegetation surveys were carried out in a 4 x 4 m plot delimited within each of the 20 grasslands between April 18^th and May 19^th 2017. Based on expert knowledge and region-specific literature [63–65], plant species were classified according to their biogeographic status into native or alien. Species introduced after 1492 (neophytes) were considered ‘alien’, while species introduced before 1492 (archeophytes) and native species were merged into the category ‘native’. Thereafter, the proportion of alien species at each grassland plot was estimated. Rao's quadratic entropy was calculated using Gower distances between species pairs based on twelve plant functional traits: plant height, specific leaf area, life form, flower color, flower class, clonal growth organs, length of dispersal unit, seed mass, leaf area, leaf nitrogen content, nitrogen fixation and mycorrhizal infection [66]. Trait data were extracted from the TRY database [67] and the BiolFlor database [68]. Additionally, the percentage of plot surface covered by moss (moss cover) and litter (litter cover) was visually estimated.

Abiotic novelty

A total number of 22 environmental variables related to abiotic novelty were assessed in the 20 selected dry grasslands plots within the CityScapeLabs (Table 1). Abiotic novelty-related parameters included urban parameters (i.e. human population density, impervious surface, floor area ratio, etc.), climate and soil parameters assumed to be linked to urbanization (i.e. heavy metals in the soil, long-term average air temperature, etc.) as well as habitat continuity (age) and size of dry grassland biotope patch.

The age of the grasslands was assessed by digitizing and georeferencing historical land use maps (i.e. Preußische Uraufnahme (1831–71) and Preußische Neuaufnahme (1927–40)) and intersecting them with the current biotope mapping of Berlin and Brandenburg. Data on size of dry grassland biotope patch, share of grassland, impervious surface, FAR, PopD, RoadD, RailwD, RDist, RailwDist, long-term air T and UrbClim (see Table 1 for parameter description) were extracted from the Geoportal Berlin of the Senate Department for Urban Development and Housing [69]. A mean value for buffer areas of multiple sizes around the biotope patch was calculated for share of grassland, impervious surface, FAR, PopD, RoadD and RailwD using QGIS 2.18.0 [70]. RDist and RailwDist were also calculated using QGIS 2.18.0 as the shortest distance from the grassland plot midpoint to the nearest road or railway, respectively.

To analyze soil parameters (Table 1), 15 subsamples of 30 cm depth were taken with a 1.4 cm stainless steel soil corer in the aforementioned delimited 4 x 4 m plot within each grassland in June 2017, air-dried, sieved (2 mm) and then mixed and homogenized into a single sample. Total soil N (TN) was analyzed by adsorption chromatography [71], whereas soil Cu, Zn, Cd, Pb and Ni were analyzed by inductively coupled plasma optical emission spectrometry (iCAP 6000 Spectrometer, Thermo Fisher Scientific, Dreieich, Germany) following soil extraction using ammonium nitrate solution [72]. Gravimetric soil water content (W_C) was determined as % of the soil dry weight by oven drying the samples at 105°C until constant weight.

At each grassland plot, air temperature (AirT) and relative humidity (AirRH) were recorded every 30 minutes at a 20 cm height throughout the growing season (between April 21^th and 26^th November 26^th, 2017) with a data logger (OMEGA, OM-EL-USB-2). Subsequently, the mean AirT and AirRH for the whole growing season was calculated. Due to logger misappropriation, it was not possible to record AirT and AirRH either partly or entirely in 6 plots, where 50–100% of the data were missing. Therefore, we collected missing data from meteorological stations located within 0.36 to 5 km to the plots [73]. At the plots where only part of the data was missing, we confirmed that AirT and AirRH from Weather Underground differed less than 5% and 12%, respectively, from the data recorded with microclimate logger data. Sky view factor (SVF) is a parameter related to the urban heat island phenomenon, with higher values indicating lower long-wave radiation emissions of built surfaces to the sky during the night [74]. SVF was estimated as the share of open sky based on the analyses of images taken using a Canon EOS 700D camera coupled to a circular fisheye (4.5 mm F2.8 EX DC HSM, Sigma) at ≈ 1.5 m height on the plots. Namely, SVF was calculated using the subversion SOLWEIG1D of the SOLWEIG model (version 2015a) [75, 76].

Biotic novelty

The degree of biotic novelty of individual plots was estimated by a biotic novelty index (BNI). The BNI is a compound measure based on a modified version of Rao´s quadratic entropy and designed to quantify the functional ecological novelty of communities [66]. The index captures the functional diversity contributed by novel species recently arrived in the community, weighted by their relative abundance. It is based on two components: the pairwise functional distance between species, see [62, 77]) and a temporal coexistence component that weighs the functional differences between pairs of species based on how long both species have been present in the region. For example, if a given species pair consists of one native and one recently arrived alien species, the trait distance between both will receive a higher weight than the distance between a pair consisting of one native and one earlier arrived alien species. This idea is based on the finding that alien species will gradually become familiar with their interaction partner(s) over time [78–81] which may lead to a decrease of novelty in the community. The temporal coexistence component of the index was calculated from species’ residence times in the Berlin area. For the calculation of the functional diversity component after Rao, we used the same method and traits as described above. More information can be found in the S1 Appendix.

Plant community aboveground biomass and soil collection for isotope analyses

Total standing aboveground biomass (AGB) of the plant community was sampled between July 24^th and August 9^th 2017 in the 20 selected dry grasslands in Berlin. As grasslands are regularly managed by the local authorities, protective bands were placed around the 4 x 4 m plot in early spring 2017 to prevent mowing prior to sampling. Since the plots had been last mowed in autumn 2016, the AGB collected in late summer 2017 biomass integrated the biomass production throughout the 2017 growing season, including both alive and dead plant parts. Replicate samples (n = 3) were collected by clipping the vegetation to the ground in 20 x 50 cm quadrats randomly selected within each plot [82]. After identification to the genus or species level, the collected plants were classified into three functional groups: legumes (N-fixing dicots), forbs (non-N-fixing dicots) and graminoids (grasses, sedges and rushes). Next, they were separated into biogeographic status classes (native vs alien). Samples were transported to the lab within 6 hours, oven dried at 70°C to constant weight and weighed thereafter (Kern EW 620, Kern & Sohn GmbH, Balingen, Germany). For each plot, the AGB of different plant groups (i.e. graminoids (AGB_G), forbs (AGB_F) and legumes (AGB_L), natives (AGB_N) aliens (AGB_N), native graminoids (AGB_NG), alien graminoids (AGB_AG), native forbs (AGB_NF), alien forbs (AGB_AF), native legumes (AGB_NL), alien legumes (AGB_AL)) as well as the whole community AGB (AGB_C) was calculated as the average of the three replicates within the plot. Following aboveground biomass collection, a soil core of 30 cm depth was extracted with a 5 cm diameter soil core sampler (Wurzelbohrer V2A, Umwelt-Geräte-Technik GmbH, Müncheberg, Germany) at every quadrat replicate (n = 3) within each grassland plot. The upper 10 cm of the soil cores were homogenized, sieved (2 mm mesh size) and oven-dried at 105°C prior to isotope analyses.

Species level gas-exchange and chlorophyll-a fluorescence measurements

In order to assess the effects of biodiversity, abiotic novelty and biotic novelty on AGB production at the species level, we selected two plant species that were relatively common in the study area and belonged to the two dominant functional groups in terms of biomass: the grass Calamagrostis epigejos (L.) Roth and the forb Plantago lanceolata L.. Ecophysiological parameters were measured during two time periods: in spring (mid-growing season), between May 8^th and May 19^th 2017, and in summer (peak of the growing season), between July 27^th and August 9^th 2017, simultaneously to the collection of AGB and soil. Whereas in spring we selected individuals that were exclusively present in the 4 x 4 m plot delimited within each grassland, in summer we selected individuals that were also outside the plot in order to obtain a minimum number of replicates. These individuals were located within a distance of 10 m from the plot, in comparable vegetation and soil conditions and thus we do not expect them to be experiencing different conditions from the ones inside the plot. Therefore, the number of plots in which the species were present varied between months, with C. epigejos being present in 13 and 17, respectively, out of the 20 selected grassland plots in spring and summer, respectively, and P. lanceolata occurring in 9 and 12 plots, in spring and summer, respectively. In total, C. epigejos and P. lanceolta co-occurred in 40–55% of the study plots, which covered a comparable gradient of plant species richness, abiotic, and biotic novelty (S1 Table).

The ecophysiological parameters were measured with a portable gas exchange fluorescence system (Walz GFS-3000, Heinz Walz GmbH, Effeltrich, Germany) equipped with a clamp-on leaf chamber of 8 cm². Gas exchange measurements were carried out mostly on clear sunny days between 09:00 and 16:30 on unshaded, fully expanded mature and healthy leaves from three individuals per species in each plot/grassland. The leaf chamber conditions were set constant to 20°C leaf temperature, 60% relative humidity, 400 ppm CO₂ concentration and saturating photosynthetic photon flux density (PPFD_sat, 1500 μmol m^-2 s^-1). Gas exchange parameters were logged after the leaves reached steady-state conditions at saturating PPFD and comprised photosynthetic rate (A, μmol CO₂ m^-2 s^-1), transpiration rate (E, mmol H₂O m^-2 s^-1), stomatal conductance (g_s, mol air m^-2 s^-1) and ratio of the intercellular to ambient CO₂ concentration (c_i/c_a, μmol CO₂ μmol^-1 CO₂). Instantaneous water-use efficiency (instant-WUE, μmol CO₂ μmol^-1 H₂O) was calculated as A/E. Chlorophyll-a fluorescence parameters were measured with a fluorescence Module (LED-Array/PAM-Fluorometer 3055-FL, Heinz Walz GmbH, Effeltrich, Germany) integrated in the leaf chamber. These parameters included steady-state fluorescence of the light-adapted leave (F), steady-state fluorescence of the dark-adapted leave (F₀), fluorescence of the dark-adapted leaf during a saturating light pulse (F_m), fluorescence of the illuminated leaf when a saturating light pulse is superimposed on the prevailing environmental light levels (F_m´, [83]). The maximum potential quantum yield of photosystem II was calculated as F_v/F_m, where F_v = F_m—F_0. The effective quantum yield of electron transport through photosystem II was calculated as ΔF/F_m′ = (F_m´ - F)/F_m´ [84]. Photosynthetic light-response curves of ΔF/F_m′ and apparent electron transport rate (ETR) were obtained by increasing PPFD intensity (i.e. 0, 100, 300, 600, 900, 1200, 1500 and 2000 μmol m^-2 s^-1) after 30 to 120 s (i.e. instant light curves, [85]). ETR was estimated as:

Carbon and nitrogen isotope composition

The oven dried leaves and AGB samples were finely ground with a Vibrator Disc Mill RS 200 and a Mixer Mill MM 200 (Retsch GmbH Haan, Germany). About 1.5 mg of the powdered material from each sample was placed into a tin (Sn) capsule (IVA Analysentechnik, Meerbusch, Germany) for analysis of carbon isotope composition (δ¹³C, ‰) nitrogen isotope composition (δ¹⁵N, ‰), total nitrogen content (N%, %) and total carbon content (C%, %). The analyses were performed at the ZALF Stable Isotope Facility, Müncheberg, Germany. C%, N%, δ¹³C and δ¹⁵N were determined using a Flash 2000 HT Elemental Analyzer (Thermo Fisher Scientific, Bremen, Germany) coupled to a Delta V isotope ratio mass spectrometer (IRMS) via a ConFlo IV interface (both Thermo Fisher Scientific, Bremen, Germany). In the case of the soil samples extracted following AGB collection, 20 mg were placed into a tin (Sn) capsule (Säntis Analytical AG, Teufel, Switzerland). The analyses were performed at the WSL Stable Isotope Research Centre (SIRC), Birmensdorf, Switzerland. C%, N%, δ¹³C and δ¹⁵N were determined using an EA1110 Elemental Analyzer (CE Instruments, Milan, Italy) coupled to a Delta Plus XL ratio mass spectrometer via a ConFlo II interface (both Thermo Fisher Scientific, Bremen, Germany).

Reported isotope ratios were calculated as:

The average δ¹³C and δ¹⁵N of each plant group in a given plot (δ¹³C_FGi and δ¹⁵N_FGi by functional group, biogeographic status or by their combination) was calculated as the C-weighted and N-weighted average of every replicate in the plot, respectively. In order to estimate the average δ¹³C and δ¹⁵N at the community level (δ¹³C_C and δ¹⁵N_C) we applied a stable isotope mixing model were δ¹³C and δ¹⁵N of each plant group was scaled by the group specific contribution to the total C and N pool. We thus assumed that the contribution of each plant group to δ¹³C_C and δ¹⁵N_C was proportional to the relative contribution of that particular group to the total community C pool and N pool, respectively, as described by [87]. That is, average δ¹³C_C and δ¹⁵N_C were calculated as the C-weighted and N-weighted average, respectively, of all the plant groups present in a plot:

As a next step, carbon isotope discrimination, which is the carbon isotopic ratio in plant tissue relative to that of the atmosphere, was calculated at the plant group, community and species level as:

To correct δ¹⁵N values for site-specific differences in background bulk soil δ¹⁵N, we estimated the ¹⁵N enrichment factor of biomass compared to soil background values [94] at the community (Δδ¹⁵N_C), plant group (Δδ¹⁵N_FGi) and species level (Δδ¹⁵N_Spi) as:

Data analysis

The relationship between species richness and AGB, iWUE or Δδ¹⁵N was assessed by regression analysis. To test if the relationship between biodiversity and AGB in Berlin grasslands was comparable to that found in multiple experimental grasslands across Europe, we extracted numerical data from plots in Fig 3 in [36] using Web Plot Digitizer v4.1 (https://apps.automeris.io/wpd/). To ensure that differences in the biodiversity-AGB relationship were not stemming from environmental differences among grassland locations, the extracted data were min-max normalized. The slopes from the corresponding biodiversity-AGB significant regression or ANOVA models were then synthesized.

Random forest (RF) analyses was used to identify the main factors driving the variability in AGB, iWUE and Δδ¹⁵N across the studied urban grasslands in Berlin at the plant group (i.e. concerning functional group or biogeographic status) and community level and to thereby assess the relative contribution of biodiversity, biotic and abiotic novelty in controlling urban grassland functioning. These factors included functional and taxonomic diversity, the BNI, environmental variables related to abiotic novelty but also moss and litter cover in the grasslands. We used a conditional inference random forest algorithm [96, 97] to deal with missing values, nonlinear associations among variables, and a number of predictors larger than the number of samples. RF was applied for each response variable to estimate the relative importance of the predictors, measured as the contribution to model fitting accuracy of R² with variable selection [98]. RF can robustly estimate the relative importance among highly correlated predictors (e.g. land use cover proportion at different buffer radii; [99–101]. As small sample size can result in instable estimates [102], we took the following strategies. Firstly, we determined the number of tree hyper-parameter (ntree = 2000) after confirming its performance stability compared with lower numbers (10, 100, and 1000). Secondly, we estimated variable importance 100 times with changing random seed for randomization and then took an average. Finally, we estimated statistical significance for each predictor, i.e. the probability that the averaged importance score is obtained just by chance (999 permutation; alpha = 0.05) [102, 103]. Thereby, significant predictors were selected and their explanatory power (R²) was estimated for each response variable. For non-significant predictors, variable importance was set to 0 [103]. Modeled association patterns between significant predictors and each response variable were visualized using partial dependence plot (PDP) approach, which averages out the other predictors’ effects to confirm the direction of the effects but not the effect size [104]. All procedures were done in R v3.5.1 [105], with ‘party’ and ‘mlr’ packages [106, 107]. The R script including the other hyperparameter values is available at [108].

It was not possible to obtain a RF model for every plant group (an overview of the RF models developed is provided in S2 Table). This was partly due to the sparse presence of some functional groups at the plots (i.e. alien graminoids and alien legumes), but also due to a relatively limited sample size and the lack of explanatory power of the predictors; we emphasize, however, that the sample size (i.e. the number of sites) was relatively large for such a field study, where all plots had to be sampled and measured within a short time frame for comparability. Likewise, the available data was not sufficient to identify the main factors driving the variability in photosynthetic parameters, iWUE and Δδ¹⁵N at the species level by using RF model analysis. Therefore, we opted for analyzing the correlation (Spearman´s rho, r_s) between AGB, iWUE or Δδ¹⁵N at the species level and species richness, SVF or BNI as representative parameters for biodiversity, abiotic and biotic novelty, respectively. Mann-Whitney or Kruskal-Wallis with pairwise multiple comparison tests were used to analyze significant differences (p < 0.05) in iWUE and Δδ¹⁵N among plants with different biogeographic status (i.e. native vs alien) or plant from different functional groups. The Mann-Whitney test was also used to analyze seasonal differences in photosynthetic parameters, iWUE and Δδ¹⁵N at the species level. The aforementioned analyses were performed with IBM^® SPSS^® Statistics 22.0.0.0 Software.

Results

Variability in aboveground biomass, intrinsic water use efficiency and ¹⁵N enrichment factor

The 20 urban grasslands studied in Berlin covered a wide range of plant species richness and functional diversity as well as biotic and abiotic novelty (Table 1). Plant species richness varied more than threefold among grasslands and the BNI (biotic novelty index) showed a tenfold variation. Likewise, soil parameters, climate and urban parameters showed wide value ranges across sites. Total plant community aboveground biomass (AGB_C) varied threefold across the studied urban grasslands in late summer 2017 (Fig 2A). Important variations were also found between replicates within plots in AGB_C (CV = 5–66%, mean 30%). Overall, native species dominated AGB_C, contributing more than 60% in 18 out of 20 plots (Fig 2B). Graminoids and forbs made up for most of AGB_C (on average >90%) and their share was similar, whereas legumes contribution was relatively small (Fig 2B). All legume and graminoid species present in the sampled quadrats were native to Berlin with the exception of Medicago x varia T. Martyn and Eragrostis minor Host, which were found in one and two plots, respectively. Indeed, the average contribution of alien legumes and graminoids to AGB_C across the 20 studied urban grasslands was 1.4% and 0.9%, respectively. However, approximately one third of the forb species were alien to Berlin, and these were present in most plots, contributing on average 13.1% to AGB_C in the study sites. An overview of the species present in the sampled quadrats is provided in S3 Table. Both iWUE_C and Δδ¹⁵N_C varied markedly among urban grasslands, having an average value of 50.84 μmol mol^-1 (SE = 0.44) and -3.43 ‰ (SE = 0.37), respectively. Plants with different biogeographic status did not significantly differ in their iWUE or Δδ¹⁵N (Mann-Whitney, p > 0.05; Fig 3A and 3B). However, iWUE_G was significantly higher than iWUE_F and iWUE_L (Kruskall-Wallis, p < 0.05; Fig 3C) and Δδ¹⁵N_G was significantly lower than Δδ¹⁵N_F and Δδ¹⁵N_L (Kruskall-Wallis, p < 0.05; Fig 3D).

With regard to the seasonal differences at the species level, in C. epigejos A increased significantly from spring to summer, whereas instant-WUE decreased (Mann-Whitney, p > 0.05; Fig 4A and 4C). In P. lanceolata, both instant-WUE, F_v/F_m and ΔF/Fm′ significantly decreased from spring to summer (Mann-Whitney, p > 0.05; Fig 4C, 4E and 4F). Whereas iWUE_C.ep significantly increased from spring to summer, iWUE_P.la significantly decreased over the growing season but no seasonal differences were found for Δδ¹⁵N_C.ep or Δδ¹⁵N_P.la (Fig 4H and 4I). Intrinsic photosynthetic capacity, as indicated by the light-response curves and their associated cardinal points, did not show important seasonal variations in the case of C.epigejos but decreased from spring to summer in P.lanceolata (Fig 5).

Key predictors for the functioning of urban grasslands

In the following, we report the RF models´ results obtained for AGB, iWUE and Δδ¹⁵N. For each ecosystem functioning proxy, we first described the RF models´ results for the whole plant community, followed by those for plant biogeographic status classes, plant functional groups and finally for the combination of the latter two. As mentioned above, it was not possible to obtain a RF model for every plant group. Therefore, the RF model description is not consistent among ecosystem functioning proxies and plant groups and we decided to only include results of those plant groups that were more consistently modelled in Fig 6. The results for all the plant groups are provided in the S1 File, S2 File and S3 File.

Aboveground biomass

The explanatory power (R²) of the AGB_C, AGB_N and AGB_A random forest models were 0.87, 0.83 and 0.76, respectively, with a corresponding validation accuracy (R²) of 0.36, 0.41 and 0.17. Out of 27 explanatory variables, six were selected for the AGB_C RF model. Species richness was the most important predictor (R²contribution: 41%), followed by size of dry grassland biotope patch (16%), moss cover (12%), functional diversity (8%), SVF (7%) and litter cover (3%). While AGB_C was predicted to markedly increase above 25 plant species and a dry grassland biotope patch size of 3500 m², AGB_C diminished with increasing moss cover, being this trend especially strong at moss cover < 10% (Fig 6A; Figure A in S1 File). The five variables selected for the AGB_N RF model partially overlapped with those selected for the AGB_C and included in descending order of importance: size patch, litter cover, relative humidity, species richness and N soil content. As for the AGB_C RF model, the AGB_N RF model predicted a monotonic increase with species richness (8%), size of dry grassland biotope patch (40%) and litter cover (22%; Fig 6B; Figure B in S1 File). In the case of AGB_A, most predictors selected for the RF model were related to abiotic novelty and specifically to urban parameters. In general terms, AGB_A increased monotonically with both PopD and FAR in buffer areas of multiple sizes around the biotope patch (R²: 59% and 11%, respectively; Fig 6C; Figure C in S1 File). Functional diversity influenced AGB_A positively as well (Figure C in S1 File), but its explanatory power was low (6%).

The explanatory power (R²) of the AGB_G and AGB_F RF models was 0.57 and 0.72, with a corresponding validation accuracy of 0.17 and 0.30. AGB_G decreased with increasing moss cover, which was the only predictor selected for the AGB_G RF model (R² contribution: 57%; Fig 6D). Species richness was the main predictor for AGB_F (53%), followed by functional diversity (17%) and litter cover (3%). For the two strongest predictors, the RF model revealed a hump-shaped relationship (Fig 6E; Figure D in S1 File). AGB_NG, AGB_NF and AGB_AF revealed a similar pattern to that observed for the overall plant groups with different biogeographic status (Figures E-F in S1 File). AGB_NG showed a strong monotonic negative dependence on moss cover (61%; Figure E in S1 File). Whereas the size of dry grassland biotope patch and species richness were the main predictors for AGB_NF (55% and 10%, respectively), abiotic novelty-related variables explained the variance observed in AGB_AF (Figures F-G in S1 File). Specifically, AGB_AF increased monotonically with PopD in buffer areas of multiple sizes around the biotope patch (39%) but declined at low Pb (11%).

Intrinsic water use efficiency

The explanatory power (R²) of the iWUE_C, iWUE_N and iWUE_A RF models was 0.84, 0.85 and 0.88, respectively, with a corresponding validation accuracy of 0.50, 0.56 and 0.45. A total of 4 variables were selected for the iWUE_C RF model. SVF was the most important predictor (R² contribution: 62%), followed by Wc, share of grassland and Cu (8%, 7% and 6%, respectively). Community iWUE increased monotonically above a SVF of 0.8 and showed a rapid increase when the share of grassland was very low, reaching a maximum at ≈ 5% (Fig 6A; Figure A in S2 File). Likewise, iWUE_C varied markedly at the lower Cu and Wc ranges, but in this case its response was negative (Figure A in S2 File). The four predictors selected for the iWUE_N mostly overlapped with those selected for iWUE_C, with SVF being also the strongest predictor for iWUE_N (53%) and having a positive effect (Fig 6B; Figure B in S2 File). Additional parameters such as Wc, TN and Cu content accounted for 30% of R² and showed a negative e relation with iWUE_N (Figure B in S2 File). RailwD_5000 (47%) was the main predictor for iWUE_A, although SVF and AirT were also relevant for (26% and 15%, respectively). iWUE_A decreased monotonically with increasing RailwD_5000 but showed a positive response to SVF and AirT (Fig 6C; Figure C in S2 File). The explanatory power (R²) of the iWUE_G RF model was 0.83, with a validation accuracy of 0.49. iWUE_G increased with increasing SVF (49%), but decreased with Zn, W_C and Cu (14%, 13% and 8%, respectively; Fig 6; Figure D in S2 File). The iWUE_NG and iWUE_AF RF models revealed a similar pattern to that observed for the overall plant groups with different biogeographic status (Figures E-F in S2 File). iWUE_NG increased monotonically with SVF (39%) but decreased with W_C and N (14% and 10%, respectively. iWUE_AF decrease with increasing RailwD_5000 (53%) and showed significantly lower values in new plots (10%) (Figures E-F in S2 File).

¹⁵N enrichment factor

The explanatory power (R²) of the Δδ¹⁵N_C, Δδ¹⁵N_N and Δδ¹⁵N_A RF models was 0.71, 0.7 and 0.75, respectively, with a corresponding validation accuracy of 0.2, 0.22 and 0.3. The 8 predictors selected for the Δδ¹⁵N_C RF model in descending order of importance when summing up the explanatory power of the predictors for buffer areas of multiple sizes around the biotope patch were FAR (R² = 0.17), RailwD (14%), PopD (14%), impervious surface (12%) and share of grassland (7%), followed by UrbClim (3%), RoadD_100 (2%) and AirRH (2%) (Fig 6A; Figure A in S3 File). Note that when considering separately the explanatory power of the predictors for different buffer areas sizes around the biotope patch, impervious surface_100 was the main predictor selected by the model. The predictors selected for Δδ¹⁵N_N and Δδ¹⁵N_A mostly overlapped with those selected for Δδ¹⁵N_C, though their explanatory power differed (Figures A-C in S3 File). The main predictor for Δδ¹⁵N_N was RailwD in buffer areas of multiple sizes (34%), followed by impervious surface_100 (14%), share of grassland_5000 (9%), FAR in buffer areas of multiple sizes (8%), PopD_500 (3%) and AirRH (3%) (Fig 6B; Figure B in S3 File). However, PopD and FAR in buffer areas of multiple sizes (32% and 22%) were the strongest predictors for Δδ¹⁵N_A, followed by impervious surface, AirT_long term, UrbClim and share of grassland 500 (8%, 8%, 3% and 2%, Fig 6C; Figure C in S3 File)). All the aforementioned predictors showed a positive monotonic relationship with Δδ¹⁵N_C, Δδ¹⁵N_N and Δδ¹⁵N_A except for share of grassland, AirRH and UrbClim, which showed an overall negative relationship with the response variables (Fig 6A–6C; Figures A-C in S3 File). The explanatory power (R²) of the Δδ¹⁵N_G and Δδ¹⁵N_F RF models was 0.45 and 0.21, with a corresponding validation accuracy of 0.14 and 0.17. Δδ¹⁵N_G increased monotonically with RailwD in buffer areas of multiple sizes around the biotope patch and impervious surface_100 (45% and 14%) but decreased with share of grassland_5000 (13%) and similarly, Δδ¹⁵N_F RF increased with FAR, PopD, impervious surfaceand RailwD when summing up the explanatory power of the predictors for buffer areas of multiple sizes around the biotope patch (21%, 17%, 16% and 13%; Figures D-E in S3 File). Note that when considering separately the explanatory power of the predictors for different buffer areas sizes around the biotope patch, impervious surface_100 was the main predictor selected for Δδ¹⁵N_F RF model. The Δδ¹⁵N_NG, Δδ¹⁵N_NF and Δδ¹⁵N_AF RF models reflected the pattern observed for the overall plant groups with different biogeographic status, with RailwD, FAR, impervious surface and PopD in buffer areas of multiple sizes around the biotope patch being the main predictors (Figures F-H in S3 File).

Biodiversity-ecosystem functioning in urban grasslands

AGB_C increased significantly with species richness in urban grasslands (Fig 7A). Namely, AGB_C increased up to a saturating value of 350 g m^-2 at 30 species. The relationship between species richness and AGB_C was not influenced by the proportion of alien species (Fig 7A). The positive biodiversity-aboveground biomass relationship in non-manipulated urban grasslands was similar to that found in experimental non-urban grasslands across Europe, in other words, the regression slope in our study fell within the slope range reported by [36], Fig 7C. Unlike in the case of AGB_C, species richness was not significantly correlated to iWUE_C or Δδ¹⁵N_C (Fig 7B and 7D).

Species level AGB, iWUE or Δδ¹⁵N

Species richness was not significantly correlated with photosynthetic parameters, iWUE or Δδ¹⁵N at the species level. The BNI did not affect photosynthetic parameters of C. epigejos and P. lanceolata, iWUE_C.ep or iWUE_P.la but it was negatively correlated to Δδ¹⁵N_C.ep in summer. In C. epigejos, both A (r_s = 0.59, p = 0.036, n = 13) and g_s (r_s = 0.61, p = 0.028, n = 13) increased with SVF in spring. No significant correlation was found between SVF and photosynthetic parameters in C.epigejos in spring or in P.lanceolata in spring and summer. iWUE_C.ep increased significantly with SVF in both spring and summer (r_s = 0.8, p = 0.001, n = 13; r_s = 0.53, p = 0.029, n = 17). No effect of SVF on Δδ¹⁵N_C.ep or Δδ¹⁵N_P.la was found.

Discussion

The relationship between biodiversity and ecosystem function has received much attention during the last decades in ecological research (e.g. [34, 36, 39, 109]). At present, accelerated urban sprawl poses the question whether urban ecosystems are able to provide services in a comparable way to their non-urban counterparts. We focused on aboveground biomass, intrinsic water use efficiency and ¹⁵N enrichment factor (AGB, iWUE and Δδ¹⁵N) as proxies for biomass production, water and N cycling, respectively, which have been acknowledged as relevant processes for ecosystem integrity [110]. As a major insight we found that biodiversity, abiotic, and biotic novelty were related differently to these ecosystem functions. Community aboveground biomass (AGB_C) was mainly explained by species richness. In contrast, community intrinsic water use efficiency (iWUE_C) and community ¹⁵N enrichment factor (Δδ¹⁵N_C) were mostly determined by light availability (depicted by the SVF) and urban parameters, respectively, which overrode biodiversity effects. In addition, our results indicate that abiotic novelty potentially favors alien plants in Berlin, mainly by enhancing their dispersal and fitness (via increased iWUE) under low water availability conditions. We found that abiotic novelty, exemplified by urban parameters, specifically affected aboveground biomass and intrinsic water use efficiency of aliens (AGB_A and iWUE_A), but this effect did not translate into effects on the community level. However, abiotic novelty did affect both ¹⁵N enrichment factor of natives and aliens (Δδ¹⁵N_N and Δδ¹⁵N_A), suggesting a broader impact of urbanization on N cycling compared to C and water cycling. At the species level, abiotic novelty, and specifically sky view factor (SVF), appeared to be the prevailing driver of photosynthetic performance and resource use efficiency over both species richness and biotic novelty.

Key predictors for the functioning of urban grasslands

Aboveground biomass

Plant invasions have been linked to declines in species diversity and increases in primary production, often attributed to the presence of highly productive alien species [47]. In the studied urban grasslands in Berlin, the relationship between species richness and AGB_C was not influenced by the proportion of alien species, that is, lower species richness was not necessarily associated to a higher share of aliens. We found that AGB_C increased significantly with plant species richness and more importantly, according to the RF models, species richness was the strongest predictor for AGB_C. Although we did not explicitly assess complementarity and selection effects [111] in our study, the fact that the increase in AGB_C was associated with increasing Rao´s Q suggests that complementarity was the main mechanisms driving the observed diversity-ecosystem functioning relationship. Functional diversity quantifies trait variability in a given community, and it has been suggested to foster niche complementarity by reducing niche overlap among species pools comprising higher trait ranges [110]. The positive effect of size of dry grassland biotope patch on AGB_C further suggests that biomass is boosted by complementarity effects, since biotope space strengthens the biodiversity-AGB relationship by enabling higher niche complementarity among species [112].

Our results additionally indicate that light availability (depicted by the SVF) limits primary production in urban grasslands in Berlin. Moss cover was also associated to reduced AGB_C, mainly via a decrease in the aboveground biomass of graminoids (AGB_G). Numerous studies have described a negative relationship (mostly attributed to competition) between bryophyte and vascular plant biomass [113]. For instance, in the arctic tundra, moss cover was associated to lower graminoids productivity via reduced soil temperature and N availability [114]. However, our results do not provide enough information to elucidate the mechanisms behind this relationship and reduced vascular plant aboveground biomass production might have also relieved light competition for bryophytes.

Our study points to differing controlling factors for the AGB of plants with different biogeographic status. Whereas aboveground biomass of natives (AGB_N) was mainly affected by plant diversity and size of dry grassland biotope patch, aboveground biomass of aliens (AGB_A) was mainly explained by abiotic novelty-related parameters, specifically PopD (population density). We hypothesize that this discrepancy stems from different predominant dispersal mechanisms of the alien and native species most frequently found in our plots. The most common aliens were Berteroa incana (L.) DC. and Conyza canadensis (L.) Cronquist, found in 9 and 4 plots, respectively (the geographic location of these plots is provided in S4 Table). Despite the fact that these species are predominantly wind-dispersed, dispersal via footwear and mowing machines have been reported as potential dispersal mechanisms for B. incana [115]. Unintended human-mediated dispersal by vehicles [116] or clothing attachment [117] was found to effectively disperse C. canadensis, which also has highly volatile seeds to the extent that may even reach the planetary boundary layer of the atmosphere [118]. Human-mediated dispersal is known to increase with PopD due to higher density of transport networks integrated by roads or walking routes [119]. By contrast, the most common native species in our plots were the perennial grasses C. epigejos and Festuca trachyphylla (Hack.) Hack (present in 9 and 11 plots, respectively). These clonal species seem to prioritize vegetative over seed reproduction once they are established [120] which reduces the probability of human-mediated dispersal. These results are in line with those from previous studies showing that alien species diversity increased with settlement connectivity and number of inhabitants along an urban-rural gradient in the vicinity of Frankfurt (Germany), which suggests that dispersal mechanisms of alien species are associated with human mobility patterns [121].

Intrinsic water use efficiency

Unlike in the case of AGB_C, species richness was not significantly correlated to iWUE_C in the studied urban grasslands. In our study, environmental factors related to resource availability and urbanization seemed to override biodiversity effects. In particular, SVF was the strongest predictor for iWUE_C followed by soil gravimetric water content (W_C). A higher SVF implies higher incoming solar radiation, which very likely boosted non-light-saturated photosynthesis, resulting in increased iWUE_C in grasslands with a higher share of open sky. The increasing trend in iWUE_C observed with decreasing W_C most likely reflected lowered stomatal conductance to avoid water loss associated to increasing water stress [122]. While a positive biodiversity-aboveground biomass relationship seems to be relatively general across ecosystems and studies, and even robust to nutrient enrichment and drought [33], the contrasting patterns how biodiversity relates to water use reported in the literature point to a more variable relationship. For instance, species mixtures were found to be more efficient in water use (estimated via water balance) in comparison with monocultures in experimental grasslands and this was attributed to complementary [50]. By means of a soil water balance modelling approach, a more complete exploitation of soil water, ascribed to higher root density, was linked to higher diversity in experimental grasslands [56]. In contrast, stable water isotopologue labeling revealed no positive effects of species richness on water uptake by plant communities in the temperate experimental grasslands of the Jena Experiment [123].

The differences in iWUE among functional groups is consistent with previous observations in natural grasslands, in which δ¹³C_G was lowest, followed by δ¹³C_F and δ¹³C_L ley [124, 125]. Higher stomatal conductance—which might be a reason for lower iWUE—in forbs compared to grasses was observed in soil-plant monoliths from the Jena experiment and attributed to water uptake from deeper and thus moister soil layers [126]. Even though plants with distinct biogeographic status did not significantly differ with regard to their iWUE, different dominant drivers were identified for natives and aliens. Intrinsic water use efficiency of natives (iWUE_N) was mainly controlled by light and to a lesser extent by water availability (i.e. SVF and W_C). However, urbanization, and specifically railway density (RailwD) in a buffer area of 5 km, was the main factor modulating iWUE_A. Railway embankments functioned as active water sources to nearby ecosystems in a Mediterranean watershed by increasing runoff [127]. This increased water supply could be one reason for aliens spreading along railways [128] and increased iWUE_A with decreasing RailwD might point to lower stomatal conductance associated with increasing water stress in areas with less dense railway networks. [129] indicated that aliens often have competitive physiological traits that increase their fitness in unfavorable environments, including the capacity to increase WUE in water-poor environments. This might at least partly explain the alien-specific response observed in urban grasslands in Berlin, where overall, W_C was fairly low.

¹⁵N enrichment factor

As in the case of iWUE_C, no effect of species richness on Δδ¹⁵N_C was observed in the studied urban grasslands. This result contrasts with findings in central European semi-natural grasslands, where increasing plant diversity was associated to lower Δδ¹⁵N at the community and species level, indicating a more closed N cycle that was attributed to efficient N uptake via complementarity [57]. In our study, urban parameters were consistently identified as the main controlling factors for Δδ¹⁵N at the community, biogeographic status class and functional group level, which illustrates the predominant role of abiotic novelty over biodiversity or biotic novelty in regulating N cycling in urban grasslands in Berlin. Even though the relative importance of different urban parameters for different levels varied, all predictors were associated to an increasing trend in Δδ¹⁵N, that is, they favored a more open N cycle. Very likely, this common trend reflects the association of PopD, RailwD, FAR and impervious surface with urban N emissions. Increased anthropogenic N inputs stimulate N transformations (i.e. mineralization and nitrification), ultimately favoring N mobilization and loss in the form of leachable nitrate or nitrogenous gas [130], which might have led to a more N open cycle.

Biodiversity-ecosystem functioning in urban grasslands

Our results on the positive relationship between total species richness and biomass illustrate the ability of urban grasslands to provide services in a comparable way to their non-urban counterparts. Previous studies have revealed the importance of urban grasslands for biodiversity conservation [24, 131] and cultural ecosystem services [132, 133]. Our findings highlight the relevance of biodiverse urban grasslands to ensure the provision of certain supporting and provisioning ecosystem services (i.e. primary production, carbon sequestration). In particular, we found that the positive relationship between species richness and AGB_C extensively described for European experimental grasslands [36] is maintained in the urban grasslands in Berlin even though they are intensively influenced by biotic and abiotic novelty (e.g. urbanization-related parameters) (Table 1). Our results agree with previous studies which reported a positive relationship ranging from linear to log linear with regression slopes that did not significantly differ from zero ([37] and references therein).

Species level AGB, iWUE or Δδ¹⁵N

At the species level, we found no effect of species richness on photosynthetic performance, iWUE or Δδ¹⁵N, which suggests that the diversity of the surrounding vegetation did neither enhance nor hamper AGB production, water and N use of C. epigejos or P. lanceolata in the studied urban grasslands. Similar to our study, plant species richness did not affect ¹³C abundances in AGB of individual species in experimental grasslands, where photosynthetic assimilation seemed to be rather determined by the species position (i.e. height) in the plant canopy and the associated availability of light, air humidity and CO₂ concentrations [134]. In accordance with the latter study, iWUE significantly increased with light availability (depicted by the SVF) in C. epigejos both in spring and summer, and at least in spring, this increase seemed to be driven by enhanced A and g_s. The lack of SVF effects on P. lanceolata might be due to its comparably limited access to direct sunlight due to the lower position of this rosette species in the canopy. Height most likely also determined the contrasting seasonal variations observed in iWUE for the two studied species. iWUE_C.ep increased over the growing season, possibly due to increasing height and therefore greater access to light, which boosted A. However, in P. lanceolata, increasing height of the grassland canopy probably resulted in stronger shading of the rosette leaves, which increased g_s leading to a reduction in iWUE throughout the season. The significant decrease in F_v/F_m and ΔF/F_m′ in P. lanceolata from spring to summer suggests a decline in photosynthetic capacity, in agreement with our above proposed explanation that the rate of AGB production decreased over the growing season due to increased shading. This is further supported by the results of the instant light-curves, which indicated an important decrease in intrinsic photosynthetic capacity (ETR_max, ΔF/F_m′ and ΔF/F_m′ at PPDF_sat) in P. lanceolata over the growing season but minimal seasonal variations in C. epigejos. The latter results, along with the lack of consistent effects of the degree of biotic novelty of the plots on photosynthetic parameters, iWUE or Δδ¹⁵N, suggest that abiotic novelty is the prevailing driver of photosynthetic performance and resource use efficiency of the selected individual species.

Conclusions

While the biodiversity-ecosystem functioning relationship is critically understudied for urban grasslands, this study untangles a range of relationships between plant species richness, abiotic and biotic novelty and three important proxies of ecosystem functioning, namely aboveground biomass, intrinsic water use efficiency and ¹⁵N enrichment factor. Even though we identified a significant impact of abiotic novelty on AGB, iWUE and Δδ¹⁵N for various plant groups, the species richness-AGB_C relationship found in the studied urban grasslands in Berlin was comparable to that described in non-urban experimental grasslands in Europe. Our results support previous evidence that conserving and enhancing biodiversity in urban ecosystems is essential to warrant certain functions and ultimately their associated services (i.e. AGB production and primary production or carbon sequestration) and stress that management of urban grasslands in Berlin and other large cities should aim at preserving plant species richness. However, our results also suggest that preserving species richness is insufficient to ensure the provision of ecosystem services associated to other functions such as water and N use, which highlights the need for further studies addressing urban ecosystem functioning and exploring potential additional management measures directed to modify urban abiotic parameters.

Supporting information

S1 Appendix [docx]
Calculation of the Biotic Novelty Index (BNI).

S1 Table [grey]

S2 Table [docx]
Overview of the random forests run for aboveground biomass (AGB), intrinsic water use efficiency (iWUE) and N enrichment factor (ΔδN) at the community, plant group (either biogeographic status, functional or combined) and species level in Berlin grasslands in 2017.

S3 Table [docx]
Species present in the sampled quadrats in the grassland plots in Berlin in late summer 2017 along with the number and percentage of plots in which they were present.

S4 Table [docx]
Geographic location of the grassland plots in Berlin in which and were present in late summer 2017.

S1 File [pdf]
Partial dependence plots for the predictors selected by random forest models for aboveground biomass.

S2 File [pdf]
Partial dependence plots for the predictors selected by random forest models for intrinsic water use efficiency.

S3 File [pdf]
Partial dependence plots for the predictors selected by random forest models for N enrichment factor.

S1 Dataset [xlsx]
Digitized data extracted from Fig 3 in [].