Introduction

The study of biodiversity variation between geographic regions is a key issue for disciplines like ecology, zoology, conservation, and management (^{Tuomisto, 2010}; ^{Anderson et al., 2010}). Here, we are especially interested in variation of benthic macroinvertebrate communities associated with different spatial scales. These animals possess desirable attributes for studying such variation. They are highly diverse and distributed in most freshwater stream habitats under different environmental characteristics (^{Hussain & Pandit, 2012}), reflecting the streams’ health. The taxon is essential for energy flux, nutrient cycling and organic matter degradation in lentic and lotic environments, being a central component of freshwater ecosystems (^{Whiles & Wallace, 1997}; ^{Allan & Castillo, 2007}). Organisms are easily collected due to their high abundance, and the community patterns can rapidly change in response to biotic and/or to abiotic characteristics, like habitat heterogeneity and water quality, and at local, regional or geographic scales (^{Brosse, Arbuckle, & Townsend, 2003}; ^{Heino, Muotka, & Paavola, 2003}; ^{Hepp & Santos, 2009}; ^{Rezende, Santos, Henke-Oliveira, & Gonçalves Jr, 2014}).

At local scale, stream’s sediment composition, allochthonous material, and canopy openness are the main variables controlling diversity of benthic macroinvertebrate communities (^{Harper, Mekotova, Hulme, White, & Hall, 1997}). Natural topographic characteristics can also affect invertebrate diversity if fine sediments are predominant in the stream (^{Vannote, Minshall, Cummins, Sedell, & Cushing, 1980}). Physical, chemical and structural features are also important for aquatic life; for instance, dissolved oxygen, electric conductivity, alkalinity, and temperature are known to influence invertebrate occurrence and distribution (^{Melo, 2009}; ^{Rezende et al., 2014}). Additionally, the communities’ composition changes and the richness increase from the headwaters to the base level in response mainly to stream’s width, depth, temperature and production/respiration ratio (Vannote et al., 1980; ^{Jacobsen, 2004}; ^{Baptista et al., 2014}).

Widespread distribution of benthic macroinvertebrates is possible because they have evolved a diverse spectrum of morphological and behavioral mechanisms for exploiting foods, which are used to classify their functional role in aquatic ecosystems (^{Cummins & Klug, 1979}; ^{Merritt & Cummins, 2006}; ^{Hamada, Nessimian, & Querino, 2014}). Feeding groups or guilds of benthic macroinvertebrates are defined as sets of taxa exploring the same class of food resources in a similar manner (^{Merrit, Cummins, & Berg, 2008}). Differences in the distribution of feeding groups in streams of a region may indicate the integrity of aquatic system (^{Fu et al., 2015}; ^{Castro, Dolédec, & Callisto, 2017}). Therefore, studying benthic macroinvertebrate feeding group distributions in aquatic ecosystems is a useful tool for management and conservation (^{Callisto & Esteves, 1998}; ^{Callisto, Moreno, & Barbosa, 2000}; ^{Merritt, Fenoglio, & Cummins, 2016}; ^{Echelpoel et al., 2018}).

