Local and regional determinants of phytoplankton communities in water reservoirs from the Cerrado biome

Abstract Aim Based on a study comprising water reservoirs within the Brazilian Cerrado biome, we seek to answer the following question: how does phytoplankton communities respond to environmental, landscape, and spatial predictors? We expected local predictors to be the main factors structuring the communities. Since phytoplankton has high dispersal capacity, the geographical distance would be of minor importance. Methods: We collected phytoplankton samples from 40 water reservoirs in the rainy season and 37 reservoirs in the dry season. We performed a partial Redundancy Analysis (pRDA) to evaluate the factors influencing the variation in the composition of phytoplankton communities. Results We found that spatially structured environmental factors were controlling phytoplankton communities at the rainy season, whereas landscape was the main predictor in the dry season. On the other hand, phytoplankton morphofunctional groups were influenced only by local predictors. Conclusions We demonstrated that phytoplankton dynamics differs between rainy and dry seasons, and that distinct predictors affect phytoplankton communities over seasons.

The metacommunity theory provides an approach that supports the investigation of the roles that dispersal processes and environmental heterogeneity play in structuring the dynamics of communities (Leibold et al., 2004). In this way, communities can be studied from a regional perspective, in which local (environmental heterogeneity) or regional (dispersal) factors can alter the structure of phytoplankton communities (Bie et al., 2012). Thus, the morphological characteristics of individuals may influence these dynamics mainly due to their relationship with dispersal capacity (Bortolini et al., 2019).
Among several factors that influence the dispersal capacity of organisms, the connectivity and the distance between water bodies are fundamental to the success of this process (Chisholm et al., 2011;Naselli-Flores et al., 2016). The management of water bodies must consider their position and size in the landscape to increase habitat connectivity (Izaguirre et al., 2016;Lamy et al., 2013;Lansac-Tôha et al., 2016) and elaborate strategies for aquatic biodiversity conservation (Guan et al., 2019;Hill et al., 2019;Zhai et al., 2019).
In this study, we evaluated the influence of local variables (physical and chemical), as well as spatial predictors (distance between water bodies), and landscape predictors (land cover and use) on the composition of phytoplankton taxonomic and morphofunctional groups (MFG) in reservoirs located at the Cerrado biome. We tested the following hypotheses: (i) both local and spatial characteristics are essential factors structuring the phytoplankton communities for both taxonomic and MFG approaches, with the predominance of local predictors; (ii) the spatial predictors is more important during the dry season than during the rainy season. The lack of active dispersal capacity of phytoplanktonic organisms and the need for vectors for their transport (water current, animal hair/feathers, and wind) between the patches of habitats cause limitations in the movement of these organisms, especially in the dry season, characterized by reduced connectivity between the reservoirs, which can increase the influence of the spatial predictor (Bortolini et al., 2017;Bie et al., 2012;Heino et al., 2015).

