Secondary Forest and Shrubland Dynamics in a Highly Transformed Landscape in the Northern Andes of Colombia ( 1985 – 2015 )

Understanding the dynamics of natural ecosystems in highly transformed landscapes is key to the design of regional development plans that are more sustainable and otherwise enhance conservation initiatives. We analyzed secondary forest and shrubland dynamics over 30 years (1985–2015) in a densely populated area of the Colombian Andes using satellite and biophysical data. We performed a land-cover change analysis, assessed landscape fragmentation, and applied regression models to evaluate the effects of environmental and geographical correlates with the observed forest transitions. Forest cover area increased during the 30 year-span, due mostly to forest regrowth in areas marginal for agriculture, especially during the first half of the study period. However, a high dynamic of both forest regrowth and clearing near urban centers and roads was observed. Soil fertility turned out to be a key correlate of both forest recovery and deforestation. Secondary forests, <30 years old represent the most fragmented component. Our findings reflect the complexity of the processes occurring in highly transformed and densely populated regions. Overall, this study provides elements for a better understanding of the factors driving land cover change near large urban areas, and raises new iideas for further research.


Introduction
As a result of the continuous growth of the world's population [1], the natural environment is suffering from massive habitat losses [2,3].Land transformation is one of the main drivers of ecosystem degradation [4,5], and biodiversity loss [6,7].Land-cover change also affects water and energy budgets, thereby altering climates [8,9].Assessing and monitoring land-cover dynamics is therefore key to understanding and managing changes in the state, structure, and function of terrestrial ecosystems [3,10].
In tropical South America, one of the regions harboring the greatest concentration of biodiversity worldwide [11], export-driven agriculture [12], subsistence cropping and road development [13,14] are the main drivers of deforestation and the major cause of land cover change [15,16].Lowland tropical forests, particularly in the Amazon Basin, have received much attention because they contain a large proportion of the world's biodiversity and have experienced unprecedented rates of deforestation [17].Yet many other tropical ecosystems such as the Brazilian Cerrado, Atlantic forests, Chilean matorrals, and tropical dry forests, have lost proportionately more natural vegetation cover than the forests of the Amazon Basin.The deforestation of these ecosystems adds up to 3.6 million km 2 [18].Indeed, between 2000 and 2012, 58% of the deforested area in South America took place outside the Amazonian region [19].
Among the most diverse and threatened tropical ecosystems are those hosted by the Andean mountain ranges, considered a major biodiversity hotspot [20].The Northern Andes in Colombia are particularly vulnerable to land cover change, as they own an extraordinary level of species diversity and endemism [11,21].Unfortunately, Andean forests have been subject to substantial forest loss, as they have been densely settled and exploited since the Spanish colonization [22].Current drivers of landscape change are cropping and cattle grazing, and to a lesser extent, to mining.In remote areas, illicit crops are also an important cause of deforestation and forest fragmentation [23][24][25].As a result, only 31% of the Andean natural forests were left at the beginning of the century [26,27].
During the last decade, Andean ecosystems have experienced considerable vegetation recovery thanks to land abandonment resulting from spontaneous or forced population migration.As an example, in the 2001-2010 period, about 65% of the total national woody net gain occurred in different Andean montane forest ecoregions in South America [28].Secondary forests provide multiple ecosystem services and serve as critical reservoirs for thousands of terrestrial species [3,29], attenuating the footprint left by anthropogenic activities.Forest regrowth and deforestation are therefore opposing forces acting on the same landscape simultaneously.In this context, one may ask how forest cover has changed in the last few decades.
Among the Andean regions that have undergone the most substantial land cover change is the high-plain around Bogotá, Colombia's capital, located in the eastern Andean mountain range.Bogotá has experienced fast economic growth over the last 60 years.The city of Bogotá reached a population of over 8 million people in 2015 from about 0.9 million in 1957 [30].Today, the metropolitan area surrounding Bogotá includes at least 15 municipalities, and is home to about 25% of the country's entire population [31], making it the second largest city in South America, after Saõ Paulo in Brazil.This substantial population growth has increased pressure on the surrounding landscape, as well as on the provision of natural resources and ecosystem services.The dynamics of the complex landscape mosaic near Bogotá (e.g., deforestation, land degradation, and local forest regrowth) are not well known.Assessing the extent of forest loss and recovery over the last few decades, and its relation with the different drivers of landscape transformation is critical to the development of effective management strategies in a region substantially transformed by human activities, but that host dozens of endemic species [32].
Following these considerations, the main objective of this study was to assess the temporal and spatial dynamics of forests and shrublands in a highly human-impacted Andean ecosystem.To this end, remote sensing, geographical information systems and statistical techniques were used to evaluate temporal changes in land-cover and landscape structure of the last three decades in the northern Bogotá high-plain.These changes were then related to environmental and landscape variables.More precisely, the specific objectives of the study were: (i) To map forest and shrub land using multi-temporal multispectral medium resolution satellite data covering 30 years.(ii) To quantify land cover change and fragmentation patterns for different types of vegetation.(iii) To identify the major environmental correlates of land-cover changes over time.