Environmental assessment of aquatic invertebrate responses to anthropogenic disturbance helps to understand which factors drive community structure changes (^{Hepp et al., 2016}; ^{Castro, Dolédec, & Callisto, 2018}; ^{Ongaratto, Loureiro, Restello, & Hepp, 2018}). Aquatic ecosystems in subtropical Brazil are presently under intense environmental pressure due to riparian vegetation suppression and allochthonous inputs from watercourses’ surroundings (^{Ferreira, Cyrino, Duarte-Neto, & Martinelli, 2012}). Inserted in this climatic region, South Braziliangrasslands encompass regions of Pampa biome extending from the southernmost part of Brazil to Uruguai and Argentina, and grasslands at highlands in Atlantic Forest biome of South Brazil (^{Overbeck et al., 2009}). The grasslands guarantee key environmental services like hydrological resource conservation and biodiversity (^{Pillar, Müller, Castilhos, & Jacques, 2009}). In the last decades, the conversion of grasslands to agriculture, silviculture, and pasture have reduced the natural grassland vegetation area by 50 % (^{Cordeiro & Hasenack, 2009}). These anthropogenic activities impact the aquatic ecosystems, leading to multiple physical-chemical and structural changes in the water courses, affecting taxa occurrence and distribution (^{Utz, Hilderbrand, & Boward, 2009}). Land use effects on freshwater communities are complex and scale dependent (^{Allan, 2004}). At landscape scale of streams’ catchments, human activities like agriculture can increase sediment and nutrient load, as well as flood frequency, and decrease groundwater recharge (^{Carpenter, Stanley, & Zanden, 2011}). At local scale, replacement of native vegetation by human activities may, for instance, change input of solar radiation, thus affecting primary productivity (^{Burrell et al., 2013}).

Based on the responses of aquatic communities to local and landscape environmental features, our objective was to compare the benthic macroinvertebrate communities in streams of three South Brazilian grassland landscapes. To do so, we measured communities’ richness, abundance and feeding group compositions, and assessed the role of variables associated with local (water chemistry and riverbed structure) and landscape (catchment land use) scales, and with geographical distance as drivers of macroinvertebrate community variations. We hypothesized that streams with more preserved habitat conditions linked to less disturbance in surrounding land have higher taxonomic richness and more diverse feeding groups.

Materials and methods

Study area: We performed the study in three South Brazilian grassland (Campos)landscapes (Fig. 1). The region of South Brazilian grasslands is characterized by mosaics of shrubland and forest patches embedded in a matrix of natural grassland. Historically, the grasslands were used for livestock, mainly sheep and cattle for beef production (^{Overbeck et al., 2007}). Recently, a large part of these natural grasslands was converted to croplands, mainly wheat, soybean, rice and maize. The landscapes were chosen based on the presence of at least 50 % of natural grassland cover after satellite imagery inspection in Google Earth and local validation. Landscapes are inserted in two large river basins: (i) Paraná river - city of Palmas (UTM 22 J 437612.67 m E - 7065595.88 m S; 1 115 m.a.s.l.) and city of Tibagi (UTM 22 J 575509.83 m E - 7282637.15 m S; 748 m.a.s.l.); and (ii) Uruguai river - city of Painel (UTM 22 J 580844.12 m E - 6901796.10 m S; 1 144 m.a.s.l.). In Palmas (PAL), the collections were carried out at the Campos de Palmas Wildlife Refuge, a full protection unit located in the South-central region of Paraná state. In Tibagi (TIB), the collections were carried out in the State Park of Guartelá full protection unit in the central-eastern portion of the state of Paraná. Painel (PAI) landscape is not part of a protection unit. Surroundings of streams in PAI are mainly composed of areas of natural grasslands typically used for cattle breeding and fewer areas of plantations and silviculture (Table 1). Silviculture area in the streams’ catchments was larger in PAL landscape, and TIB showed larger area of plantations (Table 1). The stream sediments in PAI and PAL are composed mostly by boulders (> 256 mm), gravels (4 - 64 mm), and a mixture of sands (0.06 - 1.2 mm) (^{Wentworth, 1922}). Streams also have dry leaves, branches and twigs characterizing a heterogeneous habitat. On the other hand, sediments in TIB streams are dominated by sands, characterizing a less heterogeneous habitat. All streams had riparian forest on both margins.

Fig. 1 Location of the south Brazilian region under study in South America. 1 - River basin of Tibagi/PR showing the grid of the 25 km² landscape and the location of the five selected streams. 2 - River basin of Palmas/PR showing the grid of the 25 km² landscape and the location of the five selected streams. 3 - River basin of Painel/SC showing the grid of the 25 km² landscape and the location of the five selected streams. 

