Our analysis focuses in the Eastern Plains of Uruguay, the Southeast portion of the Rio de la Plata Grasslands (Fig. 1) (Soriano, 1991). The Eastern Plains have an extension of 743,600 ha covering a latitudinal range that goes from 32°28′S to 33°54′S. The average annual temperature is around 16.4 °C and the annual accumulated rainfall is approximately 1200 mm (INUMET, 2019). This area is traditionally devoted to rice cultivation as well as to extensive cattle raising (MGAP and DIEA, 2011). The predominant natural ecosystems are grasslands, wetlands and palm groves savannas (Achkar et al., 2012).

We carried out the LULC map in two stages. First, we made a supervised classification of a Landsat 8 image from 2016 that discriminates among five categories (Section “Supervised classification”). Second, we made decision-tree classifications of two categories (bare soil and “other” uses) of each of the best images available from 2007 to 2016, in order to generate a multi-temporal cropland mask with all those pixels of the study area that had agricultural use at least on once in the mentioned period of time (Section “Cropland mask”). The final LULC map was generated by superimposing the cropland mask on the 2016 classification (Section “Final land use/land cover map”).

The spatial variability of landscape patterns was explored through the intersection of LULC map with a 10 × 10 km cell size grid, each resulting cell was considered the “landscape unit” of the analysis. The study area comprises 74 cells that have at least 70% of its surface classified on the map (Fig. 1). For each cell, we quantified the percentage of grassland, wetland and agriculture use (“Percentage of landscape”, PL) and estimated grassland fragmentation from three widely used metrics that can give good reference about grassland spatial structure and configuration, i.e., the Number of patches (NP), the Effective Mesh Size (EMS) and the mean Nearest Neighbor Distance (NND) (Table 1) (Jaeger, 2000; McGarigal and Marks, 1995). For the calculation of these metrics, we converted the LULC map from raster to vector format and eliminated all plots corresponding to urban, water bodies, and all those with an area smaller than 1,8 ha (equivalent to 20 Landsat pixels). The number of grassland patches, area and Euclidean distance were obtained with ArcGIS v.10.3. and the metrics were calculated manually in an excel spreadsheet. The PL, NP and NND are easily interpretable attributes that describe simple structural features of the landscape (Jaeger, 2000). However, they have the disadvantage of not necessarily reflecting the fragmentation processes, since two different landscape units could have the same PL of some class, but have different division or isolation degree. The EMS, on the contrary, presents a more linear response to the different stages of the fragmentation process, even when fragmentation is not causing a marked decrease in habitat, like when the landscape is partially dissected by roads or other human infrastructures (Jaeger, 2000). It is calculated as an averaged size of the focus class weighted according to the class individual patch sizes. We made a mathematical modification of the EMS proposed by Moser et al. (2007) to eliminate what the author calls the “edge problem”, which arises when considering the artificial edges that are generated in the limits of the landscape units.

To analyze the biophysical and human controls of fragmentation variability, we modeled natural grassland fragmentation considering the PL, EMS and the NND as dependent variables; and 13 independent variables (Table 2) from five datasets related to rice crop development aptitude. This comprised (Table 2): (1) anthropogenic variables, particularly the road network density (since a greater amount of roads tends to facilitates the human activities development) and cadastral parcels size (since rice activity tends to occupy the largest ones, facilitating logistic activities); (2) edaphic variables, like soil type (particularly Planosols percentage and Planosols and/or Gleysols percentage, which are the optimal soil types for rice activity in the region); (3) hydrological variables, particularly hydrological network density, water courses surfaces and edges density (due to the high water requirements of rice crops); (4) topographic variables, particularly ground elevation and elevation variations (since low and flat lands facilitates the necessary flooding and drainage process of rice crop); and (5) spatial location variables (latitude and longitude) as aggregation indicators which may be associated with the presence of some industrial foci (Table 2). The set of independent variables has different spatial resolutions (ranging from pixels of 900 m2 in the digital elevation model, to irregular units with an approximate median value of 1.500.000 m2 in the soil types cartography) so they were summarised for each grid cell of the previous section, as shown in Table 2. Correlation analysis between all explanatory variables was performed by means of Pearson's coefficients. Finally, we explored the dependences of landscape metrics to the independent variables (biophysical and anthropogenic) by a stepwise multiple linear regression model. High correlated variables (significance value upper to 0.05 and correlation coefficient value equal or upper to |0.6|) were excluded.