Study Area
The study area is located in a high-plain region north of Bogotá city, within the Cundinamarca Department.It encompasses about 136,000 ha and comprises a mountain plateau located in the eastern Andean cordillera, at an average altitude of 2600 m above sea level (a.s.l.).(Figure 1).The high-plain is crisscrossed by hills and mountain ranges that reach over 3500 m a.s.l.The average temperature is about 14 • C, while annual average precipitation varies from 600 mm in the center of the study area, increasing to 800-1200 mm towards the western and eastern borders respectively, with rainy seasons in April-June and September-November [27].The region is the most densely populated in Colombia, with a population of nearly 9 million [30].It is an important agro-industrial center of the country, with flower growing industries and dairy farming, livestock and mining being the predominant land uses [33].The study area borders with the northern section of the capital city of Bogotá, with a population of 8 million people, and it is part of the socio-economic influence region of the city [33].
Forests 2017, 8, 216 3 of 17 high-plain is crisscrossed by hills and mountain ranges that reach over 3500 m a.s.l.The average temperature is about 14 °C, while annual average precipitation varies from 600 mm in the center of the study area, increasing to 800-1200 mm towards the western and eastern borders respectively, with rainy seasons in April-June and September-November [27].The region is the most densely populated in Colombia, with a population of nearly 9 million [30].It is an important agro-industrial center of the country, with flower growing industries and dairy farming, livestock and mining being the predominant land uses [33].The study area borders with the northern section of the capital city of Bogotá, with a population of 8 million people, and it is part of the socio-economic influence region of the city [33].The high-plain hills/mountains are composed of Cretaceous and Tertiary shales and sandstone, whereas the plain is composed by Quaternary sediments of fluvial and lacustrine origin [34].The native vegetation of the area comprises a variety of vegetation types and species depending on soils and climatic conditions.Original forests of the plains were dominated by species such as Quercus humboldtii, Weinmania tomentosa and Vallea stipularis [35], and mature forests hosted many species The high-plain hills/mountains are composed of Cretaceous and Tertiary shales and sandstone, whereas the plain is composed by Quaternary sediments of fluvial and lacustrine origin [34].The native vegetation of the area comprises a variety of vegetation types and species depending on soils and climatic conditions.Original forests of the plains were dominated by species such as Quercus humboldtii, Weinmania tomentosa and Vallea stipularis [35], and mature forests hosted many species belonging to the Lauraceae family, but only small patches of intact forest remain due to logging.Nowadays, the most common and species-rich families are Melastomataceae and Ericaceae.In drier areas, one may encounter shrubs of the Morella and Myrcianthes genus.The old-growth vegetation of the hills (2700-3300 m a.s.l.) is typically dominated by trees and shrubs of the genus Weinmania, Clusia, Cedrela and Miconia [36].From 3300 m to 3500-3600 m a.s.l., the sub-páramo prevails, a transition between forest and tropical highland grasslands dominated by shrubs and grasses belonging to the Asteraceae family (genera Espeletia and Espeletiopsis) [37].

