Acidification and increase of phosphorus levels in Pampean streams after 12 years of agricultural intensification

Agricultural intensification is a process that is still undergoing in many emergent economies. In the Pampas (Central Argentina), the former low-external-input farming was replaced by a model based on genetically modified crops and the intensive use of pesticides (mainly glyphosate) in the last decades. Here, we analyzed changes in water chemistry (pH, conductivity, dissolved oxygen, nutrients, and carbonates) in 41 streams of Buenos Aires province between 2003/2004 and 2015/2016, and the impact of geology, soil type, and land use change on water chemistry. We also used the SPAtially Referenced Regressions On Watershed attributes (SPARROW) model to analyze possible drivers of stream phosphorus loads. We observed modifications at reach scale in several streams, including changes in channel morphology, riparian vegetation, and land use, and a moderate expansion of agriculture in most basins. Mean nitrate concentration did not change significantly between 2003/2004 and 2015/2016. Dissolved phosphorus concentration increased in the streams, especially in the southern fluvial regions, but contrary to our expectations, phosphorus levels were not associated with land use but with pH. The SPARROW model also supported the link between water acidification and phosphorus concentration, and indicated that the whole basin acts as a phosphorus source. Possible explanations of fluvial acidification may be related to current agricultural practices, including higher inputs of labile carbon from croplands, soil acidification by nitrogenous fertilizers, and the generalized use of glyphosate. This highlights the necessity of adopting new agricultural paradigms to reduce the use of agrochemicals in Pampean basins.