Introduction
Among several aquatic systems, reservoirs play an essential role in conservation and biodiversity maintenance (Céréghino et al., 2014), as these ecosystems compose about 30% of the surface of water bodies (Downing et al., 2006) and contribute considerably to the biodiversity of landscapes (Williams et al., 2004). Also, reservoirs provide environmental services that assist water resource management, such as sediment retention and reduced nutrient loads in the watershed (Céréghino et al., 2014;Yasarer et al., 2018).
However, the functioning and maintenance of these ecosystems depend on complex relations between several components, such as resources, abiotic conditions, community composition, and trophic interactions. Among these components, biological diversity is fundamental because it affects ecosystem processes and stability (Thébault & Loreau, 2005). In addition, species that present redundancy in their ecological functions increase the stability of ecosystems (Walker, 1992). Thus, it is important to use functional approaches in studies that aim to understand the functioning of ecosystems (Salmaso et al., 2015).
One of the main primary producers of aquatic ecosystems is the phytoplankton community (Reynolds, 2006). This photoautotrophic group is hugely diverse and has a wide range of functional characteristics [e.g., sizes, shapes, toxin production, and productivity rate; (Kruk et al. (2010)]. The dynamics of this community are structured by several factors, such as local environmental conditions (e.g., light, nutrients, temperature) and biological interactions (e.g., competition and predation) (Machado et al., 2016;Wu et al., 2018). Furthermore, recent studies have found that spatial factors related to dispersal processes are crucial for determining the composition of this community (Bortolini et al., 2017;Oliveira et al., 2020). The spatial structure of communities is affected by the high dispersal capacity of phytoplanktonic organisms, which occurs mainly passively (Naselli-Flores et al., 2016). Moreover, regional factors greatly affect the spatial structure of aquatic communities, such as connectivity and land use (Devercelli et al., 2016;Izaguirre et al., 2016). Phytoplanktonic organisms developed several adaptations to overcome the barriers that terrestrial environments present, such as mucilage production, increased cell walls, spores, akinetes, and cyst formation (Ribeiro et al., 2011). area of 10,200 km 2 and is part of the São Francisco River basin. The climate in the region is classified as Aw based on the Köppen classification, where rainfall ranges from 1,600 to 1,900 mm. year -1 and the annual mean temperature ranges between 19 to 20 ºC. It has two well-defined seasons, one of drought, which occurs from April to September, and another of rain, which runs from October to March (Alvares et al., 2013). The land use in this basin consists of highly mechanized agriculture using medium to large areas of central pivots (Borges et al., 2007). For this reason, most of the water is destined for irrigation, which comes from reservoirs (Carneiro et al., 2007). We collected samples from 40 reservoirs in the rainy season (February 2018) and from 37 reservoirs in the dry season (September 2018). We sampled the same reservoirs; however, three reservoirs were dry during the dry season and could not be sampled.

Phytoplankton community
We collected 100 mL of water at the subsurface (ca. 50 cm) in the central point of each reservoir to quantify phytoplankton and we stored it in dark amber flasks. Subsequently, we fixed the phytoplankton with Lugol solution. We estimated phytoplankton density using a Zeiss inverted microscope with 400× magnification following Utermöhl (1958). We identified the organisms to the lowest possible taxonomic level and the density was expressed in individuals per milliliter (ind.mL -1 ) (Bicudo & Menezes, 2006;Komarek & Anagnostidis, 1983;Komarek & Fott, 1983). Then, we classified the species found according to their morphofunctional groups (MFG), based on the classifications by Salmaso & Padisák (2007).

Local predictors
We measured water temperature (ºC), turbidity (NTU), pH, dissolved oxygen (DO; mg.L -1 ), electrical conductivity (mS.cm -1 ), oxi-reduction potential (ORP; mV), and total dissolved solids (TDS; g.L -1 ) using the Horiba multiparameter probe (Model U-50) in locu. We collected 500 mL of water from the subsurface (ca. 50 cm). We stored it on ice until we performed the gas chromatography method (APHA, 2015) to determine cations and anions (nitrite, nitrate, phosphate, ammonia, magnesium, fluoride, chloride, bromide, sulfate, sodium, potassium, and calcium). Total phosphorus (TP) and total nitrogen (TN) were also determined by the ascorbic acid method (APHA, 2015) and by total organic carbon (TOC) by combustion method (APHA, 2015), respectively, but only during the dry season since we did not have access to the equipment at the date of the rainy season.

Landscape predictors
We used the Digital Elevation Model (DEM) from Advanced Land Observing Satellite (ALOS), Phased Array L-band Synthetic Radar (PALSAR), from the Alaska Satellite Facility (ASF) electronic page: https://asf.alaska.edu/. The data were Radiometric Terrain Corrected (RTC) high resolution, 12.5 meters spatial resolution DEM. This DEM generated catchment areas of each sample unit that were considered the exutory of each sub-basin (catchment areas), following the theoretical approach of analyzing the riverscapes (Allan, 2004;Morley & Karr, 2002). We used the Terrain Analysis Using Digital Elevation Models (TauDEM) available on hydrology.usu.edu/taudem/ taudem5 in ArcMap® to fill sinks, generate the flow direction, flow accumulation, stream definition (>0.001 km 2 ), and catchment area (sub-basin).
After generating the catchment areas, we obtained the land cover and use (forest formation, savanna, grassland formation, agriculture, and urban occupation). We downloaded the data from the Mapbiomas Project (https://mapbiomas.org), a multi-institutional initiative to generate the annual land cover and use maps using automatic classification processes applied to satellite images. We generated the area, in percent, for each sub-basin (catchment area) from the 2017 data available in Mapbiomas. We gathered the morphometric data from Google Earth images calculating the area and then the perimeter and margin development index (MDI) for each water reservoir (Håkanson, 1981).

Spatial predictors
We obtained the spatial variables by converting the latitude and longitude geographic coordinates to the Cartesian plane through the geoXY function of the SoDA package (Chambers, 2014) in the R environment (R Core Team, 2022). Subsequently, for independent ordering on orthogonal axes, we generated distance-based Moran Eigenvector Maps (dbMEM) . We performed this analysis using the dbmem function in the adespatial package (Dray et al., 2022). Then, we used the same forward selection steps described below.

Statistical analysis
To prevent rare species from significantly influencing the analyses, we standardized the species abundance matrix following the Hellinger method, using the standardize function in the vegan package (Oksanen et al., 2013). All statistical analyses were performed in the R environment (R Core Team, 2022).
We measured the collinearity between the variables within the local predictors' matrix for parsimony and reduction of the number of local predictors. We assessed this linear dependence with the Variance Inflation Factors (VIF), removing variables with VIF values greater than 10. After this, we selected the variables based on the knowledge about phytoplankton community ecology and using the forward selection procedure with two stopping criteria (Borcard et al., 2018). First, we performed a redundancy analysis (RDA), and if the test was significant, it was possible to perform the forward selection. Subsequently, to reduce the risk of adding too many variables, the p-value was used as the first stopping criterion, followed by the adjusted R 2 (Blanchet et al., 2008). Thus, we included only variables that presented p values ≤ 0.05 and adjusted R 2 lower than that of the global model. We performed this analysis using the rda and forward.selection functions in the adespatial package (Dray et al., 2022).
After screening for collinearity and based on the forward selection analysis and our knowledge about phytoplankton ecology, we selected seven local predictors for the taxonomic classification and eight for the morphofunctional groups in the rainy season. On the other hand, in the dry season, six local predictors were selected for the taxonomic classification and seven for the morphofunctional groups. Regarding spatial predictors, for the rainy season, we selected four Moran eigenvector maps for the taxonomic classification. However, for the dry season and MFG the spatial predictor was not significantly related to the phytoplankton community.
To analyze the amount of variance explained by (a) local variables, (b) spatial variables or (c) the set of all variables, we performed a variance partitioning using partial Redundancy Analysis (pRDA) for the taxonomic classification and the functional groups. We performed these analyses with the rda function in the vegan package.
Then, to complement information from the pRDA analysis, we carried out a Multiple Regression Tree (MRT) to categorize the relationships between species, MFGs, and environmental characteristics. Thus, groups of sample units were defined by limit values of explanatory variables (De'ath, 2002). Subsequently, we combined this model with the Indicator Species analysis (IndVal) to select the species that most contributed to the explained variance of each group (Borcard et al., 2018). We performed these analyses using the mvpart function in the mvpart package (De'ath, 2014), and we identified the discriminant species by the summary function of the package mvpartwrap (Ouellette & Legendre, 2012).
During the rainy season, we tested three predictors (local, spatial, and landscape) for the global model, and only the local (p = 0.001) and spatial (p = 0.001) predictors could significantly explain the variation in the phytoplankton community. For this reason, we performed the variance partitioning between local and spatial predictors, in which none of the partitions explained the phytoplankton community variance (Table 1).
The global models for local (p = 0.001) and landscape (p = 0.015) predictors significantly explained the phytoplankton community in the dry season. However, the community was not spatially structured (p = 0.329). Regarding to variance partitioning, only the landscape predictor explained the phytoplankton community variation ( Table 2).
The MRT analysis of the rainy season explained 15% of the phytoplankton variation (Figure 2). The model separated the reservoirs into four groups; the first separation factor was the margin development index (MDI), whose node explained 5% of the split. The species Peridinium volzii Lemmermann was related to values below 1.4. The second node divided the communities with MDI above 1.4 followed by turbidity below and above 5.2 NTU, explaining 5% of this split. Six taxa (Ankistrodesmus falcatus Ralfs, Hariotina reticulata P.A.Dangeard, Monoraphidium kormakovae Nygaard, Kirchneriella lunaris Möbius, Euglena polymorpha Dangeard, and Mucidosphaerium pulchellum C.Bock, Proschold & Krienitz) were related to turbidity values higher or equal than 5.2 NTU. Finally, the last node divided the reservoirs with turbidity below 5.2 NTU according to the dam area, below and above 0.6 ha. This node explained 5% of this split. Four taxa were related to values above or equal to 0.6 ha (Cryptomonas ovata Ehrenberg, Cryptomonas marssonii Skuja, Parvodinium umbonatum Carty, and Edaphochlamys debaryana Pröschold & Darienko), and three taxa were related to values below 0.6 ha (Euglena gracilis Klebs, Closterium setaceum Ralfs and Nitzschia sp.).
In the dry season, the MRT model explained 17% of the phytoplankton community variation (Figure 3), and four groups were formed. The first discontinuity divided the reservoirs into two groups  according to pH (below or above 5.4, explaining 5.7% of the split. The species Peridinium volzii was associated with pH below 5.4. The second factor that divided the reservoirs with pH above or equal to 5.4 was conductivity (below or above 0.07 mS.cm -1 ), explaining 4.9% of the split. Only the species Tetradesmus lagerheimii M.J. Wynne & Guiry was related with reservoirs with conductivity above or equal to 0.07 mS.cm -1 . Finally, pH explained 5.7% of the last split, separating the reservoirs with conductivity below 0.07 mS.cm -1 into groups with pH below or above 6.9. The species Pseudanabaena limnetica Komárek, Monoraphidium contortum Komárková-Legnerová, Mucidosphaerium pulchellum, Scenedesmus sp. 1, and Trachelomonas hispida F. Stein were associated with reservoirs with pH above or equal to 6.9.

Morphofunctional groups
We found 25 MFG groups in the rainy season. The 2b group was dominant, with 16% of all density, followed by 9a with 10.5%. The dry season presented 24 MFG groups, and the dominant groups were 5a and 9a with densities contributions of 23.1% and 18.8%, respectively. However, there was no difference in the average morphofunctional densities between the rainy and dry seasons (t-test = 1.27, p = 0.20).
The global models for the spatial and landscape predictors were not significant and only the local predictors significantly explained the phytoplankton community in the rainy (R 2 = 0.15 p = 0.001) and dry (R 2 = 0.065 p = 0.018) seasons. For this reason, we did not perform the variance partitioning.
The MRT models using the MFG explained more of the phytoplankton community variance than the taxonomic classification. In the rainy season (Figure 4), the area of the reservoirs explained 8% of the first split and separated the reservoirs into two groups, with sizes below and above or equal 0.3 ha. Only the MFG 5a was associated with reservoirs small than 0.3 ha. The second factor that separated the group with an area above or equal 0.3 ha was area again, but now it divided samples between reservoirs bellow and above 7 ha, explaining 8% of this split. The MFGs 5d and 5c were associated with the four biggest reservoirs. The group with area below 7 ha was then divided by turbidity values above and below 11.3 NTU. This variable explained 6% of this split. Two MFG groups (2c and 1c) were related to reservoirs with higher turbidity, whereas the MFG 2b was associated with reservoirs with lower turbidity.
The MRT model for the dry season explained 24% (Figure 5), and pH was the first factor that separated the data into two groups, below and above 6.8. This variable explained 9% of this split, and the MFGs 5d and 5c were associated with values above or equal to 6.8. Secondly, DO explained 8% and divided the reservoirs with pH below 6.8 into two groups with DO values below and above 8.5 mg.L -1 . The MFGs 2c and 1c were associated with higher DO values. Lastly, pH was again responsible for the split, now dividing the reservoirs with DO values above 8.5 mg.L -1 into pH values below and above 5.4. Now, this factor explained 7% of this split. Regarding the indicator MFG groups, the 2b group was related to higher pH values, and the 5a group was related to reservoirs with more acidic waters.

Discussion
We found that local predictors were not important for the phytoplankton community in any season, but in the rainy season we had indications that environmental and spatial relationships might be structuring the phytoplankton community. For the dry season, we found a greater influence of landscape on phytoplankton. Regarding  functional groups, only local predictors structured them. Finally, we also detected changes in the environmental variables that structure the reservoirs according to the seasons.
In the rainy season, we found no evidence that the spatial or local predictors individually influenced the phytoplankton community. However, the fraction that represents the shared variation of local and spatial predictors could not be tested. This may indicate that spatially structured environmental factors were controlling the phytoplankton community. The rainy season was related to nitrate, phosphate, temperature, pH, magnesium, and ammonia. Since agricultural landscapes predominate our study area, we suggest that those nutrients that are usually found in agricultural fertilizers would have been carried out to reservoirs by surface runoff (Silva et al., 2011). Furthermore, the water temperature can also increase due to runoff draining from impermeable surfaces of agricultural areas (Doubek et al., 2015). In this way, the dynamic of environmental factors with agricultural landscapes could explain their relationship with the rainy season in our study area.
Even though the MDI, turbidity, and reservoir area were not the most important variables affecting communities' structure, they were important for grouping the reservoirs during the rainy season. The MDI and reservoir area tend to be higher in the rainy season than in the dry season, which can produce more heterogeneous microhabitats and larger species pools (Smith et al., 2005;Tundisi & Matsumura, 2011). Furthermore, in the rainy season, high turbidity values are usually related to the sediment inflow from the surface runoff (Girardi et al., 2016;Reynolds, 2006), and some colonial chlorophyte species found in our study, such as Mucidosphaerium pulchellum, and Hariotina reticulata, are tolerant to low light environments (Reynolds et al., 2002), which explains the ability of these species to indicate high turbidity reservoirs.
In the dry season, the phytoplankton community was structured primarily by landscape predictors, mainly agricultural landscapes. Those landscapes are mainly associated with increases in nutrients and temperature in reservoirs (Silva et al., 2011). In addition, cyanobacteria were the most abundant class in the dry season. These organisms are normally found in agricultural landscapes and also in nitrogen-limited reservoirs (Doubek et al., 2015). This effect occurs mainly in the dry season when the N:P ratio is reduced (Hayes et al., 2015). The main variables selected for the dry season were conductivity, pH, nitrite, TN, and TP. Nutrients such as phosphorus and nitrogen had lower values in the dry season than in the rainy season, probably due to reduced elements of allochthonous origin (Silva et al., 2011). On the other hand, higher values of conductivity are common during the dry season because of the concentration of ions caused by the smaller amount of water in the reservoirs (Devercelli et al., 2016).
Conductivity and pH were crucial in generating groups in the reservoirs during the dry season; in general, the conductivity values were low in all reservoirs except those with predominant urban land use. This phenomenon was similar to pH, where urban reservoirs had pH close to 7, and the other reservoirs had values around 6. Four species of Chlorophyceae were identified as indicators of more neutral waters, and this class is generally associated with organisms without specialized traits and tolerance to pH variations (Chakraborty et al., 2011;Kruk et al., 2010).
Our first hypothesis that local predictors would be more important in explaining the variation in the phytoplankton structure but that spatial predictors would be more significant in the dry season was partially supported, as local predictors were only significant for MFG data. This greater predictive ability of local predictors over functional groups is expected because the functional attributes of organisms are directly linked to their ability to acquire resources, trophic interactions and their niche (Litchman et al., 2012;Litchman & Klausmeier, 2008). We can also observe that the environmental variables were similar between the functional groups and the taxonomic data at both seasons. Our results showed a complementarity of functional groups and taxonomic data, reinforcing that functional approaches are a fundamental strategy for understanding the responses of phytoplankton to environmental changes and also for biomonitoring purposes (Litchman & Klausmeier, 2008;Reynolds, 2006;Kruk et al., 2010).
At last, our results evidenced metacommunities dynamics where, during the rainy season, spatially structured environmental predictors could affect the phytoplankton community due to the proximity and greater connectivity of the reservoirs, resulting in a larger influence of the space in the local predictors (Heino et al., 2015;Rocha et al., 2020). On the other hand, landscape was more influential in the dry season, mainly agricultural land uses. This predictor is usually related to environmental conditions such as water temperature and nutrients and this association could be structuring the phytoplankton community in the dry season. In this sense, we showed that, for both seasons, the local predictors had some influence. However, for the rainy season, the spatial predictor was more important, while for the dry season it was the landscape predictor.

Conclusion
We demonstrated a complex relationship between the phytoplankton community and the local, spatial, and landscape predictors. While the local predictors had a small contribution to phytoplankton structure in both seasons, the spatial and landscape predictors were important in the rainy season and dry seasons, respectively. We also revealed that nutrients such as nitrogen and phosphorus were important in both seasons. The importance of other environmental characteristics changed between the seasons.
Finally, the MFGs are a crucial complementary tool to understand the dynamics of phytoplankton metacommunities in reservoirs. We suggest the use of these groups in studies that aim to unravel the complex relations of this community with its local predictors.