Vegetation Cover Mapping
Three Landsat TM and one Landsat OLI images were used to map recent changes in land cover.They were downloaded from the USGS Glovis server [38], spanning 30 years of landscape transformation, and corresponding to 1985, 1991, 2001 and 2015 (Table 1).The uneven length of the time intervals is due to both the difficulty of finding images with a low proportion of clouds/shadows and to errors present in the Landsat 7 ETM+ SLC-off imagery.Standard pre-processing of the images involved radiometric calibration [39] and atmospheric correction, using the Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) model [40,41].Subsequently, a parallelepiped supervised classification [42] to the subset of images corresponding to the study area for each year were used.Two main vegetation classes of interest were defined: (i) forest, which includes mature primary, secondary and planted forest; and (ii) shrub land, referring to natural shrubs and early vegetation regrowth.All the land that did not enter into these two vegetated classes of interest was merged into the agricultural and built-up lands class (ABL), which includes productive land, such as agriculture and pastureland, and built-up areas.The mapping task focused mainly on the forest and the shrubland vegetation classes.All other classes (crops, pastures and built up lands) were thus grouped into a single "agricultural and built up lands" class.In order to train the classifier, a number of training samples were defined from historical aerial imagery available in the national Colombian Geographical Institute-IGAC-for the year 1985 [43], and Google Earth historical imagery.We selected, under a simple random sampling scheme, 219, 128, 235 and 200 training points for the forest class, and 222, 125, 100 and 200 points for the shrub land vegetation class, for years 1985, 1991, 2001 and 2015, respectively.Changes in land cover areas were calculated (i) using a standard post-classification approach, and (ii) integrated by uncertainty measures (confidence intervals) of error adjusted area estimates using the methods described by Olofsson et al. [44,45] (Table S3 in Supplementary Materials).
A single mask with no land cover information was also produced by merging the classified clouds and shadows for all four years.Finally, an accuracy assessment was performed using the same sources used for training the classifier.Standard accuracy metrics were calculated using 346 (1985), 394 (2001) and 293 (2015) randomly distributed validation points.No validation points were available for 1991.

Landscape Change Analysis
To avoid underestimation of the land-cover pattern quantification due to the rugged topography of the study area, a patch area base correction was performed [46,47].The DEM Surface Tools ArcGIS extension [48] and an ASTER Global Digital Elevation Model version 2 [49] were used to compute the surface area of each pixel based on the triangulation method proposed by Jennes [50].Surface area for every patch and land-cover class was obtained by summing the values of its pixels [46,47].
Using the four land cover maps, transition matrices were calculated for contiguous periods of time, and for the overall period , based on a pixel-by-pixel comparison [51].
To identify the spatial distribution of different vegetation successional stages, we derived a 'forest age' map following Etter et al. [52], by assessing the presence or absence of each forest class in the classified imagery for the four dates (see Table S1 in Supplementary Material) and by establishing different forest age classes as follows (following the uneven distribution of the available images): -A4: ≥30 years of presence.-A3: 24-30 years.-A2: 14-24 years.-A1: ≤14 years.
Additionally, to describe the temporal landscape structure variation, a landscape pattern analysis was performed by computing a set of landscape metrics for the four land cover maps and the Forest Age map [53].The analysis was carried out for each vegetation class (forest and shrub land) separately, using the Patch Analyst extension in ArcMap 10.1 (Environmental Systems Research Institute, Redlands, CA, USA).The landscape metrics selected were: number of patches (NP), mean patch size (MPS), total core area (TCA), mean core area per patch (MCAP) and core area percentage of landscape (CORLAND) [54].NP and MPS are measures of landscape fragmentation and characterization [55], and are especially meaningful when comparing the same landscape in different time periods.TCA refers to the total amount of core area or area not affected by edge effects, MCAP refers to the distribution of the core area within the class patches, and CORLAND a measure of the dominance of the core areas in the landscape.In general, core area metrics are strongly related with the edge effect and therefore are indicators of habitat quality for many species highly dependent on core area availability [55,56].

