We carried out this study in ten headwater streams (1st to 3rd order, sensu Strahler 1957) belonging to the Pirapó River basin (Upper Paraná River system). The Pirapó River presents a drainage area of 5,000 km², with approximately 168 km of extension to its mouth, in the Paranapanema River (Pagotto et al. 2012). Streams were sampled within the rural and urban areas of the municipality of Maringá, Paraná State, Brazil (Fig. 1). A mixture of temporary croplands (mainly soy, corn, and sugar cane) and urban areas predominate in the landscape, which stands out among the three most populous cities in the Paraná State. This region integrates the transition between tropical and subtropical climate and is classified as a permanently hot humid rain climate zone, Cfa (h), according to the Köppen scale (Alvares et al. 2013). Mean annual temperature varies between 16° and 20°C, with January being the hottest and wettest month and July the month with coldest and driest records. In general, the annual rainfall rate in the region exceeds 1,000 mm (Alvares et al. 2013).
Sampling was carried out in April and May 2017. We selected 30 sampling sites, encompassing three mesohabitats (riffles, runs, and pools) from five rural (Atlântico, Lombo, Queçaba, Romeira, and Roseira) and five urban (Miosótis, Maringá, Mandacarú, Morangueira, and Guaiapó) streams (Fig. 1). Before the collection day, we visited several stretches along each stream to visually select the mesohabitats according to the description by Rincón (1999), with the following characterization: (i) riffles presenting fast and turbulent waters, with a substrate composed mainly of large rocks; (ii) runs with relatively fast waters, but deeper and less turbulent than riffles; (iii) pools presenting deep and slow water, with fine sediment the most common substrate. The prior visual selection of the mesohabitats ensured the standardization of their quantity, with each stream presenting a riffle, run, and pool. Also, each mesohabitat had to be at least ten meters long, since this length was used to standardize the size of the sample units, where the following variables were recorded: current velocity (with a JDC electronic flowmeter, model Flowatch FL-K2), channel depth and width, proportions of flooded vegetation, canopy shading and substrate type (sand, rock, clay, and artificial substrate), dissolved oxygen (O2; DIGIMED, model DM-4P), pH (DIGIMED, model DM-22) and electrical conductivity (DIGIMED, model DM-32). As there was no strictly aquatic vegetation in any stream, we considered the roots, trunks, and branches of the riparian forest as flooded vegetation. Also, construction waste and trash were treated as artificial substrates.
The stream width was measured in three transects, comprising upstream, downstream, and in the intermediate portion along each mesohabitat. The other variables were measured on the right and left margins and in the middle of each of the three transects, comprising nine collection points for each mesohabitat. We quantified the proportions of flooded vegetation, canopy shading by riparian vegetation, and substrate type using a 0.25 m2 wooden square subdivided into 25 squares smaller than 0.01 m2 and estimated their values from the sum of the filled squares. After quantifying these variables, we calculated their averages to characterize the mesohabitats according to their environmental conditions.
We collected fish using electrofishing with a 2,500W alternating current generator operated at 500V and 2A through successive passages deployed for 30 minutes in each mesohabitat. To optimize sampling effort and reduce fish escape, we blocked each mesohabitat downstream using a 5mm mesh size. We anesthetized the captured fish with benzocaine as an immersive solution for at least 10 minutes or until the opercular movement ceased. Subsequently, we placed the specimens in vials containing 4% diluted formaldehyde, to be transferred to 70% alcohol 72 hours later. We counted and identified the individuals according to Graça and Pavanelli (2007) and Ota et al. (2018) to be checked according to the list of fish fauna from the Paraná State (Reis et al. 2020) and the Pirapó River basin (Pagotto et al. 2012). Testimony specimens were deposited in the Ichthyological Collection of the Research Center in Limnology, Ichthyology, and Aquaculture of the State University of Maringá (voucher numbers: NUP 20040 to NUP 20128). The survey license was granted by the Instituto Chico Mendes de Conservação da Biodiversidade (ICMBIO; nº. 25560-1).
Two morphological measures related to the body shape of the fish species were obtained, namely body height (BH) and body width (BW). BH and BW were taken from ten adult individuals of each species using a digital caliper (0.01 mm approximation), and for species that did not have such abundance, all individuals were measured. From the average of these two morphological measures for each species, the Compression Index (CI) was calculated by dividing the BH with BW. Thus, species with the highest CI values have a laterally compressed body, while the lowest values are for species that have a dorsoventrally depressed body (Winemiller 1991).
Fish diversity was assessed by species richness and the Shannon-Wiener index. The statistical differences in these assemblage attributes were assessed with a two-way ANOVA, including the interaction of the land use types and mesohabitats in the model. For these ANOVAs, we removed the pool site of the Lombo stream, as it was an outlier in these analyzes, presenting one and zero values for species richness and Shannon-Wiener index, respectively (Table S1 in Supplementary Information). Because each stream is considered three times in the analysis (three mesohabitats), stream identity was used as a blocking factor (additive factor) in ANOVAs and in all the models made below, to control its effect on the model variances
To assess the differences in species composition among mesohabitats, an abundance sites × species matrix was subjected to Hellinger standardization (Legendre and Gallagher 2001), from which the Bray-Curtis distance was calculated. From this distance matrix, a Principal Coordinate Analysis (PCoA) was used to visually assess the mesohabitat ordination regarding the differences in relative abundances of species. Subsequently, the 12 environmental variables collected were standardized for zero mean and unit variance and correlated with the PCoA scores using the envfit function of the vegan package. This function considers the environmental variables as the dependent variables that are explained by the ordination scores, and each dependent variable is analysed separately. Only the significant variables (p<0.05 based on 999 permutations) were added to the ordering graph.
The statistical differences in PCoA ordination according to land use types and mesohabitats were evaluated by Permutational multivariate ANOVA (PERMANOVA based on 999 permutations). Pairwise tests (based on 999 permutations) were then performed to assess the differences between all combinations of land use and mesohabitat factors. The assumption of homogeneity of variances was evaluated and met by permutation tests of multivariate dispersion (PERMDISP based on 999 permutations).
The indicator species of each mesohabitat were evaluated by the procedure adopted by Dufrêne and Legendre (1997) considering the probability <5% as indicator values (IndVal, based on 999 permutations) using the labdsv package (Roberts 2019). For this purpose, we regarded the combination of the levels between land use and mesohabitat factors on a matrix of species abundance. For both PERMANOVA and IndVal, rare species (with no more than 1–2 occurrences) were removed because they do not influence IndVal p-values and disproportionately influence the separation of points in PERMANOVA. Thus, seven species were removed from these analyzes, which were carried out with 19 species (Table 1).
To assess the difference in the body shape of the fish species among mesohabitats, a matrix was generated with the average values of the Compression Index (CI) weighted by the species abundance for each sampled site. For this, the CI values of each species were multiplied by the abundance matrix of the 19 species used in PERMANOVA and IndVal using the SYNCSA package (Debastiani and Pillar 2012). Then, an ANOVA was performed on the mean IC values using the interaction between the factors land use and mesohabitat as predictor variables. The assumptions of normality and homogeneity of variance were evaluated and met for all ANOVAs by the Shapiro-Wilk and Levene’s tests, respectively.
To verify whether the composition of fish species was related to the proximity among sampling sites, we performed a Mantel correlogram with the hydrological distance among all sites using the mpmcorrelogram package (Matesanz et al. 2011). The species composition matrix was generated through the Bray-Curtis distance calculated on the Hellinger standardization of the species abundance. The hydrological distance was obtained from the calculation of the distance, in kilometers, between the collection points in the QGIS program (QGIS Development Team 2018), using the stream network of the Pirapó River basin, downloaded from the website of the Instituto Água e terra do Paraná (IAT 2021), and the geographical coordinates of the collection points. We performed all analyses in the R program (R Core Team 2020). The PCA and MANOVA were performed using the stats package, while the vegan package (Oksanen et al. 2020) was used to running the ANOVAs, PERMANOVA, PERMDISP, and envfit function. All graphics were generated by the ggplot2 package (Wickham 2016).