Human pressure drives biodiversity–multifunctionality relationships in large Neotropical wetlands

Many studies have shown that biodiversity regulates multiple ecological functions that are needed to maintain the productivity of a variety of ecosystem types. What is unknown is how human activities may alter the ‘multifunctionality’ of ecosystems through both direct impacts on ecosystems and indirect effects mediated by the loss of multifaceted biodiversity. Using an extensive database of 72 lakes spanning four large Neotropical wetlands in Brazil, we demonstrate that species richness and functional diversity across multiple larger (fish and macrophytes) and smaller (microcrustaceans, rotifers, protists and phytoplankton) groups of aquatic organisms are positively associated with ecosystem multifunctionality. Whereas the positive association between smaller organisms and multifunctionality broke down with increasing human pressure, this positive relationship was maintained for larger organisms despite the increase in human pressure. Human pressure impacted multifunctionality both directly and indirectly through reducing species richness and functional diversity of multiple organismal groups. These findings provide further empirical evidence about the importance of aquatic biodiversity for maintaining wetland multifunctionality. Despite the key role of biodiversity, human pressure reduces the diversity of multiple groups of aquatic organisms, eroding their positive impacts on a suite of ecological functions that sustain wetlands. Analysis of the species richness and functional diversity among species across 72 lakes finds that both variables are positively associated with ecosystem multifunctionality, but that—for smaller organisms only—these positive relationships break down with increasing human pressure.

H uman activities are causing biodiversity to decline worldwide 1,2 , which has led to an interest in how biodiversity loss might alter the functioning of ecosystems 3 . Most studies have revealed positive and saturating effects of biodiversity on single ecosystem functions 4 . Empirical evidence suggests that species are ecologically unique and can play complementary roles in natural systems, thus varying in their contributions to different functions [3][4][5] . As a consequence, the effect of biodiversity on ecosystem functioning is stronger-and the relationship is non-saturating-when multiple functions are considered (hereafter 'multifunctionality') [5][6][7][8] . Therefore, it has been increasingly recognized that biodiversity and multifunctionality are strongly associated. This recognition has led to the prediction that as biodiversity declines in human-dominated ecosystems, their ability to sustain multiple ecosystem functions is impaired, ultimately altering the biodiversity-multifunctionality relationship 3,[9][10][11][12][13] . Current evidence supporting the anthropogenic impacts on biodiversity-multifunctionality relationships are scarce and comes mostly from experimental manipulations of single trophic levels [10][11][12][13] . It is possible that these studies under-estimate human impacts on biodiversity and ecosystem multifunctionality because natural systems are comprised of multiple organismal groups of varying trophic levels, and different trophic levels may combine to have stronger impacts on multifunctionality [5][6][7] . Further research applying a multi-trophic perspective is needed to develop a more mechanistic understanding of the consequences of human pressures for biodiversity-ecosystem functioning relationships in natural systems worldwide.
Here we used a unique dataset from 72 lakes distributed across four large Neotropical wetlands of Brazil (Amazon, Araguaia, Pantanal and Paraná) to test how the cumulative effect of multiple human pressures impacts the relationship between biodiversity and multifunctionality. These four wetlands provide a unique opportunity to test the influence of human pressures across broad spatial scales as the lakes span a 3,700,000 km 2 gradient of distinct human activities (Fig. 1). We quantified human pressure on the wetland using the Human Footprint (HFP) index 14 , which was extracted for each lake individually (Methods). The HFP is a recently developed index that incorporates eight different human pressures: (1) built environments, (2) crop land, (3) pasture land, (4) human density, (5) nighttime lights, (6) railways, (7) roads and (8) navigable waterways into a standardized cumulative index of human pressure 14 . This index provides an interesting opportunity to understand how human pressures are affecting biodiversity-multifunctionality relationships in natural to human-dominated systems.
We compiled data on the species richness and functional diversity of seven taxonomic groups, including fish, aquatic macrophytes, microcrustaceans, rotifers, phytoplankton, ciliates and testate amoebae. These data comprised 1,465 plant, animal and microbial species. Because biodiversity-multifunctionality relationships can be multidimensional 6,7 , we also used measures of multi-diversity (joint diversity of all organismal groups, both for species richness and functional diversity) 15 . Studies considering multi-diversity have found strong biodiversity-multifunctionality relationships [6][7][8] . To estimate functional diversity, we focused on a core set of independent organismal traits that mediate the species response to human pressures (Supplementary Table 1): body size, resource-use (for example, feeding groups, growth forms and mixotrophy) and mobility (for example, migration ability, propagation method and cell motility) traits. These traits are often linked to multiple ecosystem functions in wetlands. For instance, body size, feeding groups and migration ability are related to metabolism, multi-trophic biomass production and nutrient cycling 16,17 . We further quantified ecosystem multifunctionality by using a set of 11 variables that included nutrient concentrations (in situ measurements of N and P water concentrations), metabolism (daily changes in water O 2 concentration), biomass at multiple trophic levels (algae, herbivores, carnivores, detritivores and omnivores), microorganism abundance (bacterial cell densities), availability of photosynthetically active radiation (light availability underwater) and variation in habitat complexity underwater (variation in plant above bottom cover). Together, these variables measure environmental characteristics that are directly linked to ecosystem functions. A detailed rationale for each variable is provided in Supplementary Table 2. We quantified multifunctionality using three common approaches: (1) the averaging multifunctionality index, (2) the multi-threshold multifunctionality index and (3) multiple single functions. The averaging approach takes the average of the standardized values of each single function. In contrast, the multi-threshold considers the number of functions that simultaneously surpass a range of thresholds, which are expressed as a percentage of the highest observed level of functioning (here 1-99%). These three approaches are complementary, and when taken together, they provide a robust estimation of how multiple functions (averaging and multi-pillar approach) and single functions respond to biodiversity enhancement [5][6][7][8]18 . Because no studies have examined the broad-scale relationships between biodiversity and ecosystem multifunctionality across wetlands, we first established whether species richness and the functional diversity of the seven organismal groups were, in fact, related to multifunctionality as previous narrow-scale evidence suggests 17,19 . For this, we employed multiple linear mixed models considering species richness and functional diversity as predictors and multifunctionality as the response. After confirming a consistent relationship, we also used linear mixed models to determine how human pressures alter these biodiversity-multifunctionality relationships. Lastly, we used structural equation models (SEMs) to investigate the direct and indirect biodiversity-mediated pathways by which human pressure can influence multifunctionality in wetlands.

