Spatial and Temporal Dynamics of the Vegetation Cover from the Bijagual Massif, Boyacá, Colombia, during the 1986–2021 Period

1. Introduction

High mountain ecosystems in Colombia are characterized by a degree of endemicity and richness of taxonomic entities [1,2] due to the outcome of landscape and topographic evolution, as well as climatic modeling and ecological interaction [3,4]. For example, plant physiognomy, functionality, and diversity have allowed the identification of ecological associations related to different altitudinal belts, such as high Andean forest and paramo [5,6,7,8]. In fact, these ecological associations have been considered strategic [9] because they provide ecosystem services such as food and climate, air and soil regulation, and water capacity filtration and concentration [10,11]. However, for example, in the Andean forest, human transformations have shaped the native forest landscapes into a matrix of forest patches, cultivated and grazing areas, and socio-ecosystem interactions or relationships [12,13].

Due to the socio-ecosystem relationship that occurred in this region since the 19^th century [14,15,16], the vegetation cover changes over time can be identified and, hence, produce accurate information on the distribution and rates of vegetation cover chances, particularly regarding forest and paramo, which are essential for conservation strategies [17,18]. Change detection at the landscape scale can be identified by analyzing the spatiotemporal effect of vegetation reflectance [19,20,21]. Spectral variations model the vegetation reflectance using combinations of satellite image bands [17,21]. Furthermore, because Landsat satellite images can be obtained from 1970, the historical changes in plant covers can be identified over time [20,21].

In Colombia, some Andean forest and paramo areas are delimited by administrative units called complexes [22], which include the Bijagual Massif under the Tota-Bijagual-Mamapacha complex. The land use and forest deforestation history of Bijagual has been related to intensive and extensive agricultural and livestock activities [23,24,25]. However, the magnitude of human activities under the spatiotemporal effect of vegetation reflectance and, therefore, the historical, spatial, and temporal dynamics of fragmentation and ecological connectivity are unknown [26,27].

Accordingly, this work studied the spatial and temporal dynamics of the Bijagual Massif vegetation in Colombia to measure the ecological integrity of the natural area replaced by excessive human activities [25]. The supervised classification method with satellite images for a 35-year period from 1986 to 2021 was employed to obtain diversity, composition, and ecological connectivity information through landscape-scale information and coverage monitoring to suggest natural conservation strategies and territorial planning development. The results and their analyses are explained in three sections: 1. Land covers - field verification, 2. Classification accuracy, and 3. Multitemporal dimension.

2. Results

2.1. Land covers – field verification

Field verification allowed the identification of five cover types spectrally assigned to the Corine Land Cover categories (CLC; IDEAM [28]), as follows: 1. High open forest (Hof): includes the high Andean forest formation dominated by trees > 10 m tall; 2. Dense grasslands of non-wooded mainland (Dgnm): belongs to the paramo ecosystem characterized by scrubs, pajonal-chuscal, frailejonal, and sphagnal type communities; 3. Mosaic of pastures and crops (Mpc): lands whose use-occupation are pastures and crops, particularly potatoes; 4. Lagoons, lakes, and natural swamps (Lag): natural water surfaces that, for Bijagual, correspond to the La Calderona and El Pato lagoons; 5. Bare and degraded lands (Bdl): deforested area (Figure 1).