Geographical Factors and Change
To gain an understanding of the factors that may explain the spatial dynamics in vegetation change, a multinomial linear regression model was used to evaluate the effect of environmental (soil fertility, slope) and geographical (distance to roads, distance to urban) explanatory variables on forest transitions for the 1985-2015 period [53,57].A soil fertility map was created by establishing five fertility classes (very low, low, intermediate, high, very high) based on the national soil map [58].The slope map was derived using the ASTER Global Digital Elevation Model version 2 [49].Euclidean distances to roads and municipalities were calculated based on national datasets from the Planning and Land Management Geographical Information System of Colombia SIGOT [59].All geographical layers were fitted to the 30 m raster format of the land cover maps.
Multinomial regression models generalize logistic regression to multiclass variables; that is, more than two possible discrete outcomes in the dependent variable.In other words, the transition probabilities of the different possible outcomes for each land cover class (forest, shrubland and agricultural and built-up lands) were predicted.For each land cover class, there are three possible discrete outcomes (e.g., for forests: forest to forest, forest to shrub and forest to ABL).Thus, three separate multinomial regression models for each land cover class were run to estimate the transition probabilities depending upon slope, soil fertility, distance to urban centers and distance to roads (predictive variables).In the first model, the possible outcomes of the dependent variable were forest to forest, forest to shrub and forest to ABL (n = 157,522).In the second model, the possible outcomes of the dependent variable were shrub to shrub, shrub to forest and shrub to ABL (n = 233,971).Finally, in the third model the possible outcomes of the dependent variable were agricultural and built-up lands (ABL) to ABL, ABL to forest, and ABL to shrub (n = 1,074,224).
In a second step, another multinomial linear model was built with the Forest age map, to predict the probability of vegetation cover being shrub land, young secondary forest (by merging classes A1 and A2), or late secondary and mature forest (by merging forest age classes A3 and A4) against the same predictive variables (n = 1,465,717).To deal with spatial autocorrelation, each of the models was run using a randomized subsample of 10% of the entire data set.This procedure was repeated 1000 times and the mean of each estimate was calculated, as well as the one-sided 95% confidence intervals, depending on the sign of the mean of the parameter.The multinomial models were run using the R statistical package software [60].

Land Cover Change Analysis
Four land cover maps, covering the years of 1985, 1991, 2001 and 2015, were produced (Figure 2).The land cover maps had overall accuracies of 89% (1985), 81.4% (2001) and 83.6% (2015).The 1991 map accuracy was not calculated due to the absence of a validation dataset for that year (Table 2 reports all accuracy measures).The analysis of 30 years of land cover changes in the Bogotá high-plain indicated a general trend of woody vegetation regrowth across the study area.The forest class increased by 42% during the study period, from 16,100 ha in 1985 to about 22,850 ha in 2015.The shrubland class also increased by a similar amount (43.3%) from 23,300 ha in 1985 to 33,400 ha in 2015.This gain in woody vegetation covers was at the expense of other land covers previously used for crops and pasture, which showed a decrease of 16,850 ha between 1985 and 2015 (−16.5%).The intermediate period between 1991 and 2001 was characterized by the strongest trend of forest regrowth, while the period 2001-2015 showed a stabilization in land cover change (Figure 3).The most frequent transition in the whole study period was from shrubland to forest, indicating an important process of forest recovery.In contrast, the shrubland class showed a high degree of change, both toward recovery or toward loss (Table 3).Land cover error-adjusted area estimates with 95% confidence intervals are reported in Table S3 (Supplementary Materials, using methods from [44,45]).Overall land cover changes  were spatially represented by nine different land cover transitions (Figure 4.) Altitudinal changes were also analyzed.The most stable areas of shrub lands and forests were found at altitudes over 2800 m (Figure 5).The altitude belt over 2600 m was the most dynamic, showing a strong increase in forest and shrubland areas.On the other hand, around 50% of the agricultural and built-up lands (ABL) in 1985 below 2600 m remained stable over the study period (Figure 5).

Patterns of Landscape Change
About 45% of forest cover remained stable over the study period (forest age class A4).These older secondary forests are located in the highest elevations and steeper slopes, where anthropogenic disturbance is less likely and where less fertile soils are found.Young secondary forests (forest age classes A1 and A2; Table S2) were more common at lower altitudes, where land cover change is more dynamic.Forest regrowth occurred to a large extent in the western part of the study area (Figure 6).
Land cover change processes had a large influence on the spatial structure of woody vegetated areas, showing an overall increase in the number of forest patches, resulting from a scattered pattern of forest regrowth.Mean patch size (MPS) and mean core area (MCA) oscillated in size during the study period, yet forest patches showed a notable decrease in size in 2015 compared to 1985.The number of forest patches belonging to age class A4 (≥30 years of presence) was lower than that of younger secondary forests (age classes A1 to A3; i.e., ≤30 years of presence), and showed a higher mean patch size and mean core areas.For the shrubland class, the number of patches also increased over time, as well as their mean patch size and core area (Table 4).