Results and discussion
Across four hyperdiverse Neotropical wetlands, we found significant positive relationships between the diversity of single groups of aquatic organisms and the multi-diversity of all groups with ecosystem multifunctionality (Figs. 2 and 3 and Supplementary Table 3). This finding was consistent for both species richness and functional diversity (Figs. 2 and 3). Our model averaging procedure revealed that the biodiversity of organismal groups was the best predictor of multifunctionality, even after accounting for influence of other well-known drivers of multifunctionality such as space, climate (precipitation and temperature) and aquatic properties (conductivity, pH and water level; Supplementary Table 4). The positive association between aquatic biodiversity and multifunctionality persisted regardless of how the measures of multifunctionality were weighted (Supplementary Figs. 1 and 2). The multi-threshold approach provided additional evidence showing that the mean minimum threshold at which the species richness of organismal groups had its strongest effects on multifunctionality averaged 57% (range 5-92%, Supplementary Fig. 3). Similarly, the mean minimum threshold at which functional diversity had its strongest effects on multifunctionality was 91% (range 70-99%, Supplementary Fig. 4). The diversity of aquatic organism groups was also positively associated with most of the individual ecosystem functions, although each organismal group was more closely associated with specific functions (Supplementary Tables 5 and 6). Here fish diversity was strongly related to multi-trophic biomass, macrophyte diversity was most The linear association between multifunctionality and the species richness of the seven selected taxonomic groups and the composite metric of their joint richness (multi-diversity, standardized between 0 and 1) 15 in four Neotropical wetlands; n = 72 lakes. Statistical analysis was performed using linear mixed-effect models. Dashed black and solid lines are predicted (fitted) values from LMMs for overall and local trends (for each wetland ecosystem), respectively. Shaded areas show the 95% confidence interval for the overall trend. R 2 = marginal (that is, variance of the fixed effects). *P < 0.05, **P < 0.01, ***P < 0.01. The richness of microcrustaceans, testate amoebae and phytoplankton was log-transformed before the analysis. Full model results are provided in Supplementary Table 3. Multifunctionality is represented by the averaging index, which reflects changes in the average level of the 11 ecosystem functions. Very high averaging index levels (close to 1) mean that all functions reach their maximum level of performance simultaneously. By contrast, the lowest values (close to 0) mean all functions are at their minimum level of performance. Illustration credit: João Vitor Fonseca da Silva (https://gqromero.wixsite.com/lab/team-3).
strongly related to light availability and habitat complexity, whereas microorganism diversity was most related to nutrient concentrations and ecosystem metabolism (Extended Data Figs. 1 and 2). Finally, aquatic biodiversity had stronger effects on multifunctionality than other multifunctionality drivers (Extended Data Figs. 3 and 4; SEM: total effect of composite species richness on multifunctionality 0.79, total effect of composite functional diversity on multifunctionality 0.72). Collectively, our broad-scale dataset revealed strong and consistent associations between the diversity of multiple groups of aquatic organisms and ecosystem multifunctionality. These results underline the important role of multiple elements of biodiversity in driving ecosystem functioning in Neotropical wetlands 15,16,18 , as in other ecosystem types such as dry land 8 and forest 7 . The close association between biodiversity and multifunctionality suggests that biodiversity loss might impact the ability of wetlands to maintain their functioning [4][5][6][7][8] . Analysis of the relationship between HFP and biodiversity revealed a decline in species richness and functional diversity with increasing HFP (Supplementary Figs. 5 and 6). To test how this affected the relationship between biodiversity and multifunctionality, we examined how interaction of HFP × biodiversity influenced the slope of biodiversity-multifunctionality relationships. While the isolated effect of species richness on multifunctionality was positive for most organismal groups, the interactive HFP × species richness effect was negative (Fig. 4a). Similarly, the isolated effect of functional diversity on multifunctionality was positive, but the interactive HFP × functional diversity effect was strongly negative (Fig. 5a). This suggests that human pressure can alter the relationship of both species' richness and functional diversity with multifunctionality. By decomposing the effect of biodiversity on multifunctionality through low, medium and high HFP intensity, we found that the positive effect of species richness and functional diversity on multifunctionality declined from low to high HFP intensity ( Fig. 4b and Fig. 5b). In particular, the effect of the diversity of smaller organisms (such as, microcrustaceans, testate amoebae, ciliates and rotifera) on multifunctionality shifted from positive at low HFP intensity to neutral or negative at high HFP intensity (Figs. 4 and 5). By contrast, the positive effect of the diversity of larger organisms (such as fish and macrophytes) on multifunctionality was maintained despite increased HFP. These results illustrate how the ability of smaller organisms to promote multifunctionality is sensitive to human pressure and simultaneously highlight the importance of larger organisms for maintaining ecosystem functioning in a human-dominated world 20 .
The changes in the magnitude and direction of the relations between biodiversity and multifunctionality suggest that such relationships can be context dependent in wetlands 21 . This is more evident for smaller groups of aquatic organisms as their effects on multifunctionality changed from positive at low HFP intensity to negative at high HFP intensity. Using a structural equation model, we disentangled the direct and biodiversity-mediated, indirect pathways by which human pressures affect multifunctionality.

Fig. 3 | Relationship between the functional diversity of aquatic organisms and multifunctionality in Neotropical wetlands.
The linear association between multifunctionality and the functional diversity of the seven selected taxonomic groups and the composite metric of their joint functional diversity (multi-diversity; standardized between 0 and 1) 15 in four Neotropical wetlands; n = 72 lakes. Statistical analysis was performed using linear mixed-effect models. Dashed black and solid lines are predicted (fitted) values from LMMs for overall and local trends (for each wetland ecosystem), respectively. Shaded areas show the 95% confidence interval for the overall trend. R 2 = marginal (that is, variance of the fixed effects). *P < 0.05, **P < 0.01, ***P < 0.01. The richness of microcrustaceans, testate amoebae and phytoplankton was log-transformed before the analysis. Full model results are provided in Supplementary  Table 3. Multifunctionality is represented by the averaging index, which reflects changes in the average level of the 11 ecosystem functions. Very high averaging index levels (close to 1) mean that all functions reach their maximum level of performance simultaneously. By contrast, the lowest values (close to 0) mean all functions are at their minimum level of performance. Illustration credit: João Vitor Fonseca da Silva (https://gqromero.wixsite.com/lab/team-3).
We demonstrate that the direct effect of HFP on multifunctionality was consistently negative across all wetlands ( Fig. 6a and Supplementary Tables 8-10). This is consistent with the fact that the studied wetlands cover regions with intensive human activities ( Fig. 1). Most of the studied wetlands cover areas of simultaneous crops of soy and sugar cane and pasture lands grazed by cattle [22][23][24][25] , and Paraná wetland is located downstream of one of the most populated areas on the planet 22 . Consequently, multiple human pressures can jointly affect the integrity of these wetlands by decreasing biodiversity and ecosystem multifunctionality ( Supplementary Fig. 7).
Beyond the direct negative effect on multifunctionality, HFP had large indirect negative effects on the multifunctionality mediated by declining species richness and functional diversity (Fig. 6). Although indirect negative effects of human pressure were driven by the decline in the diversity of most organismal groups, these effects were strongly mediated by fish diversity (Fig. 6b,d). This is consistent with the fact that fish diversity has the greatest influence on functioning of wetlands 16,17 , and loss in fish diversity is known to impact multiple ecosystem functions 26 . The negative indirect biodiversity-mediated effects of human pressure on multifunctionality were also consistent across wetlands (Supplementary Table 11).
Combined with the fact that the positive effects of biodiversity on multifunctionality decreased with increasing HFP (Fig. 4), our results highlight that if the human pressures continue to increase 27 , preservation of biodiversity for maintaining multifunctionality will not be sufficient unless they are accompanied by a reduction of human pressures. Seen in the light of the increasing human influence on natural landscapes, our results illustrate the importance of considering multiple pathways through which human pressures can influence ecosystem multifunctionality.

Conclusion
We have provided empirical evidence of a positive broad-scale relationship between the diversity of multiple groups of aquatic organisms and the multifunctionality of wetland ecosystems. We demonstrate that a positive association between aquatic biodiversity and multifunctionality occurs for both single metrics of diversity and for those combined into multi-diversity. These positive relationships are also apparent for the seven groups of aquatic organisms, although larger organisms are more strongly linked to multifunctionality than smaller organisms. Collectively, our findings highlight the importance of aquatic biodiversity for maintaining ecosystem multifunctionality and their associated services 28 . It is imperative that biodiversity conservation be a key management priority in wetlands 29 and that ecosystem management targets the joint conservation of multiple components of aquatic biodiversity, from vertebrates to plants and microorganisms. We have also shown that human pressures degrade the positive relationship between biodiversity and multifunctionality, which occur both directly and indirectly as human pressures reduce the biodiversity needed to maintain numerous ecosystem functions. These findings demonstrate that human pressures are degrading multifunctionality through multiple pathways. Consequently, conserving the functioning of wetlands will be a major challenge as human pressures continue to increase in these ecosystems worldwide 29,30 . More broadly, reducing human pressures must be addressed urgently in wetlands  48 , including information about the diversity of seven taxonomic groups of aquatic organisms (Methods). We accounted for multiple ecosystem drivers, including distance from the equator, climate (temperature and precipitation) and aquatic properties (pH, conductivity and water level). We grouped the different categories of drivers (climate, space and water properties) into the same box for graphic simplicity; nevertheless, it does not represent latent variables. Solid black and dashed grey arrows represent significant pathways (P ≤ 0.05) and non-significant pathways (P ≥ 0.05), respectively. The thickness of the significant pathways (arrows) represents the magnitude of the standardized regression coefficient. Numbers adjacent to arrows are the standardized effect size. R 2 for component models are given in Supplementary Table 12. Significance levels are *P < 0.05, **P < 0.01, ***P < 0.001. For simplicity, we grouped the effects of ecosystem drivers (distance, HFP, climate and water properties) on the diversity of each of the seven taxonomic group in boxes. Specifically, Box A represents the effect of distance from the equator, Box B the effect of HFP, Box C the effect of climate and Box D the effect of water properties. Full model outputs and information about boxes A-D are provided in Supplementary Tables 8 and 9. CFI, comparative fit index; RMSEA, the root mean square error of approximation; MAP, mean annual precipitation; MAT, mean annual temperature b,d, The standardized indirect effects of the human footprint on multifunctionality mediated by species richness (b) and the functional diversity (d) of each organismal group used to compute the composite diversity (Supplementary Table 11). Illustration credit: João Vitor Fonseca da Silva (https://gqromero.wixsite.com/lab/team-3).
as these systems rank among the most diverse and productive ecosystems globally, providing a suite of functions and services essential for human well being.

Methods
Study sites and data collection. The study comprised the four largest South American wetlands-Amazon, Araguaia, Pantanal and Paraná-encompassing a subcontinental spatial area of approximately 3,700,000 km 2 and 72 lake ecosystems (Fig. 1). These wetlands are subject to distinct intensities of human pressure. Amazon is a global biodiversity hotspot and is more preserved than Araguaia and Pantanal that are both subject to moderate human pressure (Fig. 1). Paraná includes 150 constructed dams 31 and faces the strongest human pressure among the four wetlands. The climate ranges from subtropical to tropical, with a mean annual temperature of 16-29 °C and a mean precipitation of 1,300-2,000 mm yr −1 (ref. 32 ). The field data were collected between August 2011 and May 2012. The wetland lakes were surveyed under the Brazilian programme 'National System for Research in Biodiversity' (Sisbiota Brazil). The field surveys were designed to include lakes representing a wide range of climate, human pressure and environmental conditions. They followed a standardized sampling protocol, and the sampling effort was the same in all lakes 32 . To provide a comprehensive assessment of aquatic communities, we performed one sampling during the dry season and another during the wet season in each lake. The sampling included fish, aquatic macrophytes, microcrustaceans (cladocerans and copepods), rotifers, phytoplankton, testate amoebae and ciliates (Supplementary Methods).
Diversity measure. We quantified the species richness of the seven taxonomic groups in all 72 lakes. After identifying each individual to species level, we determined 325 fish species, 87 macrophyte species, 99 microcrustacean species, 124 rotifer species, 598 phytoplankton species, 124 testate amoebae species and 108 ciliate species. Sample coverage was equal for all wetlands, and we calculated estimated species richness as the Chao index with abundance-based data using the R package 'iNEXT' 33 , which is based on rarefaction and extrapolation of Hill numbers and provides an unbiased estimate of asymptotic species richness and enables comparisons among wetlands with different numbers of individuals. We used the Chao species richness because richness is the most commonly used and simplest metric of biodiversity [5][6][7][8] . We also measured the key functional traits for all organismal groups. We focused on the traits that are known to govern the patterns of spatial distribution and individual fitness and which also influence ecosystem processes (Supplementary Table 1) 16,34 . These traits fall into the three broad categories: (1) body size (maximum body length for animal taxa or cell volume for phytoplankton), (2) resource and habitat use traits (feeding groups for animal taxa, growth form for macrophytes, nitrogen fixation or mixotrophy for microorganisms) and (3) mobility traits (dispersal ability for animal taxa, propagation means for macrophytes and cell motility for microorganisms).
To determine the functional diversity of each organismal group, we calculated functional dispersion-that is, the mean distance in multidimensional trait space of the individual species to the centroid of all species 35 . This measure provides a robust estimate of functional diversity. Because the relationship between biodiversity and ecosystem functioning can be multidimensional on both the predictor (biodiversity) and response side (multifunctionality) 6,7 , we estimated a multi-diversity index including the diversity of the seven organismal groups 15 . We first standardized the diversity values of each organismal group between 0 and 1 (species richness or functional diversity) by scaling them to the maximum observed value, and then we average these standardized diversity values 15 . This procedure ensures that the diversity of each organism group contributes equally to the multi-diversity. We calculated separately the multi-diversity index for species richness and functional diversity. The multi-diversity index has been widely used because it reflects very well the biodiversity-multifunctionality relationships in multi-trophic ecosystems 8,11,15,17 .
Assessing ecosystem functions and properties. In each lake, 11 ecosystem variables regulated by aquatic organisms and belonging to a wide range of ecosystem functions and properties were measured (Supplementary Table 2). These functions and properties included: (1) nutrient concentrations represented by in situ measurements of total phosphorous (mg l −1 ) and total nitrogen (mg l −1 ) available in the water. Total phosphorus and nitrogen cover all fractions of these nutrients, including nitrate, nitrite, ammonia, particulate phosphate, dissolved organic phosphate and orthophosphate. We took water samples in each lake and in the laboratory; nitrogen was quantified according to Mackereth et al 36 ., while phosphorus was quantified following ref. 37 . (2) Ecosystem metabolism represented by the daily variation of dissolved oxygen in the water (mg l −1 d −1 ), which was measured from dawn to dusk in each lake using a digital oximeter portable YSI aid (Digimed). We use the mean of daily oxygen variation as it represents the change in the metabolic underwater regime 38 . (3) Multi-trophic standing biomass was represented by the biomass of algae, carnivorous fish, omnivorous fish, herbivorous fish and detritivorous fish. Algae standing biomass was quantified using biovolume (individuals per mm l −1 ) of identified algae species. Biovolume was estimated by multiplying the abundance of each species by its mean volume 39 .
Fish were classified into trophic groups using information from feeding trials and gut content analysis 16,32 . Afterwards, the fish counts within each trophic group were converted to biomass (g m −2 ) using published species-specific length-weight relationships 40 . (4) Availability of photosynthetically active radiation represented by light availability underwater (m). We quantified light availability underwater by the depth of the euphotic zone, which represents the depth (m) of the lake where there is sufficient light incidence for autotrophs. The euphotic zone was calculated as Secchi depth multiplied by 1.7, where 1.7 is a correction factor for estimating the light available underwater 32 . (5) Microorganism abundance (cells per ml) was quantified using bacterial abundance. To record the accumulative abundance of bacteria, we took water samples at the subsurface (approximately 30 cm below the air-water interface) at the central, deepest region of each lake using polyethylene flasks. Bacteria were analysed from water samples treated with a fixative solution composed of alkaline Lugol's solution, borate buffered formalin and sodium thiosulfate that was filtered through black Nuclepore filters (0.2 μm and 0.8 μm, respectively) and stained with fluorochrome DAPI (4,6-diamidino-2-phenylindole) 41 . Bacterial quantification was done with an epifluorescence microscope at a magnification of ×1,000 (Olympus BX51). (6) Variation of underwater habitat complexity was quantified based on variations in the aboveground cover of aquatic plants (m −2 ). We estimated the area of all leaves and culm of each plant species. We then summed the area of all leaves and culm to obtain the aboveground area cover by each individual. We calculated the standard deviation of the aboveground area cover between all plant species and used this standard deviation as a proxy of variation in the aboveground vegetal cover.

Pairwise correlation between ecosystem functions.
To assess the potential for a trade-off between individual ecosystem characteristics, we calculated Pearson correlation coefficients between each pair of individual standardized functions. Of the possible 45 combinations of pairwise functions, we found only seven strong correlations (r = 0.5; Supplementary Fig. 8). To remove any bias in our multifunctionality index, the highly correlated functions were down weighted in its calculation ( Supplementary Fig. 9), as described in Manning et al. 42 . Ecosystem functions were grouped into clusters according to their correlations. This weighted approach indicated three different clear clusters: (1) aboveground plant cover, (2) available N and P, light availability underwater, daily oxygen variation and algal biomass and (3) carnivore biomass, omnivore biomass, detritivore biomass, omnivore biomass and bacterial abundance. Weighted multifunctionality was then calculated as the average of all variables within each cluster. For instance, each function within cluster 2 was weighted with a weight of 0.2. These functions were then averaged into a standardized variable. We repeated the analyses of the relationship between biodiversity and multifunctionality for the weighted multifunctionality to determine whether the results differed between weighted and non-weighted multifunctionality (see ref. 42 ).

Assessing ecosystem multifunctionality.
To obtain robust and quantitative multifunctionality indexes, we used three multifunctionality approaches: (1) the averaging multifunctionality index, (2) the multi-threshold multifunctionality index and (3) the multiple single functions index 18 . To obtain the averaging ecosystem multifunctionality index, we standardized all 11 ecosystem functions between 0 and 1 (rawFunction − min(rawFunction) / (max(rawFunction) − min(rawFunction)) and then calculated their means. The averaging ecosystem multifunctionality index is the most commonly used index in the multifunctionality literature 5,18 but has the limitations in that the number of functions with high performance are impossible to obtain and it does not allow for potential trade-offs between functions. To take these limitations into account, we used the multi-threshold index. This index calculates how many functions simultaneously exceed a predefined percentage of the maximum observed value of each individual function. Because the selection of any threshold is arbitrary, analysing multiple thresholds of maximum functioning is recommended 18 . We analysed the effect of the diversity of each organismal group on multifunctionality across the full range of thresholds from 1% to 99%. We used the mean of the three largest values of each ecosystem variable across all lakes as the observed maximum to reduce the impact of potential outliers.
Assessing the HFP on wetlands. We used the global HFP map as a surrogate of the cumulative human-induced pressure on the wetlands 14 . This map is constructed from an ensemble of eight human pressures: (1) the extent of built environments, (2) crop land, (3) pasture land, (4) human population density, (5) nighttime lights, (6) railways, (7) roads and (8) navigable waterways. To facilitate comparison among pressures, each pressure was weighted (details on the weightings are provided below). The pressures were weighted according to their relative intensity 14 . For example, (1) constructed environments are areas related to urban settlements such as buildings and urban areas. The pressure of built environments was assigned a score of 10 (that is, a score of 10 is assigned if there are built environments; otherwise a score of 0 is assigned). (2) Crop land is characterized by monocultures with high inputs of pesticides and fertilizers. In terms of HFP, the crop land pressures received a score between 0 and 7, where 7 indicates intensive agriculture, and 0 indicates the absence of crop lands. (3) Pasture land includes some of the major land uses worldwide and is characterized by cattle and sheep farming. The pressure of pastures on wetlands was assigned a score of 4, which was scaled from 0 to 4 using the percent of pasture for each 1 km 2 pixel. (4) Human population is an important underlying driver of the global change of natural ecosystems. Human density was mapped using gridded population downscaled to match the 1 km 2 resolution. All areas with a population above 1,000 people km −2 were assigned a pressure score of 10. For less populated areas, the pressure score is logarithmically scaled using the following estimation: pressure score = 3.333 × log (population density + 1). (5) Nighttime lights include electric infrastructure related to more rural areas that are not part of built environments. To calculate the pressure of nighttime lights, the areas were divided into ten quantiles of increased nighttime light intensity associated with scores between 1 and 10, while areas with no lights were assigned a zero score. (6) Railways are essential human infrastructures that influence natural ecosystems. The direct pressure of railways was assigned a score of 8 for a distance of 0.5 km on either side of the railway. (7) Likewise, roads modify the landscape where they are built. The direct and indirect pressure of roads on wetlands was assigned a score of 8 for 0.5 km (direct impact), while nearby areas up to 15 km received a score value that decayed exponentially on either side of the road (indirect impact). (8) Navigable waterways act as conduits for people to access nature, resulting in impacts on wetlands. The pressure of navigable waterways was assigned a score of 4, which decayed exponentially out to 15 km away from the water banks. For full details of HFP estimation, see refs. 2,14 .
The average HFP of the 1 km 2 pixels (cell-size resolution) overlapping each lake was extracted to derive the cumulative pressure, and this average HFP ranged between 0 and 50 (cumulative sum of all individual human pressures). The average HFP was extracted using the 'raster' R package 43 through a global HFP map that was available for the year 2009. The eight human pressures are not mutually exclusive and may co-occur in the same wetland or among and within wetlands. The HFP was initially developed to represent human pressures in terrestrial systems 14 , but most of these human pressures extensively affect wetland ecosystems. For instance, Brazil has experienced rapid expansion of urban areas 44 . Along with the increase in human populations in the vicinity of wetlands, there has been an increased pressure on these ecosystems from sewage, cattle and sheep pastures, railways, roads and navigable waterways 45 . We found negative correlations between the individual pressures with biodiversity and multifunctionality, which suggest that the use of the HFP in our study is robust (Supplementary Fig. 7).

Statistical analyses.
Linking aquatic biodiversity to multifunctionality. To determine the direct link between aquatic biodiversity and average multifunctionality across four wetland ecosystems, we fitted a series of linear mixed effects models (LMMs) to the surveyed data. We tested the relationship of (1) species richness and (2) functional diversity of single organismal groups and (3) multi-diversity with the ecosystem multifunctionality. The models were run using the function lme in the 'nlme' package 46 . We included wetlands and two sampling periods as our random structure and allowed the intercept and slopes to vary by wetland. The assumptions of normality, linearity and homoscedasticity were verified using graphical diagnostics (QQ plots and residual plots). To determine the importance of other biotic and abiotic variables besides biodiversity for multifunctionality, we included other well-known drivers of multifunctionality, such as space (distance from equator), climate (temperature and mean annual precipitation) and aquatic properties (pH, conductivity and water level; Supplementary Methods). We performed a model averaging procedure that calculated all possible subset models and chose from this set those subset models with the lowest values (ΔAICc ≤ 2) of the Akaike Information Criterion corrected for small sample size (AICc). This analysis was conducted using the R package 'MuMIn' 47 .
Using LMMs, we also assessed the relationship of species richness and functional diversity of single organismal groups and multi-diversity with each of the 11 individual ecosystem functions. This allowed us to compare the multifunctionality results to the performance of individual functions. Before these analyses, we standardized all individual ecosystem functions (z-scored: mean centred and divided by the SD (standard deviation)) to better meet model assumptions. Even so, for some functions, the residuals were highly heteroscedastic. We then modelled the variance using the function varIdent, with diversity nested by wetlands as the stratum. We considered quadratic terms for some ecosystem functions to evaluate potential nonlinear relationships.
We also modelled aquatic diversity against the number of functions above a threshold using generalized LMMs, assuming a Gaussian error distribution in the 'MASS' package 48 . Because we wanted to know whether the relationships between species richness and functional diversity with ecosystem multifunctionality varied as a function of organismal group and among the four wetlands, we fitted the generalized LMM individually to each organismal group. We then extracted and plotted the linear coefficient (fitted values) of the relationship between biodiversity and each threshold level (1% to 99%; 99 thresholds) to each wetland system. This led us to examine changes in the shape of the fitted curve for each wetland at multiple thresholds.
Effect of human pressure on biodiversity-multifunctionality relationships. We conducted LMMs between HFP and biodiversity (species richness and functional diversity of single organismal groups and multi-diversity). We found strong negative effects of HFP on biodiversity ( Supplementary Figs. 5 and 6), allowing us to determine whether HFP altered the relationship between biodiversity and ecosystem multifunctionality. We then added interaction terms for HFP × species richness and HFP × functional diversity of each single organismal group and multi-diversity to the mixed effects models and measured the estimated coefficients of these interactions on ecosystem multifunctionality. Because biodiversity and HFP are both continuous variables, analysis of their interactions could result in an interaction predictor that is collinear with the main effect 49 . Thus, we centred these variables by subtracting the sample mean from all input variable values. The mean of the centred variables is zero and the collinearity is reduced. We also scaled all the variables, dividing them by their standard deviations to interpret parameter estimates from models at a comparable scale. Because HFP is a continuous covariate, there are an infinite number of values we can use to analyse the effect of biodiversity on multifunctionality. For a better interpretation of the interactive effect, we selected three values (thresholds) of the scaled HFP: a mean value (0), a value of standard deviation above the mean (1) and a value of standard deviation below the mean (−1). This is a common approach to analyse interaction between continuous predictors 50 . These three HFP values can be interpreted as three levels of HFP intensity, low intensity (below average), moderate intensity (on average) and high intensity (above average). The slopes of each relationship between HFP and species richness, functional diversity and ecosystem multifunctionality are similar among wetlands, suggesting absence of any bias in our results ( Supplementary Fig. 10).
Pathways by which human pressure affects multifunctionality. To disentangle the direct and biodiversity-mediated pathways by which HFP affects multifunctionality, we ran SEM using the R package 'lavaan' 51 . Considering that all seven organismal groups worked in combination to determine multifunctionality (Figs. 2 and 3), we used their diversity to construct composite variables in our SEM. We combined the species richness and functional diversity of the seven organismal groups to construct a composite index for species richness and functional diversity, respectively. A composite index collapses the effects of multiple related variables into a single composite effect, thus representing a good way to analyse complex multivariate relationships in SEM 52 . We accounted for six ecosystem drivers: distance from equator, climate (mean annual temperature and precipitation) and aquatic characteristics (pH, conductivity and water level) in the SEM. The SEM was fitted based on a meta-model ( Supplementary Fig. 11). We calculated the standardized direct coefficients for each pathway within the model. We also estimated the indirect effect of HFP on multifunctionality mediated by diversity (species richness and functional diversity) of single organismal groups. To do so, we multiplied the coefficient of HFP on diversity (species richness and functional diversity) of a given organism group by the standardized loading of this organism group on composite. Finally, we multiplied the above result by the coefficient of composite on multifunctionality (Supplementary Table 11). We applied multi-group analysis in the SEM to evaluate whether (1) the effects of selected predictors (HFP, biodiversity, climate, space and aquatic properties) on multifunctionality and (2) the effect of HFP on biodiversity varied across wetlands. We considered the four wetlands as the grouping variable (Amazon, Araguaia, Pantanal and Paraná). We constructed a SEM in which all parameters were free to differ between wetlands and a model in which all parameters were fixed (that is, constrained to a single value determined by all wetlands). We compared the free model with the constrained model, where non-significant differences indicated no variation in pathway coefficients by wetlands, whereas significant difference indicated that pathway coefficients varied by wetlands. Because we found significant differences between the free and restricted/constrained model for both species richness and functional diversity, our next step was to understand which pathways differed. We analysed only the differences (multi-group) of the pathways including multifunctionality and biodiversity (species richness and functional diversity; Supplementary Table 10). Differences between other pathways within the model were not analysed. We evaluated the SEM fit using the comparative fit index (CFI; the model has a good fit when CFI ≥ 0.95) and the root mean square error of approximation (RMSEA; the model has a good fit when RMSEA ≤ 0.05). For our species richness model, the CFI was 0.997 and the RMSEA was 0.041, and for our functional diversity model, the CFI was 0.998 and the RMSEA was 0.026, indicating a good model fit. All analyses were conducted in R version 3.4.4 (ref. 53 ).
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
The data that support the findings of this study are publicly available on Zenodo Digital Repository at https://doi.org/10.5281/zenodo.6406782. Source data are provided with this paper.

Code availability
The code that supports the findings and figures of this study is available on Zenodo Digital Repository at https://doi.org/10.5281/zenodo.6406786. Last updated by author(s): May 31, 2022 Reporting Summary Nature Portfolio wishes to improve the reproducibility of the work that we publish. This form provides structure for consistency and transparency in reporting. For further information on Nature Portfolio policies, see our Editorial Policies and the Editorial Policy Checklist.

Statistics
For all statistical analyses, confirm that the following items are present in the figure legend, table legend, main text, or Methods section.

n/a Confirmed
The exact sample size (n) for each experimental group/condition, given as a discrete number and unit of measurement A statement on whether measurements were taken from distinct samples or whether the same sample was measured repeatedly The statistical test(s) used AND whether they are one-or two-sided Only common tests should be described solely by name; describe more complex techniques in the Methods section.
A description of all covariates tested A description of any assumptions or corrections, such as tests of normality and adjustment for multiple comparisons A full description of the statistical parameters including central tendency (e.g. means) or other basic estimates (e.g. regression coefficient) AND variation (e.g. standard deviation) or associated estimates of uncertainty (e.g. confidence intervals) For null hypothesis testing, the test statistic (e.g. F, t, r) with confidence intervals, effect sizes, degrees of freedom and P value noted

Data analysis
All the analyses were conducted using the free R software and language.
For manuscripts utilizing custom algorithms or software that are central to the research but not yet described in published literature, software must be made available to editors and reviewers. We strongly encourage code deposition in a community repository (e.g. GitHub). See the Nature Portfolio guidelines for submitting code & software for further information.

Data
Policy information about availability of data All manuscripts must include a data availability statement. This statement should provide the following information, where applicable: -Accession codes, unique identifiers, or web links for publicly available datasets -A description of any restrictions on data availability -For clinical datasets or third party data, please ensure that the statement adheres to our policy The data that support the findings of this study are available at https://doi.org/10.5281/zenodo.5764910. The R code used to calculate the all statiscial analysis are publically archived at zenodo: https://doi.org/10.5281/zenodo.5765674.
Policy information about studies involving human research participants and Sex and Gender in Research.

Reporting on sex and gender
No human were used in this study

Population characteristics
No human were used in this study

Recruitment
No human were used in this study

Ethics oversight
No human were used in this study Note that full information on the approval of the study protocol must also be provided in the manuscript.

Field-specific reporting
Please select the one below that is the best fit for your research. If you are not sure, read the appropriate sections before making your selection. All studies must disclose on these points even when the disclosure is negative.

Study description
Here, we used a large-scale field study to provide a novel test of the linkages between aquatic biodiversity (including taxonomic and functional diversity of seven organismal groups) and multifunctionality in Neotropical wetlands. We also tested how human pressure affect the relationship between aquatic biodiversity and multifunctionality.

Research sample
We used a total of 72 lakes distributed across four Neotropical wetlands in this study. Our data included sampling for seven aquatic organisms (Fish, macrophytes, microcrustaceans, rotifers, ciliates, testate amoebae, and phytoplankton) and 11 variables representing ecosystem functions Sampling strategy Our sampling was explicitly designed to assess aquatic biodiversity and ecosystem functions at the small shallow lake level, which are the most common freshwater system worldwide. At each lake, standardized samples were realized during a dry period and a rainy period. In these samplings were collected seven organismal groups: fish, macrophyte, microcrustaceans, rotifers, ciliate, testate amoebae, and phytoplankton, and 11 variables presenting ecosystem functions. Each organismal group was collected following common sampling methods for these groups. A total of 144 samples were analyzed in this study. We calculated lake-level estimates of aquatic biodiversity and ecosystem functions as explained in the Material and Methods of our paper.

Data collection
A total of 144 samples were realized in this study. Our field study included information on ~1,465 aquatic species, whereby 325 fish, 87 macrophyte, 99 microcrustacean, 124 rotifers, 598 phytoplankton, 124 testate amoebae, and 108 ciliate. Moreover, data included 11 ecosystem variables regulated by aquatic organisms and belonging to a wide range of ecosystem functions and properties: nutrient concentrations, ecosystem metabolism, multitrophic standing biomass, availability of photosynthetically active radiation, microorganism abundance, variation of underwater habitat complexity.
Timing and spatial scale Sampling of the wetland lakes occurred on a scale of 3,700,000 km² between the years 2011 and 2012

Data exclusions
No data were excluded from the analyses.

Reproducibility
Detailed description of the methods are provided, and all code and data necessary to reproduce the findings will be freely available on Zenodo.

Randomization
Surveys were conducted around 72 Neotropical wetland lakes: 16 in Amazon, 18 in Araguaia, 18 in Pantanal, and 20 in Paraná. Lakes of similar size were chosen haphazardly before the study. Initial underwater surveys were started at a random point adjacent to each lake. Replicate surveys (2011-2012) were conducted around the same lakes, and surveys were started at same points in 2011 and 2012.
Blinding N/A Did the study involve field work?
Yes No