The vegetation of the Bijagual Massif belongs to Andean mountain formations that characterize the high mountain ecosystems in Colombia. According to Gil-Leguizamón et al. [29,30], the strip of high Andean forest (Hof) shows a wide species turnover at short distances, with trees and shrubs of encenillos (Weinmannia tomentosa L. f., W. rollottii Killip, W. fagaroides Kunth, W. balbisiana Kunth, and W. reticulata Ruiz & Pav.), gaque (Clusia alata Planch. & Triana, C. elliptica Kunth, and C. multiflora Kunth), susque (Brunellia propinqua Kunth, and B. cf. comocladifolia Bonpl.), garrocho (Viburnum tinoides L. f., and V. triphyllum Benth.), uva camarona (Macleania rupestris (Kunth) A.C. Sm.), tunos, charne, and siete cueros (Miconia cataractae Triana, M. cundinamarcensis Wurdack, M. jahnii Pittier, M. ligustrina (Sm.) Triana, M. theizans (Bonpl.) Cogn., Axinaea scutigera Triana, Bucquetia glutinosa (L. f.) DC., and Tibouchina grossa (L. f.) Cogn.), and manos de oso (Oreopanax bogotensis Cuatrec., and O. mutisianus (Kunth) Decne. & Planch.), among other species identified in field verification visits and corroborated with the collections of the UPTC Herbarium of Universidad Pedagógica y Tecnológica de Colombia.

In the dense grassland of non-wooded mainland or paramo coverage (Dgnm), there are associations of scrubs, pajonales-chuscales, frailejonales, and sphagnals included in the Linochilo phylicoidis-Arcytophyllion nitidi Rangel & Ariza 2000 nom. mut. nov. (all.) alliance described by Pinto Zárate [31] for the Eastern Cordillera of Colombia. Shrub communities include species such as Bejaria resinosa Mutis ex. L. f., Escallonia myrtilloides L. f., Miconia cleefii L. Uribe, Geissanthus andinus Mez, Ageratina theaefolia (Benth.) R.M. King & H. Rob., Diplostephium phylicoides (Kunth) Wedd., Pentacalia pulchella (Kunth) Cuatrec., Monnina salicifolia Ruiz & Pav., Hesperomeles obtusifolia (Pers.) Lindl., and Symplocos theiformis (L. f.) Oken., among others. The pajonal-chuscal formations are comprised of Agrostis perennans (Walter) Tuck., Anthoxanthum odoratum L., Calamagrostis effusa (Kunth) Steud., and Chusquea tessellata Munro., among others. Frailejonal formations include Espeletia argentea Bonpl., E. boyacensis Cuatrec., and E. murilloi Cuatrec., and the sphagnals formation comprises a complex of Sphagnum species [29,32].

Endemic and endangered species are found in the Bijagual Massif, including forest formation species, such as Greigia stenolepis L.B. Sm., Hieronyma rufa P. Franco, Tillandsia pallescens Betancur & García Nestor, Symplocos venulosa Cuatrec., Diplostephium oblongifolium Cuatrec., and Dunalia trianaei Dammer [30]. Furthermore, Espeletia cayetana (Cuatrec.) Cuatrec., E. murilloi Cuatrec., and Puya goudotiana Mez. [32] were found in the paramo ecosystem. In addition, species widely distributed in the high Andean forests and paramos in Boyacá were recorded, including angelitos (Monochaetum myrtoideum Naudin, M. uribei Wurdack, and M. meridense Naudin), chites (Hypericum laricifolium Juss., H. mexicanum L., H. myricariifolium Hieron., and H. lycopodioides Triana & Planch.), and cortadera (Cortaderia nitida (Kunth) Pilg.) [33,34].

The species composition in the Bijagual Massif is related to a high biodiversity potential [29]. Nevertheless, vegetation is limited by the pastures and crops cover (Mpc) mosaics extending from south to north in the study area. This cover type is dominated by potato (Solanum tuberosum L.), mora (Rubus urticifolius Poir.), and beans (Phaseolus vulgaris L.), as well as grazing cattle and sheep, considered potential threats to the ecological integrity of the verified ecosystems [35].

2.2. Classification accuracy

During the field verification, 276 control points were generated. Producer and user accuracy fluctuated between 98 and 100% of the covers. The Kappa coefficient established a thematic reliability of 95% (gray diagonal boxes; Table 1).

2.3. Multitemporal dimension