Invertebrate sampling: We surveyed 15 streams (2^nd to 4^th order), five in each landscape (Table 1). Each stream was sampled once. Inside each stream we chose an easy access point to start the sampling and go through 150 m downstream. Along that distance, we collected benthic fauna in 20 backwater points. At each point, the researcher kicked 10 times the sediment inside a sweeping net, method called kick and sweep net (K&S), which is routinely used for assessing biodiversity (^{Tubic et al., 2017}). The sediments were placed inside plastic bags and milk pails with absolute alcohol. Streams were sampled between June and November 2015.

In the laboratory, the sediments were washed in 1.0-, 2.0- and 3.0-mm-mesh sieves. The invertebrates from the 3.0-mm-mesh sieve were sorted by visual inspection. The material of the other two sieves were fixed in 70 % alcohol. The great amount of sandy sediments in each sample made extensive sorting of macroinvertebrates unviable. Therefore, we used the collector’s curve asymptote criterion to set the maximum number of samples of sandy sediments from each river. We started arbitrarily with a sample of 7.5 g of sediments and counted the number of morphotypes. We repeated the process until the collector’s curve reached the asymptote for each river (^{Magurran, 2013}). After 10 samples (75 g) all streams’ curves reached the asymptote and we stopped the macroinvertebrate sorting. Invertebrates were identified to the lowest taxonomic level according to ^{Hamada, Nessimian, and Querino (2014}), ^{Dominguez and Fernandez (2009}) and ^{Souza, Costa, and Oldrini (2007}). The taxa were classified in functional feeding groups of (i) shredders/detritivores - which consume leaf litter or other coarse particulate organic matter, (ii) predators - which feed on other consumers, (iii) scrapers - which consume algae and associated material, (iv) collectors/filters - which collect fine particulate organic matter from the water column using a variety of filters, and (v) gatherers/collectors - which collect fine particulate organic matter from the stream bottom (^{Merrit, Cummins, & Berg, 2008}; ^{Hamada et al., 2014}).

Environmental variables: We measured eight environmental variables inside and in the surroundings of the catchments (Table 1). Locally, variables were mean water velocity (three measurements in the beginning, middle and final points of the 150 m), mean stream’s width (same three points), pH, dissolved oxygen, and water temperature. We set a 500 m radius buffer initiating in the starting point of macroinvertebrate sampling and, inside each buffer, we measured the relative cover of forests, grasslands, plantations and silviculture (in ha) in Google Earth.

Data analyses: We performed a principal component analysis (PCA) with all environmental variables to characterize the three landscapes. We transformed all variables to zero mean and unit variance before analysis. We used the first and second principal components - PC1 and PC2 of PCA as predictor variables of macroinvertebrate community responses. Richness (number of taxa) and abundance (number of individuals) of macroinvertebrates were highly correlated (R= 0.812). It means that any analysis of possible correlations between environmental variables and macroinvertebrate richness would be biased by abundance. Therefore, we performed a simple linear regression with abundance (X) and richness (Y) and used the vector of residuals as the new richness of macroinvertebrates. To assess whether richness and abundance differed among landscapes, we used two permutational analyses of variance (^{Wheeler & Torchiano, 2016}). The test criterion of permutational ANOVAs is the sum of squares (SS). To identify which predictor variables affected the response of richness and abundance, we performed two multiple linear regression with sequential entering of predictor variables, first PC1 and then PC2. Abundance of macroinvertebrates was log transformed before multiple regression to meet normality assumptions.

We explored the similarity of benthic macroinvertebrate feeding group compositions among streams with a non-metric dimensional scaling ordination (NMDS) based on Bray-Curtis dissimilarity of feeding groups abundance data. We performed a permutational multivariate analysis of variance (PERMANOVA) after NMDS to check for differences among landscapes based on the feeding group compositions. Additionally, we tested the effect of space and environment on feeding group compositions by means of partial Mantel correlation tests using the dissimilarity/distance matrices. First, we decomposed the variation explained by space and environment into community composition variation among landscapes. Second, we analyzed within landscape composition variations. Analyses were performed in RStudio v. 3.5.1 (^{R Core Team, 2016}).