Introduction
Fluvial ecosystems are strongly linked to the terrestrial systems where they are embedded, hence they are very sensitive to human activities in the surrounding landscape (Lowe and Likens 2005). Land use is considered the main driver of biodiversity change in streams due to increases in the input of nutrients, sediments, and contaminants (Sala et al. 2000). Among different land uses, agriculture probably has had the widest impact globally, producing declines in water quality and biological communities, and thus altering the ecological integrity of stream ecosystems (Allan 2004;Graziano et al. 2020;Young et al. 2021). Agricultural intensification developed in Europe and North America during the twentieth century, and resulted in widespread nutrient enrichment and eutrophication of freshwaters (Moss 2008;Weigelhofer et al. 2018). However, in some emerging countries such as Argentina, this process only began some decades ago. The Pampas (Central Argentina) is one of the most fertile regions of the world for crop production. Until the onset of the twentieth century, the region was covered by native grasslands, which were then replaced by an agricultural scheme that alternated cropping with extensive cattle breeding. This scheme was characterized by limited use of pesticides and fertilizers (Viglizzo et al. 2001). However, since the 1990s, a model of large-scale agriculture has been imposed in the region, based on the introduction of genetically modified crops (mainly soybean). This led to the adoption of new agricultural practices (introduction of no tillage, wheat-soybean rotation, etc.), changes in land use and cover (expansion of crop area and displacement of cattle breeding to marginal zones or feedlots), the reduction of the diversity of crops, and a strong increase in the use of fertilizers and pesticides (Viglizzo et al. 2001;Ferreyra 2006;Aizen et al. 2009;Manuel-Navarrete et al. 2009).
In Buenos Aires province, the core of the Pampean region, crop area increased by 53% between 2001/2002 and 2015/2016. In the same period, the percentage of crop area for genetically modified crops (soybean and corn) increased, while the area for traditional crops in the region (wheat and sunflower) was reduced (Table 1; MAGyP 2022).
Additionally, fertilizer consumption in Argentina increased from 1.8 to 3.2 million tons between 2001(CIAFA 2022. Given the magnitude of transformations of farming practices in the Pampean region, we expect that the impact of agriculture on fluvial systems may have increased in the last few years. In a previous study, we sampled 41 streams in Buenos Aires province between 2003 and 2004, when the process of agricultural intensification in the region was at an early stage. Sampling was made to establish a baseline water quality for Pampean fluvial systems; therefore, we selected streams that showed a relatively low level of anthropogenic disturbance. Comparison of these data with new information may allow to document the impact of agricultural intensification on the streams, and may provide further knowledge for other regions of the world that are undergoing the same process. This is especially relevant in the context of climate change that will further alter the ecological integrity of streams and rivers. In this study, we aimed to compare water chemistry (pH, conductivity, dissolved oxygen concentration, nitrogen and phosphorus levels, and carbonates) of 41 streams of Buenos Aires province between 2003/2004 and 2015/2016. Moreover, we analyzed the impact of soil type and land use change on water chemistry during this period. In the case of phosphorus, we also used a model to analyze the influence of different characteristics of the basin (land use, type of soil, etc.) on phosphorus load in the streams. We hypothesize that dissolved nitrogen and phosphorus levels are strongly linked to land use in the basins. Therefore, we predict that soluble reactive phosphorus (SRP) and nitrate will increase from 2003/2004 to 2015/2016, and that this increase will be associated with the expansion of agriculture. We also expect that the impact of land use change will be higher in regions of the province where the agricultural cover was lower in 2003/2004.

Study region
The study was conducted in the Pampean region (Buenos Aires province), a vast plain that covers central Argentina. The climate is temperate humid or subhumid. Mean annual precipitation ranges from 1100 mm in eastern Buenos Aires province to 700 mm in the west, and it is evenly distributed throughout the year, with slight increases at the end of summer (February) and at the onset of spring (October).
Soils are fertile and composed of a sequence of sand and loess deposited during the Tertiary and the Quaternary (Sala et al. 1983). The original vegetation is grassland and autochthonous trees are almost absent, except for two species (Celtis tala Gill. ex Planch and Salix humboldtiana Willd.) that are restricted to areas with particular soil conditions. The original grassy vegetation was rapidly converted to pastures and croplands from 1880 to 1930. Until the early 1990s, prevailing agricultural practices followed a low-external-input model with low contamination related to nutrient overloading (Viglizzo et al. 2001). However, in the last decades, the region has undergone a process of agricultural intensification with the introduction of genetically modified crops, the displacement of cattle breeding to marginal zones (mostly riparian areas), and the massive use of pesticides (Manuel-Navarrete et al. 2009;Song et al. 2021). For more details on the characteristics of the study region, please see Feijoó and Lombardo (2007).
Pampean streams are characterized by low flow and water velocity, and abundant macrophyte communities. They usually originate in vegetated depressions and are mainly fed by groundwater. Streambeds have hard superficial tuff layers rich in calcium, covered by fine sediments (silt and clay), without stones or pebbles (Feijoó and Lombardo 2007). Water nutrient levels are generally high when compared with streams from other regions of the world (mean annual concentration of phosphate and nitrate are > 0.15 mg/L and > 1.5 mg/L, respectively) (Feijoó and Lombardo 2007;Amuchástegui et al. 2016). In previous studies, nitrate concentration was related to agricultural practices, while phosphorus concentration has been associated with phosphorusrich parent materials of the region (Feijoó and Lombardo 2007;Amuchástegui et al. 2016;Torti and Portela 2020).

Field survey
To analyze changes in water quality in Pampean streams during the last years, we sampled 41 streams in Buenos Aires province between 2015 and 2016. These streams were previously sampled between 2003 and 2004 to determine baseline water quality for Pampean fluvial ecosystems. Our sampling covered streams from the four fluvial regions proposed by Frenguelli (1956) (Fig. 1): Region 1: Salado River and its tributaries. Region 2: Vallimanca Stream. Region 3: Tributaries of the Paraná and Río de la Plata rivers.
Region 4: Tributaries of the Atlantic Ocean.
In the previous study (Feijoó and Lombardo 2007), the streams were selected because they represent the whole range of variation in stream morphology, substrate, and flow of the region, and because they show low disturbance levels (Hughes et al. 1986). We used this approach because pristine watersheds cannot be found in the region after decades of agricultural activities. Criteria to select the streams included the absence of riparian forest or morphological channel alteration, the presence of natural vegetation or extensive cattle breeding in the proximity of the sampling point, and good water quality (dissolved oxygen concentration > 6 mg/L, ammonium levels < 0.1 mg/L) (Feijoó and Lombardo 2007). Sampling was done in May (autumn), November (spring), and February (summer) at baseflow conditions.
Between 2015 and 2016, we repeated the sampling performed in 2003/2004 at the same streams and seasons, using the same methodological approach. Briefly, we registered channel width and depth at each stream. We measured water velocity with a multiprobe anemometer (Schiltknecht Mini-Air20), and these data were then used to determine stream flow (in L/s) by the velocity-area method (Gordon et al. 1992). We also took pictures at each stream to register land use and other features in the surrounding landscape.
Photosynthetically active radiation (PAR) in air and water was measured with a radiometer LI-250 A + quantum sensor (Li-Cor Inc., Lincoln, Nebraska). We also recorded temperature, pH, conductivity, and dissolved oxygen (DO) concentration with a Hach multiprobe (HQ40D). We collected two water samples in 250 mL acid-washed polyethylene bottles at each stream. One sample was filtered through Whatman GF/C filters in the field and it was immediately preserved with 0.1% chloroform (Gardolinski et al. 2001) to determine nutrient concentration. The other water sample was not filtered and was then used to estimate alkalinity and chloride concentration. Samples were stored in an ice chest and transported to the laboratory for analysis.

Sample analysis
We estimated alkalinity within 12 h after collection. Nutrient samples were preserved at 4 °C and analyzed within 1 week from sampling following the methods of APHA (2017). We determined soluble reactive phosphorus (SRP) by the ascorbic acid method, and ammonium by the indophenol blue method. Nitrate and nitrite concentrations were measured with a FUTURA Autoanalyzer (Alliance Instruments,  Frenguelli (1956): 1 Salado River and its tributaries, 2 Vallimanca Stream, 3 tributaries of the Paraná and Río de la Plata rivers, and 4 tributaries of the Atlantic Ocean. Black circles: streams sampled in this study Frepillon, France) through a reaction with sulfanilamide with a previous Cu-Cd reduction for nitrate. In the unfiltered sample, we determined chloride concentration by the argentometric method and phenolphthalein alkalinity (carbonates and bicarbonates) (APHA 2017).

Watershed characterization
The catchment area was estimated upwards from the sampling point in each one of the 41 streams. We delimited catchment boundaries by ArcMap 10.3 (ArcGIS ® by Esri) software from digital elevation models (DEM). DEM was extracted from the Shuttle Radar Topography Mission (SRTM) data available in Global Land Cover Facility (GLCF; www. umd. edu).
Physiographic parameters (catchment area, mean slope, stream length, drainage density, and stream order) were determined using different ArcMap tools based on the generated DEMs, watershed boundaries, and the corresponding drainage network derived from them.
Land cover was determined using TM Landsat 5 images (available at the National Institute for Space Research site of Brazil; www. gov. br/ inpe/ pt-br) obtained between December 2010 and November 2011. Landsat images were projected to the Transverse Mercator Reference System and then georeferenced to the SRTM images using ArcMap 10.3 software. Images were geometrically corrected by applying the nearest neighbor method (mean square error for the transformation < 2 pixels; Chuvieco 1996).
Land use was defined in Landsat images at the 453 (RGB) multiband composition with the ArcView GIS 3.2 software considering several visual criteria (color, size, form, texture, etc.), temporal cover variation between images, and complementary information from Google Earth. Polygons from the same land use type were digitalized manually on the screen. We considered the following land use classes: cropland, cattle breeding, natural vegetation (grassland), and other covers, which included urban land (small villages), forest, bare land, mining, and ponds. Cattle breeding use was associated with larger geometric lots compared with agricultural ones. Livestock lots presented a variation of shades in the range of red, orange, yellow, green, and light blue, and included patches of black water and grassy or forest vegetation with irregular shapes. Natural vegetation cover was differentiated from cattle breeding use because the former showed irregular shapes that appear rough due to the change of tonalities (between reddish, dark orange, and brownish color) and because it was usually associated with water bodies or depressed areas. We calculated the percentage of each land use upstream of the sampling point. Land use was verified by visiting some of the studied catchments.
Soil type was determined at the soil order level with the soil taxonomy charts of the National Institute of Agricultural Technology (INTA; Argentina). Charts were used at different scales (1:50,000 and 1:500,000) according to the availability per basin and were projected and georeferenced following the same procedure for Landsat images. The percentage for each soil type (Alfisols, Entisols, Mollisols and Vertisols; Soil Taxonomy 1975) was calculated for each watershed.
The geology of each basin was determined from the 1999 Geological and Mining Map of the Province of Buenos Aires (Argentine Geological Mining Service-SEGEMAR) at a scale of 1:750,000. Geological categories present in each basin were digitalized manually according to the map references and the percentage of the area they occupy was calculated.

Phosphorus source apportionment using the SPARROW model
The SPAtially Referenced Regressions On Watershed attributes (SPARROW) model was used to further evaluate which features of the basins may determine SRP load in the streams. SPARROW is a hybrid mechanistic-statistical technique that estimates pollutant sources and describes contaminant transport within a detailed network of stream reaches and corresponding subwatersheds (Alexander et al. 2002a, b). The SPARROW model was expressed as (based on Alexander et al. 2002a, b;and Wellen et al. 2012): where the subscript i refers to each basin; µ i is the (log-transformed) mean annual load of SRP measured at sampling station i (kg/year); N is the total number of SRP sources and n is the index for each source; β n is the estimated contribution of source n; S n is the area (in km 2 ) of source n in the basin i; Z (delivery factor) represents landscape features (rain, runoff, etc.) that influence SRP transport from land to water and α is the associated coefficient; H s i is the fraction of SRP originating in the basin that remains at station i after nutrient attenuation processes (in-stream first-order loss function); H i P is the fraction of SRP originating in the basin remaining at station i as a function of first-order loss process modulated by pH (see below); ε i represents a random error, which is independent and equally distributed among all basins; and b i is a term to account for spatially correlated error (see below).
In-stream SRP loss process was expressed as a first-order loss function: where k s is the first-order loss coefficient (km −1 ) and L i is the reach length (km). SRP loss process depending on pH (H i P ) was defined as: where k p is the first-order loss coefficient and pH i is the mean pH at station i.
SPARROW was applied using data from 39 stream reaches sampled in 2015/2016 (streams A27 and C39 were not included due to anomalous data, see below). Mean annual load (µ i ) was estimated as the product of SRP concentration and flow at each stream. The areas contributing SRP to the streams were delineated considering different information: land use (cropland, cattle breeding, natural vegetation, others), soil types (Mollisols, Alfisols, Vertisols, and others), and three classes of geologies: (a) G2: deposits of gravel, sands and silts, piedmont levels. (b) Buenos Aires formation: loessoid silts, relativelyclayey, homogeneous, without stratification with concrete calcareous nodules, epigenic. (c) Other geologies.
The land-to-water delivery factor (Z in Eq. 1; m/year) was calculated as the area-weighted runoff produced at each basin. Since runoff may reflect the effect of climate, soil properties, lateral flow paths, and water storage on nutrient export from the land phase (Alexander et al. 2002a, b), it was used as a proxy to describe the main characteristics involved in the terrestrial transport of nutrients.
The model included an additional error term b that accounts for spatially correlated errors, potentially associated with an unexplained spatial pattern in the data (Qian et al. 2005;Wellen et al. 2012).
SPARROW was executed in the WinBugs 1.4.3 environment (WinBUGS 1996-2007: Imperial College andMRC, UK). We ran the model under different scenarios, considering different delineations of the land areas contributing SRP using the land use, soil type, and geology class information in combination or alone, and including (or not) SRP loss process associated with pH (H i P ). The calibration of the models was based on the minimization of the mean root square error between observed and estimated SRP loads in the different streams, using the Markov Chain-Monte Carlo (MCMC) methodology for the exploration of the parametric space in a Bayesian inference framework. Parameters α, β, k s k p , and b i (Eq. 1) were estimated by the model. Observed loads were estimated as the product of flow and SRP concentration at each stream.

Statistical analyses
Generalized linear models (GLM) were used to evaluate differences in-stream nutrient concentrations between years and among fluvial regions, along with possible continuous explanatory variables (PAR, land use, soil type, flow, etc.).
(3) H P i = exp −k p pH i L i GLM are an extension of linear models that allows to deal with nonnormal distribution of variables and the lack of homogeneity of variances among groups (Cayuela 2009). The best models were selected considering the Akaike information criterion (AIC). Two streams were not included in the GLMs because they showed very high levels of nutrients (stream C39) or an anomalous distribution of land use (stream A27). Differences in current velocity among periods were tested using a nonparametric test (Mann Whitney U-test) because this variable could not be normalized through transformation.

Visual changes at the local scale
We compared pictures taken in 2003/2004 with those taken in 2015/2016 to visually assess changes in channel morphology, riparian condition, and land use in the proximity of the sampling point. Five of the 41 streams showed morphological channel alterations and disconnection from the floodplain. We also observed afforestation by invasive or implanted trees in the riparian zone of five streams. Finally, land use changes (mainly the replacement of natural grassland and pastures by crops and the removal of riparian vegetation by cattle grazing) were detected in 11 streams (Fig. 2).

Geology and soil type
The dominant geology in all regions was Buenos Aires Formation. Mean basin cover of this Formation was > 75% in Salado, Vallimanca, and in Paraná and Río de la Plata regions, and similar to 50% in the Atlantic region. Other geologies were also well represented in specific basins in some regions: deposits of gravel, sands and silts, piedmont levels in Vallimanca and Atlantic regions, and other geologies in the Salado and Atlantic regions.
Among soil types, Molisols showed a higher mean basin cover in all regions (> 50%). Alfisols or Vertisols were dominant soil types in some basins in Salado region and Paraná and Río de la Plata region.

Land use and cover
Agriculture was the main land use (mean 69.5%), followed by cattle breeding (16.4%), natural vegetation (6.9%), and other uses (1.7%). Agricultural land use increased in three of the four hydrological regions between 2001/2002 and 2010/2011. The exception was the region of tributaries of the Atlantic Ocean, where the agricultural area remained almost the same. Cattle breeding use was reduced in all the regions, while natural vegetation showed a slight decrease in Vallimanca and Paraná and Río de la Plata regions, and a small increase in the tributaries of the Atlantic Ocean (Fig. 3).

Stream water characteristics
Sampled streams were characterized by oxygenated and alkaline waters, high conductivities and relatively elevated levels of SRP and nitrate, and low ammonia concentration (database available at https:// doi. org/ 10. 6073/ pasta/ 9db85 b43ed 25cc9 3f3aa 41eff 3e140 e9). Stream flows were generally low (mean 0.20 m/s), and showed no significant differences between 2003/2004 and 2015/2016 (p = 0.795), suggesting that hydrological conditions were similar in both periods. Moreover, flow differences could not be detected at a regional scale (Fig. 4). We used GLM to evaluate possible drivers of nutrient concentrations in stream water. The best model for nitrate (AIC 859.5; p < 0.000001) included region and year as categorical variables and flow as continuous variable, but only the variable region was significant. In the case of SRP, the best model (AIC −303.7; p < 0.00001) included year, region, and pH (Table 2), and all these predictors were significant. This indicates that SRP increased significantly between 2003/2004 and 2015/2016, varied among regions, and was related to pH levels. Surprisingly, when we included agricultural land use in the GLM as a continuous variable, it was not incorporated into the best model, suggesting that agriculture cover in the stream basin is not a significant driver of SRP concentration.
Given that pH seemed to explain the variability of SRP, we also explored whether agricultural activities have an impact on pH. The best GLM for pH included flow, region, and the interaction year/region, but not agriculture (AIC 308.7; p < 0.000001), suggesting a dilution effect of flow on pH. Therefore, pH may be mainly regulated by hydrological conditions and not by land use in the basin. Finally, we evaluated the possible influence of increasing SRP on dissolved oxygen concentration given the known relationship between phosphorus level and aquatic metabolism (Hoellein The GLM for dissolved oxygen (as a dependent variable), which considers region, year, and SRP as explanatory variables, selected the three variables (plus region/year) in the best model, being SRP the variable with higher significance (p < 0.00001). This suggests a relationship between dissolved oxygen and SRP concentrations, which is spatially (among regions) and temporally (between years) modulated ( Table 2).

Sources of phosphorus to streams with SPARROW
We first ran the model several times, considering different delineations of source land areas (alone or combined) and without including the process of SRP attenuation related to pH (H P ). We observed that the contribution of SRP (β in Eq. 1) did not differ among source areas, independently from how they were defined. This suggests that there are no zones within the basins that contribute differentially to SRP load due to different land uses, soil types, or geology classes, and that the whole basin area acts as a homogeneous source of SRP, at least at the level of precision that this analysis allows.
The model showed a good agreement between observed and estimated loads (R = 0.99), but b (correction factor of the spatial autocorrelation; Eq. 1) explained a relatively high proportion of the estimated load (mean 12.2 ± 11.5%).
In addition, the model overestimated the lower SRP loads (when b is negative) and underestimated the higher SRP loads (when b is positive). The parameter b was negatively related to pH (R 2 = 0.42), suggesting that high pH may reduce SRP load in the streams and that its inclusion in the model would compensate for the observed overestimation (Fig. 6).
For this reason, we ran the model including SRP loss due to pH (H P ) (and considering the whole basin area as an SRP source), and we observed that the mean proportion of estimated load by basin explained by b was reduced to 5.1 ± 7.7%. Consequently, results from SPARROW indicate that the whole basin acts as a source of SRP, with a contribution of 92.4 ± 7.2 kg SRP/km 2 year. The model also shows that pH is a main driver of SRP load variability in the streams (see Supplementary Information).

Discussion
After 12 years of agricultural intensification, we observed modifications at reach scale in several streams, including changes in channel morphology, afforestation, and land use modification. We also observed changes in water chemistry, with a decrease in pH levels and a rise in SRP concentration, especially in the southern fluvial regions (Vallimanca and tributaries of the Atlantic Ocean). On the other hand, mean nitrate concentration was related to differences in cropland use (in percentage) among regions, and it did not change significantly between 2003/2004 and 2015/2016.

Land use change
The replacement of natural grassland and pastures by crops and the removal of riparian vegetation by cattle grazing were the main impacts at the reach scale. The same pattern was observed at the basin scale. Except for the region of tributaries of the Atlantic Ocean, where mean agricultural cover was already high (> 80%) by 2001/2002, crop area increased while cattle breeding area and natural cover decreased in the basins by 2015/2016. In the Pampas, the cultivated area grew at the expense of areas used for either cattle breeding, or cattle and grain rotations, and this process was associated with a dramatic increase in the use of pesticides (Ferreyra 2006;Manuel-Navarrete et al. 2009). Glyphosate is the most applied pesticide, representing 62% of total commercialized herbicides in 2014 (CASAFE 2014). Under current practices, glyphosate application rates are higher than dissipation rates, leading to a widespread presence of this herbicide in soils and waters. For this reason, glyphosate and aminomethylphosphonic acid (AMPA, its main metabolite) are considered "pseudo-persistent" pollutants in the region (Primost et al. 2017).

Nutrient levels in the streams
The best GLM for nitrate included the variables "region," "year," and "flow," but only region was significant. This indicates that nitrate levels were significantly different among regions, increasing from regions with lower agricultural cover (Salado and Vallimanca) to regions with higher agricultural cover (tributaries of the Paraná and Río de la Plata rivers and the Atlantic Ocean) (Fig. 5). In addition, the inclusion of flow as predictor variable in the best model suggests an influence of hydrological conditions on nitrate levels. Baseflow in Pampean streams is maintained by groundwater flow, which provides nitrate to stream water; however, during storm events, nitrate is diluted by the input of runoff from surrounding lands (Sala et al. 1983;Feijoó et al. 2018). Agricultural cover was not included in the best GLM for nitrate, but it appears as a significant variable in the second best model. The relationship between nitrate concentration and agriculture is well known (see, for instance, Monteagudo et al. 2012;Shupe 2013), and it was also observed in a previous study in 23 of these streams (Amuchástegui et al. 2016). Therefore, our results suggest that nitrate concentration in Pampean streams depends on the agricultural cover and is modulated by hydrological conditions. However, we did not detect an increase in stream nitrate concentration between 2003/2004 and 2015/2016, despite the rise of fertilizer use. This may be explained by the observed negative nutrient balance in Pampean soils, which is produced by the export of phosphorus and nitrogen in grains that is not compensated by fertilization under current agricultural practices (Díaz de Astarloa and Pengue 2018). In the case of phosphorus, the best GLM model included "region," "year," and "pH," and all of them were significant. This means that SRP concentration differed among regions and years and that it was regulated by water pH. SRP concentration increased between 2003/2004and 2015. In 2003/2004, we observed clear differences among regions, with lower SRP values in the southern regions (Vallimanca and tributaries of the Atlantic Ocean) and higher values in northern regions (Salado and tributaries of Paraná and Río de la Plata rivers). However, these differences were attenuated in 2015/2016 mainly due to the significant increase in SRP concentration in the Vallimanca region. In Pampean streams, most phosphorus enters via superficial or subsurface flow (Rodríguez Castro et al. 2016;Feijoó et al. 2018); thus, we expected a clear association between SRP levels and land use. However, among all the possible explanatory variables considered (including flow, agricultural cover, carbonate concentration, etc.), the best GLM included pH. The proportion of arable and pasture lands has been associated with phosphorus concentration in fluvial systems (Ferrier et al. 2001;Jones et al. 2001;Gorgoglione et al. 2020); however, neither agricultural nor cattle breeding area appeared as significant variables in the GLMs tested for SRP. The role of pH as the main driver of SRP concentration in the streams is also supported by the results from the SPARROW model. The inclusion of a term representing in-stream SRP loss associated with pH (H i P ) in the model halved the spatial error. Interestingly, the model also suggests that the whole basin area supplies SRP to the streams, without the differential contribution of some areas due to land use, soil type, or geological classes. Additionally, the model calculated a contribution of 92.4 ± 7.2 kg SRP/km 2 year, which is similar to the mean SRP load estimated by Arreghini et al. (2005) in a Pampean stream after a storm event (115.5 kg SRP/km 2 year).
Abiotic mechanisms such as phosphorus adsorption to sediments and coprecipitation with calcium carbonates may be important drivers of SRP concentration, especially in calcareous streams (Jarvie et al. 2006;Jalali and Peikam 2013;Corman et al. 2016). But of the two processes, the former seems to be more relevant in Pampean streams (Martí et al. 2020). As it was mentioned before, streambeds of fluvial Pampean systems are covered by rich calcium tuff layers and fine sediments, where adsorption of phosphate ions on calcium carbonates leads to apatite formation (Millero et al. 2001). Acidification of stream water may displace the carbon-carbonate equilibrium, producing the dissolution of carbonates and the simultaneous SRP desorption from sediments (Reddy et al. 1999). In our study, the reduction of pH and carbonates concentration jointly with increasing SRP levels in the streams between 2003/2004 and 2015/2016 supports this explanation.

Stream acidification
If SRP concentration depends on pH, which factors could explain stream water acidification in the region? As in the case of SRP, the reduction of pH levels may not be related to land use and cover. The best GLM for pH included flow as significant variable, suggesting a dilution effect on pH during high-flow conditions as it was observed in other Pampean streams (Arreghini et al. 2007;Rodríguez Castro et al. 2016) Freshwater acidification has been generally attributed to the atmospheric deposition of strong acids (sulfur dioxide and nitrogen oxides) in industrialized regions of the world and to acid mine drainage (Hasler et al. 2018), although other factors have also been indicated. These include the rise of atmospheric CO 2 , conversion of terrestrial ecosystems to intensive agriculture (which increases the transport of labile carbon to freshwaters), enhanced terrestrial productivity and respiration because of global warming and growing atmospheric CO 2 (which intensifies CO 2 input by runoff), and nitrogen fertilization that produces elevated rates of nitrification (Jacinthe et al. 2004;Phillips et al. 2015;Hasler et al. 2016Hasler et al. , 2018Raymond and Hamilton 2018).
Carbon dioxide-induced acidification might be a possible explanation for the pH decline in Pampean streams, considering the regional extent of the acidification process. However, some punctual estimations we made using data from four tributaries of the Atlantic Oceans suggest that the rise of atmospheric CO 2 between 2003 and 2016 (28.43 ppm; Mauna Loa Observatory) could only explain a small fraction of the increase of CO 2 dissociation in water. Consequently, it seems that climate change is not the main driver of stream acidification.
Other possible explanations for stream acidification may be related to current agricultural practices in the region. Agricultural land use alters the amount and composition of dissolved organic carbon that is exported from watersheds (Jacinthe et al. 2004;Graeber et al. 2012;Messetta et al. 2023). The high input of labile carbon from croplands intensifies respiration by aquatic microbes, increasing CO 2 levels in water (Graeber et al. 2012;Hasler et al. 2018). A second possible explanation could be the use of nitrogenous fertilizers that induces soil acidification (particularly when nitrification rates are elevated), which in turn acidifies freshwaters via runoff (Nohrstedt 2001;Raymond and Hamilton 2018). Nitrification is a process mediated by bacteria that converts ammonium to nitrate in oxic environments, changing the acid-base balance of soils. Nitrification in soils was described as a significant factor for acidification in medium alkaline Danish streams (Rebsdorf et al. 1991). Acidification is also occurring in Pampean soils due to the removal of cations by crops and the use of acid fertilizers as ammonium salts or urea (Cruzate and Casas, 2003;Vázquez 2011;Vázquez et al. 2012). A third factor that could explain the acidification of Pampean streams is the widespread use of glyphosate. This herbicide is a strong acid (pH in 1% solution = 2; Water Quality Monsanto; https:// monsa nto-ag. co. uk/) that can lower the pH of the medium when in solution (Schweizer et al. 2019). Moreover, the glyphosate molecule contains one atom of phosphorus, and its application represents an input of anthropogenic phosphorus to freshwaters (Hébert et al. 2019). Some heterotrophic bacteria can degrade glyphosate, liberating phosphorus that is then used by aquatic biofilms (i.e., bacteria and fungi) (Vera et al. 2010;Carles et al. 2019;Carles and Artigas 2020). The use of glyphosate in Argentina, which was marginal in the early 1990s (2.4 million liters by 1991), raise from 220 to 300 million liters between 2003 and 2016 (Pengue 2016), and it has been detected in groundwater, shallow lakes, rivers, streams, and even rainfall in the Pampean region (Caprile et al. 2017;Primost et al. 2017;Vera-Candioti et al. 2021;Castro Berman et al. 2022). The combined effect of glyphosate as an acidifying agent and source of phosphorus could explain the lowering of pH and the rise of SRP in Pampean streams between 2003/2004 and 2015/2016. It seems plausible that all these factors associated with agricultural practices may contribute, to a greater or lesser extent each one, to the acidification of Pampean streams.

Conclusions
We hypothesized that dissolved nitrogen and phosphorus levels are strongly linked to land use in the basins. However, we could not associate agricultural cover with SRP concentration in the studied streams, while for nitrate we detected an influence of both hydrological conditions and land use. Most importantly, streams of the Pampean region are undergoing two coupled processes: eutrophication due to the increase of SRP concentration, and acidification. Our results indicate that these changes are not related to the expansion of cropland area in the basins, which was not substantial after 12 years, suggesting that factors related to changes in agricultural practices (such as the large increase in pesticides consumption) may play a role. This is a concerning point because if the process of acidification continues, displacement on the carbonate equilibrium may imply enough loss of carbonates to have calcite redissolution and release of phosphorus from the sediments, further intensifying water eutrophication.
In addition, the scenario of climate change for the region predicts an increase in precipitation (Barros et al. 2014). As most dissolved phosphorus enters the streams via superficial or subsurface runoff, this may also contribute to a higher phosphorus load to the streams. Consequently, it is mandatory to apply management practices that attenuate the input of phosphorus, as well as to adopt new agricultural paradigms to reduce the use of agrochemicals in Pampean basins.