The Bijagual Massif has a total landscape area (TA) of 8,574.1 ha. In 1986, Dgnm (paramo) occupied 42.6% of the total area, followed by Mpc (32.8%) and Hof (24.5%). In 2000, the spatial composition was maintained; however, the productive systems and the forest increased (43.7% and 28.1%, respectively), while the paramo decreased (28%). In 2021, Bijagual showed paramo extensions restricted to the northeast in the municipality of Viracacha and isolated patches in the center and south (2,160 ha; 25.2% of the total area). Mpc occupies 52.9% (4,534.1 ha) and Hof 21.7% (1,861.8 ha; Table 2; Figure 1).

Between 1986 and 2000, Dgnm lost 1,248 ha (34%), while Mpc and Hof increased 936.6 ha (25%) and 309.2 ha (13%), respectively. From 2000 to 2021, Dgnm and Hof lost area (241 ha, equivalent to 10% and 548.7 ha or 23%, respectively). Furthermore, Mpc increased (785.6 ha or 17%). In 35 years, paramo and forest areas decreased (1,489 ha or 41% and 239.5 ha or 11%, respectively), while pastures and crops increased (1,722.2 ha or 38%) and occupied more than 50% of the land extension of the Bijagual Massif (Table 2; Figure 1).

2.3.1. Spatiotemporal changes 1986-2000

From 2,101.3 ha Hof registered in the year 1986, 130.6 ha changed to Dgnm (6.2%) and 92.1 ha to Mpc (4.4%) until 2000. Likewise, Dgnm was replaced in 1,182.8 ha for Mpc (32.4%) and 342.5 ha for Hof (9.4%). Similar to Hof, the coverage of Mpc changed by 191.4 ha to Hof (6.8%), 146.8 ha to Dgnm (5.2%), and 2,473.6 ha (88%) remained unchanged (Table 3). Between 1986 and 2000, of the total area of 8,574.1 ha, 6,484.3 ha maintained the same coverage (75.6%), and 2,089.9 ha changed to another cover type (24.4%).

2.3.2. Spatiotemporal changes 2000-2021

From the 2,410.5 ha Hof registered, 447.3 ha changed to Dgnm (18.6%) and 521.9 ha to Mpc (21.7%). Further, from the 2,401 ha Dgnm recorded, 257.9 ha changed to Hof (10.7%) and 740.6 ha to Mpc (30.8%). In addition, 3,748.5 ha of Mpc changed during the 2000-2021 period to another cover type; 161.6 ha changed to Hof (4.3%), 310.3 ha to Dgnm (8.3%), and solely 5.1 ha to Bdl (0.1%) (Table 3). Between 2000 and 2021, only 6,128.6 ha maintained its cover (71.5%), and 2,445.5 ha changed from one to another (28.5%).

2.3.3. Diversity

In 1986, 76 patches (TNP) corresponded to natural or transformed areas, with a density (PD) of 0.88 patches/ha. In 2000, 75 patches (PD: 0.87 patches/ha) were recorded, and in 2021, 396 patches (PD: 4.6 patches/ha) were found.

2.3.4. Composition

From the total number of patches in 1986, 25 belong to Mpc (PD: 0.3 patches/ha), but these decreased in 2000 (21 patches/PD: 0.24 patches/ha) and increased by 2021 (74/PD: 0.86). For Mpc in 1986, the largest patch index (LPI) was 18.4%, rising to 22.7% in 2000 and 48.7% in 2021. In Dgnm, the number of patches between 1986 and 2000 increased (from 25/PD: 0.3 to 31/PD: 0.36); however, the number of patches was higher in 2021 (168/PD: 1.9). Conversely, the LPI values of this coverage decreased in 35 years (1986: 21%; 2000: 8.1%; 2021: 6.7%). The abundance of Hof patches decreased between 1986 and 2000 (from 23/PD: 0.25 to 20/PD: 0.23) but increased in 2021 (150/PD: 1.7). The LPI index for Hof shows a decrease over time (1986: 8.5%; 2000: 8.2%; 2021: 6.5%; Figure 2).