Results

A total of 1 058 benthic macroinvertebrates were collected divided in 53 taxa (Digital Appendix). Bivalvia sp. (N = 363) was the most abundant taxon in all streams, followed by Aegla sp. (N = 175) and Lymnaeidae sp. (N = 147). In PAL Bivalvia sp. was the most abundant taxon (N = 306), whereas in PAI it was Lymnaeidae (N = 147) and Aegla sp. (N = 133). Finally, the streams of TIB showed greater abundance of Limnocoris sp. (N = 21). Abundance was significantly different among landscapes (SS = 25660; P = 0.038; d.f. = 2.12), however the residual richness was not different among landscapes (SS = 678.9; P = 0.888; d.f. = 2.12). TIB streams showed nearly one third lower mean abundance and around 10 % lower mean richness than streams of PAL and PAI (Table 1). Abundance responded significantly to PCA 1 predictor (R² = 0.568, P = 0.006) (Fig. 2). Residual richness did not respond to environmental variables (R²= 0.168, P = 0.330).

Fig. 2 Relationship between abundance (log transformed) of benthic macroinvertebrates and PCA axis 1 of environmental variables. Data collected in 15 streams of three landscapes in South Brazilian grasslands in Winter and Spring 2015. PAL = Palmas, PAI = Painel, TIB = Tibagi. Response model is Abundance = 3.72 - 0.60*PCA1. 

The first two axes of PCA explained 51.61 % of total variation in the matrix of environmental variables. Table 1 shows stream scores and variable loadings in the first and second PCA axes. The first axis varied from TIB streams to PAI and PAL streams. TIB was characterized by higher temperature and area of forest in the catchment and lower pH and mean width. PAI and PAL were characterized by lower temperature and area of forest and higher pH and width. Second axis varied from two PAL streams with higher plantation area in the catchment to all the rest of PAL and PAI streams associated with higher grassland area and water velocity.

Regarding the compositions of functional feeding groups, the main variation was associated with first axis of ordination (Fig. 3), which divided the streams of TIB from the ones in PAI and PAL. Results of PERMANOVA showed that the three landscapes were different regarding macroinvertebrate feeding group compositions (F_2,14 = 7.45; P = 0.001). TIB streams had more predators like Hemiptera and Odonata. PAI and PAL streams showed more diverse communities of shredders/detritivores, gatherers/collectors, collectors/filterers, and scrapers (Fig. 3). Additionally, along the second axis of NMDS, PAI streams were better described by the abundance of scrapers, compared to PAL streams (Fig. 3). PAL streams showed higher abundance of collectors/filterers Bivalvia, Hydropsychidae and Simuliidae, and the gatherers/collectors Elmidae. PAI streams also had higher diversity compared to TIB and showed the presence of the shredder/detritivore Aeglidae and the gatherers/collectors Ceratopogonidae.

Fig. 3 First and second NMDS axes based on Bray-Curtis dissimilarity distances between feeding group abundances in 15 streams of three landscapes in South Brazilian grasslands. PAI = Painel, PAL = Palmas, TIB = Tibagi. Feeding groups are Co.Fi = collectors/filterers, Ga.Co = gatherers/collectors, Sc = scrapers, Sh.De = shredders/detritivores, Pr = predators. Stress of NMDS = 0.059. 

Among landscapes, partial Mantel tests showed that the highest correlation was between feeding guild composition and space (R_M = 0.67; P = 0.001), compared to environment (R_M = 0.38; P = 0.002). Additionally, space and environment were significantly correlated (R_M = 0.32; P = 0.001). It means that environmental variation was linked to geographical distance, which is expected, given the extent region of the study. Within each landscape, the dissimilarity in macroinvertebrate compositions in TIB were correlated with the environmental distance (R_M = 0.79; P = 0.025).