Spatially and conceptually detailed data acquisition results essential to quantify subtle land cover changes in highly transformed landscapes with high biodiversity potential (Mallinis et al., 2011). Despite the fact that grasslands are among the most threatened and transformed biomes in the world, and the large-scale development of rice agriculture in the Eastern Plains, most of previous maps covering the study area did not attempt to discriminate between natural grasslands and other forage resources (Baeza and Paruelo, 2020; Graesser et al., 2015; Volante et al., 2015; Blanco et al., 2013; Clark et al., 2012; Friedl et al., 2010; Hansen et al., 2000). This is due to different objectives and scale analysis of these works, but also to the difficulty of discriminating natural grasslands from other herbaceous plant formations with satellite images. Mapping natural grasslands in the highly agriculturalised context of the Eastern Plains is not an easy task with the usual remote sensing techniques. Rotating system of rice production in Uruguay generates a wide variety of post-agricultural stages of vegetation as a consequence of different management practices, from sown pastures to different successional states of spontaneous vegetation growing after rice harvesting (i.e. crop weeds and other ruderal species). In spite of their structural differences (composition, diversity, percentage of alien species, etc.) these vegetation stages have similar phenological behavior to that of natural grassland which generates an overlapping of spectral signals and impedes natural grassland discrimination (Baeza et al., 2014). In our work these difficulties were solved with a simple and replicable technique that discriminates natural grasslands based on a combination of a broad classification of forage resources and the identification of plots that had agriculture in the recent past. The reconstruction of recent agricultural history has proven to complement the information provided by the current LULC classifications (Schmidt et al., 2016; Yan and Roy, 2014; Maxwell and Sylvester, 2012). Considering that paddocks dedicated to rice cultivation remain with very low or nil vegetation coverage most of the year, the agricultural history of the previous ten years was reconstructed from the detection of bare soil. Contrary to what happens with grasslands, bare soil before planting the crop or after harvesting has a very clear spectral signal, useful to discriminate areas that have had agriculture at some point and that allowed us to exclude them from grassland class in the final map. In a previous work, Maxwell and Sylvester (2012) combine intra- and inter-annual time series with spectral and temporal metrics to detect patterns in permanent cultivated lands and thus achieve accurate classifications. Other works combine time series with object-based analysis to generate automatic computational methodologies to distinguish cultivated areas (Schmidt et al., 2016; Yan and Roy, 2014). These techniques could not be applied in the study area due to the high clouds cover (even with years without a single useful image). Recent regional initiatives such Mapbiomas (Souza et al., 2020), particularly MapBiomas Pampa (https://pampa.mapbiomas.org/) has implemented the available technology to discriminate natural grasslands for other LULC. While their results remain to be validated, they are very similar to those found in our work.

Our results reveal how agricultural transformation has fragmented the natural grassland in the Uruguayan Eastern Plains, a fact observed in many regions of the Rio de la Plata Grasslands (Baeza and Paruelo, 2020; Clark et al., 2012; Baldi and Paruelo, 2008; Baldi et al., 2006; Guerschman et al., 2003). Only a 21% of the study area is covered by natural grassland remnants, which places the Eastern Plains as the most transformed region in the country so far. For comparison purposes, Baeza et al. (2019) mapped the natural grassland communities of the main livestock areas of Uruguay and reported that the percentage of natural grassland varied according to the region considered, from 39% in the South Central Region and 75% in the Basaltic Region (Baeza et al., 2019).

According to the conceptual scheme proposed by Forman and Gordón (1986) and modified by Jaeger (2000), natural grasslands in the Eastern Plains are in the most advanced stage of fragmentation process (i.e., shrinkage or attrition), being poorly represented in the landscape, and their remnants are highly divided into small patches. At regional level, Baldi and Paruelo, 2008 analyzed different regions of the Rio de la Plata Grasslands through the EMS and reported the most critical scenarios in the Flat Inland Pampa and Rolling Pampa, with values very close to those obtained in our work (EMS: 150 and 334 ha, respectively). The results obtained for the Eastern Plains are worrying, the degree of natural grasslands fragmentation is comparable to that of the historically most cultivated and transformed regions of the Rio de la Plata Grasslands (Baeza and Paruelo, 2018, 2020; Paruelo et al., 2005). These results contrast with those of Mittermeier et al. (2003), which consider the study area as one of the last wild areas in the world. Differences may be due to an increase in agricultural area in the first decades of the 21 st century and/or differences in the conceptual resolution in their base maps, which do not discriminate natural grasslands from other forage resources. The evidence found in the Eastern Plains does not differ from the critical scenarios reported for other temperate grasslands in the world (Watson et al., 2016; Bond and Parr, 2010; Neke and Du Plessis, 2004) which once again places them as one of the most threatened natural biomes.

The results obtained from the multiple linear regression analysis support the hypothesis that factors that define the establishment of rice activity determine the spatial variation of grassland fragmentation in the Eastern Plains. The same has already been indicated in numerous works for other crops in Uruguay (Arbeletche et al., 2012; Arbeletche and Gutiérrez, 2010; Baeza et al., 2014; Hoffman et al., 2013; Paruelo et al., 2006) and in the region (Baldi et al., 2006; Baldi and Paruelo, 2008; Clark et al., 2012; Guerschman, 2005; Viglizzo et al., 2001). But unlike most of these studies, our results quantitatively demonstrate the relationship between agricultural suitability and the loss and fragmentation of natural covers. In our work, three variables turned out to be the main drivers of grassland replacement in the Easter Plains: the cadastral parcels density (therefore their size), the road network density and the topographical variation. Interestingly, topography explained 41% of the spatial variability of the EMS in the study area.

The places where natural grassland is more fragmented are those with flatter topography, with better suitability for rice cultivation due to easier access, management and final disposal of water (Battello et al., 2013; DIEA and MGAP, 2003, 2017). The areas with the lower percentage of grassland are concentrated in the North and Centre of the study region. In these places the average size of cadastral parcel is higher than in the South, which explains their lower density. A recent study shows that 80% of the rice planted in the Eastern Plains is on leased land, being more profitable the use of larger parcels (DIEA and MGAP, 2017). Since the main agro-industrial centers of rice are located there (Alegre et al., 2015; DIEA and MGAP, 2017), these areas have denser road networks that facilitate the movement of both production and machinery needed to implement the crop. In addition, the topographical variation is lower than in the South. Isolation of grassland patches is associated with similar variables; the NND variation is also explained by the size of cadastral parcels (finding more isolated patches in areas with larger parcels) and by the position in the landscape (closer to rice logistic foci). The density of hydrological network also explains the degree of isolation of grassland remnants, as more isolated patches are found in areas where there is a greater presence of water courses, probably due to the high water requirements of rice cultivation.

Several studies have affirmed that landscape composition in the Rio de la Plata Grasslands, and indirectly grassland fragmentation, are not random processes. There is a complex interaction between different variables that restrict agriculture expansion, relegating grassland to places where unfavorable characteristics for the development of this activity are combined (Guerschman et al., 2003; Viglizzo et al., 2001). Baldi et al. (2006) found that soil drainage explained the spatial distribution of agriculture and grassland surfaces in other regions of the Rio de la Plata Grasslands. This was very logical given that the greatest crops restrictions in their study sites are water stress and water excess (Hall et al., 1992). Unlike what was initially expected, the soil type does not explain the spatial variation of the fragmentation in our work, probably due to the scale of analysis (i.e., 10 × 10 km units). However, at smaller spatial scales, paddock edaphic characteristics are among of the main determinants of rice yield in the Eastern Plains (Tseng et al., 2021) and therefore in the choice of crop location. So probably, at a finer scale of analysis, other variables linked to suitability for rice crop will become more relevant in explaining spatial variation in fragmentation.

Our work highlights the high degree of grasslands lost in the Eastern Plains due to the replacement by crops or implanted pastures, a situation that is repeated in temperate and subtropical grasslands around the world (Watson et al., 2016; Bond and Parr, 2010; Neke and Du Plessis, 2004). This replacement of natural habitats and the advanced fragmentation degree of grassland remnants imply a decline in native species diversity (Deák et al., 2021; Lampinen et al., 2018; Staude et al., 2018; Fahrig, 2003; Sala et al., 2000; Pimm and Raven, 2000). As patch size decreases, species richness may be reduced (Bennett and Saunders, 2010) because patches have become smaller than the minimum area required to support populations or individuals of species with larger range requirements (Fahrig, 2003). In addition, the loss of habitat corridors, increases the patches isolation, so a reduction in species diversity could be expected (Deák et al., 2021; Lampinen et al., 2018; Fahrig, 2003).

The current scenario of high grassland fragmentation in the Eastern Plains makes it of great importance to develop prompt conservation measures for the remaining natural areas, as well as restoration strategies for degraded sites. These actions are in line with international commitments taken by the country, such as the RAMSAR agreement for wetland conservation and the Aichi Biodiversity Targets for 2020 (https://www.cbd.int/aichi-targets/). Due to be a region where conservation and productive interests converge (Brazeiro et al., 2015; DIEA and MGAP, 2003) the compliance with the Aichi Targets, in particular point 7—which promotes the sustainable management of areas of agricultural use, guaranteeing the conservation of biological diversity – and point 15—which commits to increasing the resilience of ecosystems through conservation and restoration – assume particular relevance.

According Buisson et al. (2021), a critical step in the restoration process of grassland ecosystems is to have reliable information on the distribution of different grassland types, particularly pristine ones for use as reference ecosystems. The current LULC map generated here will allow us to describe and characterize the grassland plant communities in the remnant sites, their states and possible transitions to others, as well as evaluate their potential for post-agricultural restoration. The results of the landscape metrics obtained in our work also allow the recognition of sites with a higher and lower degree of fragmentation, which will be extremely useful for defining priority areas for conservation. In fragmented landscapes, conservation efforts have typically focused on the preservation of remaining large habitat patches that are intact and well connected (Fischer et al., 2009), but small patches may play a significant role in conserving remaining vegetation (Fischer and Lindenmayer, 2002). Particularly, plant communities present in grasslands remnants have proved to be an important source of natural diversity (Burkart et al., 2011). Increasing the quality of a focal habitat patch by ensuring proper management strategies can effectively compensate the losses in the total area (Otsu et al., 2017). Appropriate management of natural remnants (such as field margins, river and stream edges, and areas without optimal characteristics for cultivation) or suitable novel habitat patches (such as clearings under power lines, grasslands at dams) can also increase the amount of grassland and landscape connectivity (Lampinen et al., 2018; Bátori et al., 2020; Baum et al., 2004; Uezu et al., 2008; Tulloch et al., 2015).

Finally it is important to note that not all environmental changes can be detected through land cover transformation (Li and Wu, 2004; O’Neill et al., 1997). It is necessary to integrate this information with other field studies and link them with functional landscape properties to assess the state of its ecological resources (Li and Wu, 2004) and to devise management strategies to prevent further grasslands loss and promote their conservation.