2.3.5. Connectivity

Of the three fractions of the dPCk probability index, only dPCintra and dPCflux contribute to connectivity, while the contribution of the dPCconnector fraction is limited (Figure 3-A).

In 1986 and 2000, the Hof and Dgnm covers showed patches with habitat quality and availability for biota (connectivity within the patch) and act as dispersal flows (connectivity between patches); however, these ecological characteristics decreased drastically in 2021 for the two covers. Both Hof and Dgnm patches provide different connectivity values, evidenced in the dispersion of the standard deviation data concerning the average of each fraction by cover type and time (Figure 3-A).

Habitat availability contributions and flow between patches in 1986 and 2021 were low (dPCintra <10 and dPCflux <5; Figure 3-B). In the spatiotemporal dynamics, changes in the connectivity contributions of each patch were identified for the dPCintra and dPCflux fractions. In this way, in 1986, Hof showed two patches that contribute to habitat quality and availability. In 2000, one patch contributed to habitat quality and availability, two contributed to habitat quality and ecological flow, and the other two to ecological flow (Figure 3B). In 2021, only one patch contributed to habitat quality and availability and five to ecological flow (Figure 3B,C). For its part, Dgnm in 1986 showed one patch with contributions to habitat quality, habitat availability, and ecological flow (or dispersion), another with contributions to habitat quality and availability, and solely one patch with contributions to ecological flow. In 2000, only one patch registered contributions to habitat quality and availability (Figure 3B), and in 2021, one patch contributed to habitat quality, habitat availability, and ecological flow, and another solely to ecological flow (Figure 3B,C).

3. Discussion

3.1. Land cover – field verification

The Bijagual Massif corresponds to one of the six high mountain natural areas in Boyacá with broad floristic diversity [22,34]. The field verification process and the consultation of herbarium biological collections [30,32] allowed estimating a richness of up to 327 species for the Hof cover and up to 120 species for Dgnm. The results of the spatiotemporal dynamics showed that this diversity is vulnerable to the effects of the pasture and crop matrix since it has caused modifications in the vegetation [36]. According to Gil-Leguizamón et al. [29,30], the presence and dominance of some species of Chusquea Kunth in forest interiors and of Cenchrus L., Brachiaria (Trin.) Griseb., and Holcus lanatus L. in paramos indicate that the composition and physiognomy of these covers have been modified, enhancing the loss of richness and forming more heterogeneous plant communities. Even so, the research carried out by Gil-Leguizamón et al. [29,30] and Carrillo et al. [32] relate some species of ecological importance, including trees, shrubs, rosettes, and grasses that characterize the physiognomy of the Andean forest, such as C. multiflora, C. elliptica, C. alata, W. rollottii, B. comocladifolia, V. triphyllum, and E. myrtilloides, and then the species Espeletia murilloi Cuatrec., Paepalanthus columbiensis Ruhland, Puya goudotiana Mez, Hypericum lycopodioides Triana & Planch., Blechnum auratum (Fée) R. M. Tryon & Stolze, Calamagrostis effusa (Kunth) Steud., and Rhynchospora ruiziana Boeckeler typical of paramos. These records agree with the altitudinal delimitation of the vegetation proposed by Cuatrecasas [37].

Research is required for these groups to allow inferring biogeographic relationships, particularly of the Andean forest, which, through distribution records, are not only indicators that corroborate land cover but also allow interpreting the influence of the environmental gradients in the distribution patterns of the flora, its endemisms, and affinities with other biogeographic units, as has been done for the non-vascular flora of Bijagual [38].