Discussion

The lack of communities’ local responses to environmental variables measured at higher scale does not immediately translate into a lack of streams’ surroundings effect since high unexplained variation is common in many ecological studies, partly due to missing variables in the analyses (^{Genner et al., 2004}; ^{Hepp, Landeiro, & Melo, 2012}). On the other hand, although the landscapes are predominantly grassy habitats, we could observe marked differences mainly in the taxonomic and functional feeding groups of streams in the northernmost TIB landscape.

Landscapes studied here encompass a wide geographical scale. Latitudinal gradients in aquatic invertebrates are ambiguous, so patterns of variation are study-dependent. Therefore, local environment features can supposedly override historical and climatic influences on diversity (^{Hillebrand, 2004}; ^{Heino, 2009}). Hence, the type of substrate could partially explain the lower diversity found in TIB streams. Texture, compaction, and grain size can affect the composition and abundance of benthic macroinvertebrates (^{Nakamura & Kikuchi, 1996}; ^{Fidelis, Nessimian, & Hamada, 2008}; ^{Baptista et al., 2014}). Sandy soils are characterized by low heterogeneity and availability of physical habitats (^{Allan, 2004}). Additionally, these habitats are subjected to disturbances due to the stream’s current. Despite the greater coverage of riparian vegetation in TIB, allochthonous materials are drifted and not deposited in the stream bed due to the absence of larger substrates. Contrasting features were found in PAL and PAI, leading to a more heterogeneous environment. The occurrence of pebbles and coarse fractions of substrate increases the number of shelters (^{Bücker, Sondermann, Frede, & Breuer, 2010}; ^{Jun et al., 2011}; ^{Hepp et al., 2012}; ^{Rezende et al., 2014}), translating into more taxonomic and functional diversity. TIB water temperature was higher than in PAI and PAL, which may be related to the lower altitude of this landscape and the sampling period at the beginning of the warmer months of the year.

TIB benthic macroinvertebrate diversity was different from those of the other landscapes. The presence of predators was represented particularly by Gomphidae (Odonata) and Naucoridae (Hemiptera). The presence of Gomphidae may be associated with its fossorial attribute, which enables the organisms to use inorganic substrates, such as sand and gravel (^{Assis, Carvalho, & Nessimian, 2004}; ^{Worthen, Gregory, Felten, & Hutton, 2004}; ^{Worthen & Horaceh, 2015}). Hemiptera are common predators in lentic and lotic habitats, with good swimming capacity and adapted to various habitats, including sand sites, which may explain their presence in TIB (^{McCafferty, 1981}; ^{Hamada et al., 2014}; ^{Reynoso-Velasco & Sites, 2019}). In addition to these taxa, we recorded the presence of the family Odontoceridae, known for the construction of tubular shelters made of grains of sand glued with silk, attribute that facilitates the occupation of different niches compared to other species of caddisflies. Moreover, this taxon could search more actively for food (^{Guevara, Reinoso, & Villa, 2005}; ^{Vásquez-Ramos, Guevara-Cardona, & Reinoso-Florez, 2014}) in the sandy habitats found in TIB streams.

Despite the high abundance of scrapers such as Lymnaidae in PAI streams, we also found representants of the EPT group (Ephemeroptera, Plecoptera, and Trichoptera), which are the first to disappear in highly disturbed streams due to their sensitivity to pollution and thus generally considered bioindicators of good water conditions (^{Barbola et al., 2011}; ^{Júnior, Conceição, Lobo, Santos, & Sardinha, 2019}). Another sensitive taxon that appeared in all PAI and PAL streams was Aegla sp. Aeglids live preferentially in well-oxygenated and clean water habitats (^{Bond-Buckup & Buckup, 1994}; ^{Dalosto & Santos, 2011}; ^{Santos et al., 2017}). On the other hand, the presence of Oligochaetes in PAL landscape could be signaling the beginning of water deterioration, because these organisms can survive under low dissolved oxygen conditions (^{Piedras et al., 2006}, ^{Rosa, Rodrigues, Oliveira, & Alves, 2014}). Moreover, changes in water quality associated with organic contamination lead to the increase of suspended particles, favoring filter-collector organisms (Bivalvia) (^{Feio & Dolédec, 2012}).