Patterns of Landscape Change
About 45% of forest cover remained stable over the study period (forest age class A4).These older secondary forests are located in the highest elevations and steeper slopes, where anthropogenic disturbance is less likely and where less fertile soils are found.Young secondary forests (forest age classes A1 and A2; Table S2) were more common at lower altitudes, where land cover change is more dynamic.Forest regrowth occurred to a large extent in the western part of the study area (Figure 6).
Land cover change processes had a large influence on the spatial structure of woody vegetated areas, showing an overall increase in the number of forest patches, resulting from a scattered pattern of forest regrowth.Mean patch size (MPS) and mean core area (MCA) oscillated in size during the study period, yet forest patches showed a notable decrease in size in 2015 compared to 1985.The number of forest patches belonging to age class A4 (≥30 years of presence) was lower than that of younger secondary forests (age classes A1 to A3; i.e., ≤30 years of presence), and showed a higher mean patch size and mean core areas.For the shrubland class, the number of patches also increased over time, as well as their mean patch size and core area (Table 4).S1 and  S2 for age class definitions and extensions.

Landscape Change Correlates
The regression models for land cover transitions during the study period showed clear opposite trend processes related to land clearing (forest or shrubland to agricultural and built-up lands) and those related to woody vegetation regrowth (agricultural and built-up lands to forest or to shrubland).
As expected, distance to roads and to urban centers showed a strong effect on land cover dynamics.The nearer to urban centers and roads, the more forest areas were likely to become shrubs or agricultural and built-up lands by 2015.Likewise, shrubs were more likely to become forests or ABL lands than to remain stable.In contrast, ABL lands in 1985 were more likely to remain stable when they were near urban centers and roads (Table 5).
Environmental variables also appeared to affect land cover change.Steep slopes favored persistence and gain of woody vegetation cover.The steeper the slope, the less likely were forests to become shrubs and shrubs to become cleared; and agricultural and built-up lands were more likely to become shrubs or forests.Fertile soils showed higher dynamic.Here, forests were more likely to become shrubs or ABL lands, than to remain forests.Likewise, shrublands were more likely to become forests or ABL lands.In contrast, less fertile soils were more likely to show shrublands and forest regrowth (Table 5).
The forest age model showed that all woody vegetation (forest and shrubland classes) was more likely to be present far from urban centers and roads, on steeper slopes, and on fertile soils.These effects were stronger for forest age classes A3 and A4.In addition, shrublands and young secondary forests (<24 years old) were scarcer in fertile soils (Table 6).