Besides generating thematic reliability in land cover layers, field verifications allow identifying species that, due to their phenological, reproductive, and dispersal characteristics, can become a potential for ecological restoration processes in Bijagaul [35,39]. Within these, Myrcianthes rhopaloides (Kunth) McVaugh, Pentacalia pulchella (Kunth) Cuatrec., Diplostephium floribundum (Benth.) Wedd., Vallea stipularis L. f., Weinmannia fagaroides Kunth, Macleania rupestris (Kunth) A.C. Sm., Hypericum lycopodioides Triana & Planch., Bucquetia glutinosa (L. f.) DC and some species of Viburnum L. and Cestrum L. are highlighted.

3.2. Classification accuracy

From the 276 control points obtained during the field verification and the overlap of these points with the land cover raster layer (spectral response obtained from the supervised classification for the 2021 image), a thematic reliability of 95% (Kappa 0.95) was obtained. According to cartographic quality standards, the data are associated with a "very good spatial concordance" [40]. Of these 276 points, nine did not coincide with their respective classification cover due to the similarity of the digital value of a pixel in the image associated with another cover. This was evidenced for Hof with six points corresponding to forest cover in the field verification and that in the supervised classification algorithm were assigned to Dgnm (4 points) and Mpc (2 points). A similar case was obtained for Lag with only two points and Dgnm with one (Table 1). The remaining 267 control points maintain a concordance between verification and classification, which allowed producer and user accuracy to be higher than 98%, and the errors of omission and commission were not higher than 2% (Table 1). For this reason, the pre- and post-processing images were considered reliable.

3.3. Multitemporal dimension

Changes in plant cover due to changes in land use affect the properties and functioning of ecosystems [41,42]. The vegetation and natural covers in the Colombian Andean relief have been modified as a result of human colonization processes. For 35 years in Bijagual, extensive and small areas of forest and paramo were substituted for agricultural and livestock activities, which currently affect habitat loss for biodiversity [23,43]. The landscape structure [44], species dynamics [5], and the ecological network conformation of the natural space of the Tota-Bijagual-Mamapacha complex [22], currently threatened by human activities, have caused [45] a decrease in forest area and an increase in the natural distance of patches, and, therefore, the loss of this ecological continuity corridor [45,46,47] (Figure 1; Table 2).

Forest fragmentation has been more evident since 2021 due to the increase of Dgnm and Hof in the number of patches, but also due to area reduction and patch isolation. For 35 years, 1,923.4 and 614 ha of Dgnm and Hof, respectively, were converted to Mpc (Table 3), supported by the increase in the number of patches, the area, and the ones that become connected (Figure 1 and Figure 2).

Areas around the massif, which have been proposed as transition zones, were found to be affected by crops and livestock activities long before the 1980s. According to Rodríguez-Eraso et al. [48], the loss of natural area in the 80s and 90s was centered in the eastern Cordillera in Boyacá and Cundinamarca, caused by deforestation patterns of the high Andean forest and paramo due to food demand. In addition, Sánchez-Cuervo et al. [24], in the departments of Santander and Boyacá, pointed out that potato cultivation increased from 380 to 482 km² between 2006 and 2008, reducing forest covers.

Potato cultivation has a high per capita consumption in Colombia, and, therefore, expansion has occurred under native forest areas, increasing the continuous degradation of Andean forests and the paramos. Furthermore, another factor influencing forest cover is climate change, which has been enhancing during the last few years and has affected species distribution [49,50]. Since the Andean high mountain ecosystems are seasonal, the increase in temperature would prolong the intensity of the dry periods, and the paramo would reduce the capacity to capture and regulate water [36]. In the Colombian eastern mountain range, only 36.4% of the total paramo area (322,925 ha) and 25.3% of the Andean montane forest (427,675 ha that include the sub-Andean and Andean forest belts) are protected. However, of these areas, only 41 and 45% maintain their original extension [46]. The remaining areas have experienced land use changes with subsequent biodiversity loss [36].

