Human pressure drives biodiversity– multifunctionality relationships in neotropical wetlands

Dieison Moi (  dieisonandrebv@outlook.com ) UEM https://orcid.org/0000-0002-7946-9260 fernando Lansac-Tôha Universidade estadual de Maringá Gustavo Romero University of Campinas Thadeu Sobral-Souza Federal University of Mato Grosso Bradley Cardinale Penn State University Pavel Kratina https://orcid.org/0000-0002-9144-7937 Daniel Perkins University of Roehampton https://orcid.org/0000-0003-0866-4816 Franco Teixeira de Mello Universidad de la República Erik Jeppesen Aarhus University Jani Heino Finnish Environment Institute https://orcid.org/0000-0003-1235-6613 Fábio Lansac-Tôha Universidade Estadual de Maringá Luiz Velho Universidade Estadual de Maringá https://orcid.org/0000-0001-8111-4955 Roger Mormul Programa de Pós-graduação em Ecologia de Ambientes Aquáticos Continentais, Núcleo de Pesquisas em Limnologia, Ictiologia e Aquicultura, Universidade Estadual de Ma https://orcid.org/0000-0001-90204784


Introduction
Human activities are causing biodiversity to decline worldwide 1,2 , which has led to interest in how biodiversity loss might change the functioning of ecosystems 3 . To date, most studies have revealed strong positive effects of biodiversity on the ecological processes that underlie ecosystem functioning, particularly the multiple ecological processes behind the so-called 'multifunctionality' of ecosystems [4][5][6][7] .
The strong relationship between biodiversity and ecosystem multifunctionality that has been documented experimentally and/or in semi-natural ecosystems has led to the prediction that the relationship between biodiversity and multifunctionality may be weakened in human-dominated ecosystems 3,8 . However, this prediction has yet to be adequately tested 9 .
Here we used a unique dataset from 72 lakes distributed across four 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 wetland multifunctionality. These four large wetlands provide a unique opportunity to test the in uence of human pressures across broad spatial scales as the lakes span a 3,700,000 km² gradient of distinct human pressures (Fig. 1). We quanti ed human pressure on the wetland ecosystems using the Human Footprint (HFP) index (see Methods). The HFP is a recently developed index that incorporates eight human pressures into a standardized cumulative index of human pressure 10 . Thi index provides an interesting opportunity to understand how human pressures are affecting biodiversity-multifunctionality relationships in natural to human-dominated systems.
HFP was estimated independently for each lake on the global map of HFP 10 , as was the species richness and functional diversity of seven organismal groups: sh, aquatic macrophytes, microcrustaceans, rotifers, phytoplankton, ciliates, and testate amoebae. These organismal groups comprised 1,465 plant, animal, and microbial species. We further quanti ed ecosystem multifunctionality in each lake 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 under water (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.
Because no studies have examined the large-scale relationships between biodiversity and ecosystem multifunctionality on wetlands, we rst established whether species richness and the functional diversity of seven organismal groups were, in fact, related to multifunctionality as previous small-scale evidences suggests 11,12 . After con rming the consistent relationship, we then determined if biodiversitymultifunctionality relationships changed as a function of the human pressure. Lastly, we used structural equation models (SEMs) to investigated the direct and indirect biodiversity-mediated pathways by which human pressure in uenced multifunctionality in wetlands. As we show, the normally positive relationship between biodiversity and ecosystem multifunctionality changes both quantitatively and qualitatively when humans put increasing pressure on biodiversity.

Results And Discussion
Biodiversity of seven organismal groups was strongly and positively associated with ecosystem multifunctionality of wetlands (Supplementary Table 1). The biodiversity-multifunctionality relationships were consistent among the four wetland and occurred for both species richness and functional diversity ( Fig. 2). This nding persisted regardless of how the measures of multifunctionality were weighted ( Supplementary Fig. 1). The multithreshold approache provided additional evidence showing that the mean minimum threshold at which species richness of the seven organismal groups had its strongest effects on multifunctionality averaged 57% (range 5-92%, Supplementary Fig. 2). Similarly, the mean minimum threshold at which functional diversity had its strongest effects on multifunctionality was 91% (range 70-99%, Supplementary Fig. 2). Finally, biodiversity of organismal groups were highly and positively correlated with most individual ecosystem variables related to standing biomass, nutrient concentrations, ecosystem metabolism, habitat complexity and light availability ( Supplementary Fig. 3).
Our large-scale dataset revealed strong associations between the biodiversity of multiple organismal groups and ecosystem multifunctionality. Notably, the effect of biodiversity over multifunctionality is high in wetlands. These results support the hypothesis that biodiversity of multiple aquatic organismal groups has a signi cant effect on the ability of the wetlands to provide more functions working at high performance [11][12] . Consequently, the strong biodiversity-multifunctionality relationships suggest that losses of biodiversity might reduce the ability of wetlands to maintain their functioning [4][5][6][7] . To test such prediction, we examine how the strength of biodiversity-multifunctionality relationships changed across a gradient in human pressure. As HFP increased, the strength of biodiversity-multifunctionality relationships not only decreased in magnitude but changed fundamentally in direction. The slope of the relationship between species richness and multifunctionality shifted from being positive to negative with increasing HFP intensity (Fig. 3c, Supplementary Table 2; LMM estimate 0.83 species richness to -0.22 species richness x HFP). Likewise, the slope of the relationship between functional diversity and multifunctionality was positive at low HFP and negative at high HFP (Fig. 3d, Supplementary Table 2; LMM estimate 0.77 only functional diversity to -0.15 functional diversity x HFP).
The results displayed in Figure 3 demonstrate that the relationship between biodiversity and ecosystem multifunctionality breaks down when human pressures increase. To determine why this breakdown occurs, we disentangled the direct and biodiversity-mediated, indirect pathways by which human pressures affect multifunctionality. The direct effect of HFP on multifunctionality was negative in all wetlands (Fig. 4, Supplementary Table 4). This nding is consistent with fact that wetlands usually cover regions of moderate-high HFP intensity (Fig. 1). For example, the study wetlands cover simultaneously crops of soy and sugarcane, pasturelands grazed by cattle, and high human population density [13][14][15] .
Consequently, cumulative of multiple human pressures cause direct impacts on wetland multifunctionality, driven by decline in standing biomass and habitat complexity, and changes in nutrient cycling and metabolism of ecosystem 12,16 .
In addition to direct effect, HFP had a large negative effect on species richness and functional diversity, which ultimately resulted in negative effects on multifunctionality via biodiversity-mediated pathways (Fig. 4). Importantly, the indirect negative effects of human pressure also were consistent across the four wetlands (Fig. 4b). Because biodiversity had the strongest direct effect on multifunctionality relative to all other environmental and climate covariates ( Fig. 4; SEM: direct effect of species richness on multifunctionality 0. 329, direct effect of functional diversity on multifunctionality 0. 415) our ndings suggest that the loss of biodiversity with increasing human pressure will leads to a net loss of multifunctionality. Combined with the fact that positive biodiversity-multifunctionality relationship is breaks down with increasing HFP (Fig. 3), these results highlight that, if the human pressures continue to increase 17 , preservation of biodiversity for maintaining ecosystem multifunctionality will not be su cient unless they are accompanied by reduction in human pressures. In the light of the increasing human in uence on natural landscapes, our results illustrate the importance of considering multiple pathways through which human pressures can in uence ecosystem functioning.

Conclusion
We have provided the rst empirical evidence for a positive large-scale relationship between the diversity of multiple organismal groups and the multifunctionality of wetland ecosystems. The consistent positive association of species richness and functional diversity with multifunctionality highlight the importance of conserving biodiversity to maintain multifunctionality and their associated services [2][3][4][5][6][7]18 , and indicates that biodiversity conservation should be a key management priority in wetlands. However, 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 results illustrate that human pressures are causing multifunctionality to decline through multiple pathways, indicating biodiversity conservation alone will likely not be enough to sustain multifunctionality if underlying human pressure is not reduced. This will be a major challenge for the management of ecosystems worldwide, as human pressures are increasing 19,20 . More broadly, reducing human pressures must be addressed urgently, particularly in wetlands 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² 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 21 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 year-1 22 . The eld data were collected between August and May 2011 and 2012. The wetland lakes were surveyed under the Brazilian program "National System for Research in Biodiversity" (Sisbiota Brazil). The eld 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 21 . In order 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 sh, aquatic macrophytes, microcrustaceans (cladocerans and copepods), rotifers, phytoplankton, testate amoebae, and ciliates. A detailed sampling protocol for each taxon is available in the Supplementary Methods. Diversity measure. We quanti ed the species richness of the focal organismal groups in all 72 lakes. After identifying each individual to species level, we determined 325 sh 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, but the locations differed in total number of individuals present. Therefore, we calculated estimated species richness as the Chao index with abundance-based data using the R package iNEXT 23 , which is based on rarefaction and extrapolation of Hill numbers and provides and 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 4,7 . 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 tness, and which also in uence ecosystem processes 24 . These traits fall into the three broad categories: (i) body size (maximum body length for animal taxa or cell volume for phytoplankton), (ii) resource and habitat use traits (feeding groups for animal taxa, growth form for macrophytes, nitrogen xation or mixotrophy for microorganisms), and (iii) mobility traits (migration ability for animal taxa, propagation method for macrophytes, and cell motility for microorganisms). The literature sources used for functional classi cation and the predicted impact of each trait on ecosystem functions can be found in Supplementary Table 1. In order to determine the functional diversity (FD) of each organismal group, we calculated functional dispersion -i.e., the mean distance in multidimensional trait space of the individual species to the centroid of all species 25 . This measure provides a robust estimate of functional diversity.
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 (see Supplementary Table 2). These functions and properties included: (i) 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 quanti ed according to Mackereth et al. 26 , while phosphorus was quanti ed following 27 . (ii) Ecosystem metabolism represented by the daily variation of dissolved oxygen in the water (mg L -1 day -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 28 . (iii) Multitrophic standing biomass was represented by the biomass of algae, carnivorous sh, omnivorous sh, herbivorous sh, and detritivorous sh. Algae standing biomass was quanti ed using biovolume (individuals per mm L -1 ) of identi ed algae species. Biovolume was estimated by multiplying the abundance of each species by their mean volume. The algae volume was obtained from geometric models similar to three-dimensional shapes 29 . Fish were classi ed into trophic groups using feeding trials and gut contents 19 . Afterwards, the sh counts within each trophic group were converted to biomass (g m -2 ) using published species-speci c length-weight relationships 30 . (iv) Availability of photosynthetically active radiation represented by light availability under water (m). We quanti ed light availability under water by the depth of the euphotic zone, which represents the depth (m) of the lake where there is su cient 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 under water 22 . (v) Microorganism abundance (cells mL -1 ) was quanti ed 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 asks. Bacteria were analyzed from water samples treated with a xative solution composed of alkaline Lugol's solution, borate buffered formalin, and sodium thiosulfate that was ltered through black Nuclepore lters (0.2 and 0.8 μm, respectively) and stained with uorochrome DAPI (4,6-diamidino-2-phenylindole 31 . Bacterial quanti cation was done with an epi uorescence microscope at a magni cation of ×1000 (Olympus BX51). (vi) Variation of underwater habitat complexity was quanti ed based on variations in the above-ground cover of macrophytes (m -²). We estimated the area of all leaves and culm of each macrophyte species. We then summed the area of all leaves and culm to obtain the above-ground area cover by each individual. We calculated the standard deviation of the above-ground area cover between all macrophyte species and used this standard deviation as a proxy of variation in the above-ground vegetal cover.
Trade-offs between functions. To assess the potential for a trade-off between individual ecosystem characteristics, we calculated pearson correlation coe cients between each pair of individual standardized functions. Of the possible 45 combinations of pairs of functions, we found only seven strong correlations (r > 5; Supplementary Fig. 5). To remove any bias in our multifunctionality index, the highly correlated functions were down-weighted it its calculation ( Supplementary Fig. 6), as described in Manning et al. 32 . We then repeated the analyses also for the weighted multifunctionality to determine whether the results differed when using weighted vs non-weighted multifunctionality (see ref. 32 ).
Assessing ecosystem multifunctionality. To obtain robust and quantitative multifunctionality indexes for each lake, we used three independent multifunctionality approaches: (1) the averaging multifunctionality index, (2) the multithreshold multifunctionality index, and (3) the multiple single functions index 32 . To obtain the averaging ecosystem multifunctionality index, we standardized all 10 individual ecosystem variables 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 4,33 , but has the limitations 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 and to determine whether the effect of biodiversity on multifunctionality changed across wetlands systems, we used the multithreshold index. This index calculates how many functions simultaneously exceed a prede ned percentage of the maximum observed value of each individual function. Because the selection of any threshold is arbitrary, analyzing multiple thresholds of maximum functioning is recommended 33 . Thus, we analyzed 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. The number of functions surpassing the thresholds was calculated with the R-package "multifunc" 33 .
Assessing the Human-Footprint on wetlands. We used the global Human Footprint (HFP) map as a surrogate of the cumulative human-induced pressure on the wetlands 10 . This map is constructed from eight variables of human activities: (i) the extent of built environments, (ii) crop land, (iii) pasture land, (iv) human population density, (v) night-time lights, (vi) railways, (vii) roads, and (vii) navigable waterways. To facilitate comparison across pressures, each human pressure is placed at a scale varying from 1 to 10 (not all pressures range across the full 0-10 scale and details on the weightings are provided below). The pressures were weighted according to estimates of their relative intensity 10 . For example, (i) constructed environments are areas related to urban settlements such a buildings, paved land, and urban areas. The pressure of built environments was assigned a score of 10 (i.e., a score of 10 is assigned if there are built environments in the wetland area, otherwise a score of 0 is assigned). (ii) Crop land is characterized by monocultures with high inputs of pesticides and fertilizers on wetland lakes. 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. (iii) Pasture land included 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 per cent pasture for each 1 km 2 pixel. (iv) Human population growth 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² were assigned a pressure score of 10. For less populated areas, we logarithmically scaled the pressure score using the following estimation: Pressure score = 3.333 x log (population density + 1). (v) Night-time lights include electric infrastructure related to more rural and suburban areas that are not part of built environments. To estimate the pressure of night-time lights, the areas were divided into 10 quantiles of increased night-time light intensity associated with scores between 1 and 10, while areas with no lights were assigned a zero score. (vi) Railways are essential human infrastructures that in uence 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. (vii) 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). (viii) Similar to roads and railways, navigable waterways act as conduits for people to access nature, resulting in impacts on natural ecosystems. The pressure of navigable waterways was assigned a score of 4, which decayed exponentially out to 15 km away from the water banks. Full details of HFP estimation see refs 2,10 .
The HFP map is available as 1 km² cell-size resolution using ArcGis 10.1 in a global extension showing the HFP in 2009 (most recent global estimation). For any grid cell, the HFP can range between 0 and 50 (sum of all human pressures). These pressures are not mutually exclusive, and many co-occur in the same wetland may also vary strongly between lakes within the same wetland. Even though the HFP was initially developed to represent human pressures in terrestrial systems 10 , all these human activities also extensively affect aquatic ecosystems. For instance, Brazil has experienced a rapid urban expansion 34 , with human settlements being constructed near wetlands 16 . Along with the increase in human populations in the vicinity of wetlands, there has been an increased pressure on these ecosystems from cattle and sheep pastures, built environments, railways, roads, and navigable waterways 35 .

Statistical analyses
The relationship between aquatic diversity and multifunctionality. First, to determine the direct link between aquatic biodiversity and average multifunctionality across all wetland ecosystems, we tted a series of linear mixed effects models (LMMs) to the survey data. Speci cally, we tested the relationship between (i) species richness and (ii) functional diversity of seven organismal groups and ecosystem multifunctionality. We assumed a Gaussian error distribution using the function lme in the 'nlme' package 36 . We nested the two sampling periods within each wetland as our random structure to account for random variations between seasons and lakes within each wetland. The assumptions of normality, linearity, and homoscedasticity were veri ed using graphical diagnostics (QQ plots and residual plots). To improve normality, some species richness values were log-transformed prior to the analysis. We also performed Spearman correlations between species richness and the functional diversity of each group of organisms and each ecosystem variable. We did this to test the relationship between biodiversity and single ecosystem functions.
We also modeled aquatic biodiversity against the number of functions above a threshold using generalized linear mixed effects model (GLMMs), assuming a Gaussian error distribution in the MASS package 37 . Because we wanted to know whether the relationships between species richness and functional diversity with ecosystem multifunctionality changed as a function of organismal group and among the four wetlands, we tted the GLMM individually to each organismal group. We then extracted and plotted the linear coe cient ( tted 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 tted curve for each wetland at multiple thresholds.
Effect of human pressure on relationship between biodiversity and multifunctionality. To determine whether Human Footprint (HFP) altered the relationship between biodiversity and ecosystem multifunctionality, we added interaction terms for HFP × species richness and HFP × functional diversity to the mixed-effects models and measured the estimate coe cients for these interactions. Considering that species richness and functional diversity of all organismal groups had a positive relationship with multifunctionality (Fig. 2), we decided to use a synthetic index re ecting the full species richness and functional biodiversity of the wetland lakes 38 . To estimate the overall species richness, the species richness values of each organismal group were standardized by scaling them to the maximum observed value, and then the averages of these standardized richness values were calculated 38 . We used this mean value as a diversity index because it encompasses the richness of all organismal groups simultaneously. We did not sum species richness values to calculate diversity because this would give a higher weighting to species-rich groups. Likewise, the overall functional diversity was calculated as the mean standardized FD values across all organismal groups. The overall diversity index (species richness and functional diversity) ranges between 0 and 1. This index has been widely used in studies of the relationships between biodiversity and ecosystem functioning 7,38-39. Pathways by which human pressure affects multifunctionality. Finally, to disentangle the direct and biodiversity-mediated pathways by which Human Footprint (HFP) affects multifunctionality, we conducted a structural equation model using the R package piecewiseSEM 40 . SEM was constructed using species richness and functional diversity, accounting also for multiple predictors such as human pressure, climate (mean annual temperature and precipitation), and environmental covariates (pH, conductivity, and water level; Supplementary Fig. 7). Temperature, pH, conductivity, and water level were all measured in situ using an analytical thermometer, portable digital potentiometers (Digimed), and a water level ruler, respectively. Mean precipitation data were extracted from the WorldClim database (http://www.worldclim.org) at a resolution of 1 km² (30 arcsec). To evaluate whether the strength of the human pressure effects on multifunctionality varied across wetlands, we constructed the SEM using multi-level analysis considering wetlands as groups (Amazon, Araguaia, Pantanal, and Paraná). Speci cally, we implemented a model-wide interaction and tested each path interaction within the model. A signi cant path interaction indicates that the path (relationship) differs among wetlands. The aim of this analysis was to identify which paths are consistent (or differ) across wetlands when including the key drivers of multifunctionality. We tested multicollinearity for each endogenous predictor by calculating the variance in ation factor (VIF), where VIF > 3 indicates possible collinearity. We found that all the predictors had VIF < 3, indicating that there were no signi cant correlations. To reduce the number of predictors and to build more parsimonious models, we performed a model selection using Akaike's information criteria corrected for small sample size (AICc). In this selection, the full model (including all predictors) is compared with the reduced models (with some predictors removed) using AICc (AICfullmodel -AICreducedmodel). We considered traditional ΔAICc > 2 units to distinguish the complete from the reduced models and thus permanently remove the predictor. Notably, the full and reduced nal model differed in at least ΔAICc = 73 units (Supplementary Table 7). The lack of direct or indirect effects on multifunctionality was used as a criterion to remove these predictors. We used this criterion to remove pH, conductivity, and temperature from the full model. We calculated the standardized coe cients for each path within the model. Also, we estimated the indirect effect of human pressure on multifunctionality via multiplication of human pressure coe cients on biodiversity by the coe cient of biodiversity on multifunctionality. In this approach, only signi cant paths are considered, i.e. whether human pressure has signi cant effects on biodiversity and whether biodiversity has signi cant effects on multifunctionality, and we therefore assume that human pressure has an indirect signi cant effectmediated by biodiversity -on multifunctionality. Path signi cance was obtained by maximum likelihood and model t using Shipley's test of d-separation through Fisher's C statistic (p > 0.05 indicates no missing path and the model is tted). All analyses were conducted in R version 3.4.4 41 .

Data availability
The data that support the ndings of this study are publicly available on Zenodo Digital Repository at https://doi.org/10.5281/zenodo.5764910

Code availability
The code that supports the ndings and gures of this study is available on Zenodo Digital Repository at https://doi.org/10.5281/zenodo.5765674 <p><strong>Intensity of the human footprint (HFP) across Brazil and the four major neotropical wetlands (Amazon, Araguaia, Pantanal, and Paraná)</strong>. Activity data maps of the wetlands (built environments, crop land, pasture land, human population density, night-time lights, railways, roads, and navigable waterways) used in the human footprint analysis in this study were extracted from (refs<sup>10</sup>). The HFP data varied from 50 to 0 according to the pressure of multiple human activities. The HFP pressures in the four wetlands differed. The HFP data ranged from 0 to 50 according to the pressure from multiple human activities. However, the HFP data on the four focal wetlands <p><strong>The relationship between the biodiversity of aquatic organisms and multifunctionality in neotropical wetlands. </strong>The linear link between multifunctionality and the diversity of the seven selected organismal groups (including species richness and functional diversity); n = 72 lakes. Solid black and dashed colored lines are extracted from linear mixed-effect models for overall and local trends (to each wetland ecosystem), respectively. Trends from LMM show that the links between species richness and functional diversity with average multifunctionality are mostly positive for all groups of aquatic organisms and wetlands. *P &lt; 0.05, **P &lt; 0.01, ***P &lt; 0.01. R²<sub>c</sub> = conditional (variance of the xed plus mixed effects) and R²<sub>m</sub> = marginal (i.e., variance of the xed effects). &nbsp;Full model results are provided in Supplementary Table 3.</p><p><br></p> Figure 3 <p><strong>Effect of Human Footprint (HFP) on the relationship between biodiversity and ecosystem multifunctionality in neotropical wetlands. a,b, </strong>Links of multifunctionality with richness <strong>(a)</strong> and functional diversity <strong>(b)</strong> in wetland lakes. Importantly, biodiversity attributes were created based on a synthetic index re ecting the total ecosystem species richness and functional diversity (i.e., across seven organismal groups; Allan et al. 2014). Solid black and dashed colored lines are extracted from linear mixed-effect models to display overall and local trends (to each wetland ecosystem), respectively, and the shaded areas represent 95% CIs of the overall trend. <strong>c,d,</strong> Estimate coe cients for the interaction term between HFP and species richness <strong>(c)</strong> and functional diversity <strong>(d)</strong> from the same models for multifunctionality. Points represent estimates, thick lines represent 75% CIs, and thin lines represent 95% CIs. Complete model output is provide in Supplementary Table 4.</p><p><br></p> Figure 4 <p><strong>Link between human footprint and biodiversity-multifunctionality relationships</strong>. <strong>a</strong>, Using a SEM, we aimed to analyze the direct relationship between the human footprint and the combined species richness and functional diversity of seven groups of aquatic organisms and their links with multifunctionality. We also accounted for environmental drivers (temperature, precipitation, pH, conductivity, and water level) in the SEM model. Information on model selection and simpli cation steps using Akaike Information Criteria (AIC) are provided in Supplementary  Table 7 and full model results are provided in Supplementary Table 6. The nal model tted the data well (<em>P </em>= 0.646, Fisher's C = 0.875). Solid black arrows represent signi cant paths to overall model (<em>P</em> ≤ 0.05 piecewiseSEM). Arrows for non-signi cant paths (<em>P</em> ≥ 0.05) are in dashed gray. The colored coe cients indicate multigroup analysis effects for each individual wetland. The thickness of the signi cant paths (arrows) represents the magnitude of the standardized regression coe cient. Numbers adjacent to arrows are indicative of the standardized effect size of the relationship. R<sup>2</sup>s for component models are given in the boxes of variables. Signi cance levels of each predictor are *<em>P</em> &lt; 0.05, **<em>P </em>&lt; 0.01, ***<em>P</em> &lt; 0.001. <strong>b</strong>, Standardized indirect effects of the human footprint on multifunctionality mediated by species richness and functional diversity. Effects are derived from the structural equation model and displayed relative to overall effects (black bar) and to each individual wetland (multigroup analysis; colored bars).</p>

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download. MoietalSupplementaryMatterial.docx