Discussion
Although this landscape has been highly transformed for a long time [26], it continues to show strong dynamics driven by an apparent land use dis-intensification and land abandonment in the mountains and hills of the region.Overall, this study shows two important aspects: (i) a general increasing trend in woody cover, and (ii) although land cover changes strongly varied between periods and between altitudinal belts, the period 1991-2001, and above 2800 m, was dominated by woody cover increase.The abandonment of agricultural and built-up lands and shrublands followed a successional process, resulting in increasing areas of young or intermediate secondary forests across the landscape.These findings support previous studies showing that forest cover is increasing in the southwest of the Bogotá high-plain region [27], and is consistent with the trend of woody vegetation recovery reported for Colombia at the national level during the last decade [28].Indeed, net gains in forest cover have also been observed in several other areas of Latin America and the Caribbean (e.g., [61][62][63][64]).These patterns have been associated with rural-to-urban-migration dynamics, which reduce the pressure to convert forest patches, and promote forest regrowth through land abandonment [65,66].By providing increased labor opportunities and access to higher education, urban centers have become very attractive to rural populations.
Despite the fact that 45% of the forested areas remained stable during the study period, the study found considerable dynamics in many forest and shrubland patches.Young secondary forests increased in area, but mostly in a pattern of small and disconnected patches, mainly in lower elevations and accessible areas.Forests that were already present at the beginning of the study period represent only a small fraction of areas covered by woody vegetation in the study area.However, these patches have the largest mean core area and mean patch size, and represent remnants of the forests covering the whole area in the past.This is remarkable as most Andean slopes in this region were strongly deforested a century ago as a result of cropping activities, cattle grazing [26], and because wood was a critical source of cooking energy during the early 20th century [67].If these remnant forests indeed correspond to old-growth forests, they should harbor considerable levels of biodiversity, and could therefore be used as core areas to implement a protected areas network that can integrate the increasing areas of regrowth forests to shield against the high conversion risks [28,68].
During the study period, the number of patches of forest vegetation increased, but both mean patch area and mean core area decreased.These seemingly contradictory trends may be the result of two processes occurring simultaneously.First, the proliferation of patches indicates that vegetation regrowth is taking over after land abandonment which may be associated with soil productivity loss, economic marginalization and/or agriculture intensification [4,65].Second, the decline in mean patch and core area probably results from the process of fragmentation associated with land conversion, which is still underway, with some smaller areas converted to productive lands.Size reduction in forest patches may also be the result of activities related with residential construction and industrial development [69], often occurring in peri-urban areas [70].Future studies evaluating vegetation stand dynamics at the local scale can help clarify this issue.
Regression analyses showed that landscape variables had an important effect on forest transitions.First, the study found that there is an important dynamic towards both forest regrowth and clearing near urban centers and roads.Other theoretical and empirical studies have shown that deforestation processes are more likely to occur in accessible areas, while forest regrowth typically takes place in less accessible ones [71][72][73][74][75][76].Patterns of recovery near urban areas and roads may be explained by the increase of low-density suburban settlements of secondary homes of people living in Bogotá, well represented in the area, who are promoting forest regeneration within their properties.Another possibility is that access to roads is positively related to land cover change dynamics.In other words, transitions to forests or to shrublands near roads or urban centers could be the result of previously cleared lands where vegetation is now restoring, but that is likely to be become agricultural or built-up land in the near future.In 2015, however, woody vegetation was more likely to be found far from urban centers and roads.Together, these findings suggest that, although areas close to urban centers may be associated to some extent of forest regrowth, the larger fraction of forest cover, especially of older forest age classes, is found far from human settlements.
Along with human settlements, fertile soils also turned out to be a key driver of both forest recovery and degradation simultaneously.Although mature remnant forests were found in more fertile soils, which is somewhat counterintuitive, they are associated also with steeper slopes and are far from roads and settlements.In turn, regrowth forests occur on less fertile soils, which are more likely to be abandoned.Regrowth also occurs on steeper slopes and further away from roads and settlements.Indeed, soil fertility is widely recognized as a driver of either deforestation or forest succession [77,78], as it acts as a land use regulator by determining the agricultural value and economic potential of an area.For instance, forests occurring on fertile soils are more likely to undergo conversion towards farming, while converted lands on less fertile soils are more likely to be abandoned, allowing succession to result in new secondary forests [65].Also, some forest patches may have remained intact for several decades as they are often located at the top edge of mountains and along steep cliffs.Another possibility is that of an increased consciousness of the local populations about the biodiversity value and the ecosystem services provided by forest cover [32].In marginal areas with steep slopes, forests were more stable, and agricultural and built-up lands were more likely to become shrubs or forests.This is in agreement with Mendoza and Etter [27], who found that more than 70% of the remnant forests were located in areas with slopes over 25%.
Finally, these results raise new research questions that need to be further investigated at both local and regional scales.In particular, it is critical to evaluate forest resilience.Although the ecological value of secondary forests has grown in the last decade as a result of a global decrease of primary forests [79], floristic reassembly is a very slow process, and the original species composition may never return to its original state [80,81].For instance, many secondary forests in the region are dominated by Ulex europaeus, an invasive species prone to fire, which is very likely to replace native species, and is difficult to eradicate [82].Also, shrublands are an important component of the current landscape, which raises the question whether these correspond to vegetation which is gradually recovering from strong past disturbances, or whether they represent an arrested successional stage owing to soil degradation and/or seed limitation [83].In this context, future research should involve more detailed studies in the Bogota high-plain, comparing the successional dynamics among shrubland and early and late secondary forests in different fragmentation scenarios.Andean forests hold key ecosystems with the potential of being important biodiversity reservoirs and of providing important ecosystem services related to water regulation and carbon sequestration [84].Historically, much attention has been given to tropical lowland rainforests.Here, the importance to broaden ecological research to these less studied tropical ecosystems is stressed, where our knowledge about carbon and water budgets is still limited [84].
Some limitations of the research derive from the absence of validation points for the 1991 land cover map.Additionally, measures of error adjusted areas and relative confidence intervals should have been calculated also for the land cover change transition matrix (1985-2015) [44,45], which was not possible due to the absence of an ad-hoc change validation dataset.Finally, the study covered only part of the northern region neighboring Bogotá; further analysis to corroborate the generality of our conclusions should involve research in additional areas of the Bogotá high-plain.