The paramo cover (Dgnm) decreased in area between 1986 and 2021 (1,489 ha), while the forest expanded its area between 1986 and 2000 (309.2 ha) and decreased between 2000 and 2021 (548.7 ha). The forest area that decreased in the second period is larger than the area that expanded in the first (with a difference of 239.5 ha) (Table 2). Instead, the pasture and crop cover increased by 1,722.2 ha in 35 years. Therefore, it is necessary to conserve the Dgnm and Hof areas that still exist in the municipalities of Viracachá, Cienaga, and Tibana (4,034.9 ha [15,18,47]) and include structure, functionality, and ecological diversity attributes in the Bijagual extension [6,35,39], particularly in the municipality of Ramiriquí, which is the most affected by the expansion of the agricultural and livestock frontier (Figure 1).

Spectral image response outcomes show that natural succession has occurred around the paramo and forest areas (Figure 1; Table 3) [23,51,52]. Nevertheless, this process seems to be slow. Lequerica et al. [53] and Prado-Castillo et al. [35] suggest that succession is slow in areas previously used for cultivation due to degradation, and the compaction of soil produced by agriculture can limit the development of propagules and modifies the germinable or viable seed bank [54]. Hence, natural succession needs to be assisted by ecological restoration processes to improve the increase of paramo vegetation.

The increase in the total number of patches (TNP 76 to 396 between 1986 and 2021), density (PD), and percentage of patch occupation (LPI) are evidence of the intensive and extensive land use and the changes from one coverage to another. While in 1986, one patch was registered for each hectare, in 2021, four or five patches (PD) were recorded for each hectare. In 2021, Mpc coverage presented the lowest number of patches compared to Dgnm and Hof (NP 74, 168, and 150, respectively; Figure 2). However, Mpc presents patches with a larger occupation area compared to Dgnm and Hof (LPI 48.7, 6.7, and 6.5%, respectively; Figure 2). In 35 years, the LPI metric describes greater spatiotemporal occupancy of the Mpc patches and lower occupancy of the Dgnm and Hof patches (Figure 2).

The cover and topological position of the Bijagual landscape per patch allow the identification of different contributions to ecological connectivity. Each patch can act as a space in which there is connectivity, i.e., a larger area with a suitable habitat condition or quality that may have higher connectivity [55,56]. The functionality of some patches is evident, even more over those that show higher habitat availability (Figure 3 B,C) since these can act as refuge, feeding and reproduction areas or are habitat units that produce or receive dispersal flows of individuals to other units or patches [56]. In fact, these areas can be analyzed in the metapopulation dynamic to establish possible unknown connections with the biotic complex of Tota to the north and Mamapacha to the south, which, together, comprise a biological corridor.

The spatial configuration of the Hof and Dgnm patches has contributed to the ecological connectivity in Bijagual from 1986 to 2021, evidenced by the contributions of the dPCintra and dPCflux fractions; however, during this period, also natural area loss has occurred, affecting the ecological diversity [29,57,58]. Ecological connectivity in Bijagual decreased noticeably from 2000 to 2021 (Figure 3A). The massif may experience a fragmentation process in the near future [15,45,46], and the current landscape of this natural area shows a loss in habitat quality and availability [45,55,56].

As a result of the distancing and isolation of the Dgnm and Hof patches, the continuity of the ecological corridor has been interrupted (Figure 1). This, in turn, poses difficulties for the movement of species and genetic exchange [56,59] and is evident in the reduced contribution of the dPCconnector fraction (Figure 3A) that determines limited dispersion fluxes. However, the Bijagual Massif still has paramo and Andean forest patches that can act as critical habitat areas for maintaining ecological connectivity (Figure 3C), where habitat loss or deterioration would have an even more negative impact on landscape connectivity. These areas must be prioritized and considered reference ecosystems for ecological restoration processes.

This research states that the Bijagual Massif requires the inclusion of sustainable management practices and territorial planning with agroecological, forestry, environmental education, ecotourism, biogeographical, and climate research proposals, according to the natural environment considering the strategic and vulnerable ecosystem context where the State interacts through its territorial entities (government and environment), the peasant communities (through citizen participation actions of the associated communities), agricultural companies, and the academy as knowledge generators.