Different benthic macroinvertebrate functional groups explore different resources and are influenced by habitat features, food, and refuge availability (^{Callisto, Moreno, & Barbosa, 2000}; ^{Benstead &Pringle, 2004}). Therefore, the more variation in these environmental features the greater the diversity and the resistance of benthic macroinvertebrate community to disturbance (^{Hawkins & MacMahon, 1989}). Thus, there could be two different but not excluding forces driving the overall lower diversity of TIB streams. First, the low heterogeneity is driven by the presence of one dominant sandy substrate. Second, the anthropogenic land use of agriculture and extensive cattle breeding in the region that is ongoing and started before the creation of the conservation unit.

In Palmas landscape, the land use alterations also started before conservation unit creation, extending to today. So, these past environmental alterations could have led to riparian forest and landscape structural changes, factors that affect the benthic macroinvertebrate richness and composition (^{Hepp & Santos, 2009}; ^{Sensolo, Hepp, Decian, & Restello, 2012}; ^{Hepp et al., 2016}). The replacement of specialist organisms by generalists occurs as a response to changes in land use (^{Petsch, 2016}; ^{Castro et al., 2018}). Additionally, the spatial pattern of distribution of the populations depends on factors reflecting the spatial organization of habitat patches (^{Vinatier, Tixier, Duyck, & Lescourret, 2011}). Therefore, high variation in the local context of streams could lead to high dissimilarity of taxonomic and functional compositions within the landscapes.

There was a marked structural difference among the three landscapes at local and landscape scales. Studies conducted with different organisms have shown that both local environmental conditions and geographical distances shape community structure (^{Heino et al., 2010}; ^{Santos, Silva, Branco, & Huszar, 2015}; ^{Monteiro do Amaral, De Almeida Gonçalves, Da Silveira, & Gama Alves, 2019}). At regional level the observed variations in aquatic biota were partially associated with landscape changes, but the main driving force was geographical distance. Campos in PAL and TIB landscapes were more fragmented compared to PAI. The PAI landscape, despite not being included in a protected area, had greater coverage of natural grasslands, which can be related to the traditional grazing activity (e.g. cattle, horses) that was almost exclusive in PAI landscape. Given this condition, our initial hypothesis that diversity would be higher in more conserved areas was partially corroborated.

Intensive agriculture can change soil physical and chemical properties, compaction, erosion and water regime (^{Vieira & Overbeck, 2015}). However, other environmental variables should be analyzed (e.g. heavy metals, phosphorus, nitrogen) to better understand the diversity patterns found in the three Campos landscapes, mainly in regions of agriculture and livestock farming and silviculture practices, where the water courses are subjected to pesticide and fertilizer residues from the land. Besides, benthic invertebrate communities can respond to other unmeasured variables either at local or at landscape scales. We believe it is important to deepen the studies on natural grassland regions, a historically neglected biome in Brazil (^{Overbeck et al., 2009}), to assess the conservation status of landscapes and to understand the distribution of aquatic diversity at local, landscape and geographical scales. The results point to the need for discussions about public policies addressing water conservation and management, which should include a mixed collegiate body in Advisory Councils (e.g. community, researchers, authorities). These actions may focus on discussing the role of the grassland landscapes studied here. Although the streams are in protected areas, they still suffer from anthropic activities, which can impact the local and regional benthic macroinvertebrate diversity.

Ethical statement: authors declare that they all agree with this publication and made significant contributions; that there is no conflict of interest of any kind; and that we followed all pertinent ethical and legal procedures and requirements. All financial sources are fully and clearly stated in the acknowledgements section. A signed document has been filed in the journal archives.