Conclusions
Evidence of an important forest recovery process in a surrounding area of the city of Bogotá, a large urban region of the Eastern Colombian Andes, is provided.Although we may be optimistic about the future of forest cover recovery in the Andean region, our findings also reflect the complexity of the processes occurring in such a highly transformed and densely populated region, where the same factor can be a driver of both deforestation and recovery simultaneously.Overall, our results bring noteworthy elements to the understanding of the factors driving land cover change in large peri-urban areas, which greatly benefit from the environmental services that surrounding forest cover provides.A deeper understanding of the drivers of land cover change in these ecosystems is crucial for the design of regional environmental plans that can direct more sustainable settlement models, enhance conservation initiatives and improve the delivery of ecosystem services.

Figure 1 .
Figure 1.Location and elevation map of the study area; major roads (red lines), water bodies (blue), and villages (black dots).

Figure 1 .
Figure 1.Location and elevation map of the study area; major roads (red lines), water bodies (blue), and villages (black dots).

Figure 3 .
Figure 3. Changes in the area covered by each land-cover class between 1985 and 2015 (ha).

Figure 3 .
Figure 3. Changes in the area covered by each land-cover class between 1985 and 2015 (ha).

Figure 3 .
Figure 3. Changes in the area covered by each land-cover class between 1985 and 2015 (ha).

Figure 4 .
Figure 4. Land cover transitions for the 1985-2015 period.Forest to Agricultural and built up lands (ABL) (F-A); Shrubland to ABL (S-A); Stable Agricultural and built up lands (A-A); Forest to shrub land (F-S); Stable shrubland (S-S); ABL to shrubland (A-S); Stable forest (F-F); Shrubland to Forest (S-F); ABL to Forest (A-F).

Figure 4 .
Figure 4. Land cover transitions for the 1985-2015 period.Forest to Agricultural and built up lands (ABL) (F-A); Shrubland to ABL (S-A); Stable Agricultural and built up lands (A-A); Forest to shrub land (F-S); Stable shrubland (S-S); ABL to shrubland (A-S); Stable forest (F-F); Shrubland to Forest (S-F); ABL to Forest (A-F).

Figure 6 .
Figure 6.Forest age map of the study area (1985-2015).A1 to A4: forest age classes.See TablesS1 and S2for age class definitions and extensions.

Table 1 .
Landsat images used in the present study.

Table 2 .
Land-cover maps accuracy assessment.

Table 4 .
Changes in landscape pattern indicators during the study period in forests and shrubland, and for each forest age class.

Table 4 .
Changes in landscape pattern indicators during the study period in forests and shrubland, and for each forest age class.
Figure 6.Forest age map of the study area (1985-2015).A1 to A4: forest age classes.See Table S1 and S2 for age class definitions and extensions.

Table 5 .
Parameter estimates of the model predicting forest transitions between 1985 and 2015 from environmental and landscape variables.The 95% confidence intervals (CI) obtained from 1000 randomized subsamples of the entire dataset are in parentheses.The parameter estimate is significant if the 95% CI does not overlap with zero.

Table 6 .
Parameter estimates of the model predicting forest age.The 95% confidence intervals (CI) obtained from 1000 randomized subsamples of the entire dataset are in parenthesis.The parameter estimate is significant if the 95% CI does not overlap with zero.