4. Materials and Methods

4.1. Study area

The Bijagual Massif is located in the Andes eastern mountain range over the Department of Boyacá (Colombia) between 5°26′0.0″-5°16′0.0″ N and 73°14′0.0″-73°22′0.0″ W [29]. The vegetation is composed of paramo (between 2,990 and 3,460 m a.s.l.) and high Andean forest (between 2,682 and 3,268 m a.s.l.) in an extension of 8,574.1 ha (Figure 4). Its borders delimit to the north with Lake Tota and the paramos of Toquilla, Sarna, Suse, and Los Curíes and to the south with the Mamapacha Massif. Together, they shape the Tota-Bijagual-Mampacha paramo complex [19]. Precipitation values range between 1,000 and 2,500 mm/year, and temperatures range between 13 and 17°C [32].

4.2. Data source

Satellite data interpreted with landscape metrics (i.e., diversity, composition, and configuration) were used to analyze the spatial and temporal dynamics of the Bijagual Massif vegetation. The vegetation cover thematic maps were edited according to the method proposed by The Nature Conservancy [60] and Posada [61] at a scale of 1:100,000. Three images from the Landsat 5 Thematic Mapper (TM; year 1986), Landsat 7 Enhanced Thematic Mapper (ETM+; year 2000), and Landsat 8 Operational Land Imager (OLI; year 2021) sensors with a spatial resolution of 30 m, Path/ Row 7/56, and cloudiness < 10% (source: Earth Explorer of the USGS - U.S. Geological Survey) were employed. Band combinations in each image were made according to the specifications of each sensor to differentiate covers. For the 1986 and 2000 images, the combination was Red: 4, Green: 5, and Blue: 3, and for 2021, Red: 5, Green: 6, and Blue: 2. Pre- and post-processing was performed with the programs ERDAS Image ver. 2016 and ArcGIS ver. 10.7 (academic licenses) described below.

4.3. Geometric correction

One hundred and fifty control points (in images 1986 and 2021) were generated from a third-order polynomial model and a georeferenced image (image from 2000; [17,40]); in the 1986 image, 21 automatic points and 129 manual points were projected (standard error-RMSE 0.59; Std. Dev: 0.218); in 2021, 10 automatic and 128 manual points (RMSE 0.55; Erros Std. Dev: 0.208) were projected. These were considered reliable because their values were < 1.0 [40]. Subsequently, the study area was cut from the Bijagual polygon.

4.4. General contrasts (spatial enhancement)

The images were smoothed in their spatial contrasts similar to the digital Level of each pixel with respect to neighboring pixels using low-pass filters (Kernel 3X3) to reduce spatial variability and visualize less sharp, more blurred, and homogeneous profiles.

4.5. Land covers (thematic categorization) – field verification

A statistical digital classification was performed for each image to group the pixels by cover, i.e., going from a continuous measurement (digital levels of each pixel) scale to a categorical one (vegetation coverage in the Corine Land Cover system – IDEAM [28]) using the supervised classification method and the selection and delimitation of training areas (grouping of pixels of the same category).

The supervised classification method was supported with field verification activities, which involved obtaining geographic coordinates (control points in the WGS84 system) and floristic and vegetation characterizations in each natural cover type (published in Gil-Leguizamón et al. [29]). In addition, for endemic and threatened species, Herbarium collections and online catalogs were consulted [30,32].

4.6. Classification accuracy

The reliability of the training areas (spectral separability of the categories) was statistically evaluated through an error matrix that contrasted the real classes (the truth of the terrain) and the cartographic units (cover). Additionally, accuracy measurements were used for each cover to determine the reliability of the supervised classification, including:

4.6.1. Producer accuracy

The omission error measure indicates the probability that a sample point is correctly classified. It is calculated by dividing the number of pixels correctly classified into a category by the total number of pixels in that category according to Equation (1).

$$\mathit{P}\mathit{P}\mathbf{\%}=(\mathit{X}\mathit{i}\mathit{i}/\mathit{X}+\mathit{i})\mathbf{\ast }100$$

(1)

Where PP% is producer accuracy in %, Xii is the value of the column diagonal, and X+i is the total marginals of column i.

4.6.2. User accuracy

The probability that a sampling point corresponds to the category assigned in the field (commission error). It is calculated by dividing the total number of correct pixels in a category by the total number of pixels that, in fact, belong to this category according to Equation (2).

$$\mathit{P}\mathit{U}\mathit{\%}=(\mathit{X}\mathit{i}\mathit{i}/\mathit{X}\mathit{i}+)\mathbf{\ast }100$$

(2)

Where PU% is user accuracy in %, Xii is the diagonal value of row I, and Xi+ is the total marginals of this row.

4.6.3. Kappa coefficient

It measures the difference between what is observed between classes and what is evidenced in the classification according to Equation 3. The veracity corresponds to values close to or equal to 1.

$$\mathit{K}=(\mathit{O}-\mathit{E})/(1-\mathit{E})$$

(3)

Where O is the correct proportion of the classes identified in the image, and E is the correct observed percent change estimate.

4.7. Multitemporal dimension (landscape metrics)

Statistics to identify the spatial and temporal dynamics of the coverages for 1986, 2000, and 2021 were obtained according to the spatial pattern theory [62]. This incorporates the landscape as a mosaic of patches or discrete covers, which explain the distribution of geographic objects, their patterns, and processes in time [62,63]. The statistics used quantified the landscape with three parameters: diversity and composition with Fragstats [62] and connectivity with Conefor [55,56,64], as follows:

4.7.1. Diversity

It was analyzed through the total landscape area (TA, ha), the total number of patches (TNP), patch richness (PR), and patch richness density (PRD).

4.7.2. Composition

Includes the area of each cover or class (CA, ha), the number of patches (NP), patch density (PD), landscape percentage (PLAND, %), and the largest patch index (LPI, %) to analyze the composition.

4.7.3. Connectivity

The topo-ecological index of connectivity probability (dPC_k) was used and interpreted according to three fractions that quantified how the patches (as habitat or corridor) contribute to the global connectivity of the Bijagual landscape, according to Equation 4, to obtain connectivity information.

$${\mathit{d}\mathit{P}\mathit{C}}_{\mathit{k}}={\mathit{d}\mathit{P}\mathit{C}\mathit{i}\mathit{n}\mathit{t}\mathit{r}\mathit{a}}_{\mathit{k}}+{\mathit{d}\mathit{P}\mathit{C}\mathit{f}\mathit{l}\mathit{u}\mathit{x}}_{\mathit{k}}+{\mathit{d}\mathit{P}\mathit{C}\mathit{c}\mathit{o}\mathit{n}\mathit{n}\mathit{e}\mathit{c}\mathit{t}\mathit{o}\mathit{r}}_{\mathit{k}}$$

(4)

Where k is a landscape element (patch), intra is the contribution of patch k in terms of the habitat area (or quality) available within it (intrapatch connectivity), and flux corresponds to the dispersion flux (weighted by the area) received or originated through the connections of patch k with the rest of the patches (it reflects how well-connected k is with the rest of the existing habitat in the landscape). Connector evaluates the contribution of k as a connecting element or bridge patch between the rest of the patches (estimates to what extent k facilitates dispersive flows and which do not have their origin or destination in k but are enhanced and pass through k) [56,65]. Furthermore, an intermediate dispersal distance (0.5 probability) and an approximate average distance between patches of 50 m were assumed. The results are shown as bar and scatter plots obtained with the ggplot2 package [66] included in the open-source RStudio Desktop version software of 2003 for Windows 10/11.