Introduction

Protected areas are designated to conserve forest ecosystems and their associated biodiversity. These areas not only preserve species and biodiversity but also ensure the regeneration of natural resources that depend on forested environments. Over time, some of these areas develop into enduring landscapes that contribute to ecological stability, while others require adaptive management to address ongoing changes and threats-such as deforestation or degradation-thus ensuring the long-term preservation of forests and their ecological functionality.

Anthropogenic intervention is one of the main causes of the decline in native flora and fauna species, as well as forest fragmentation resulting from road construction, which disrupts habitat connectivity and reduces the functionality of ecological niches (Comisión Nacional de Áreas Naturales Protegidas [CONANP], 2023).

Monitoring protected areas plays a critical role in supporting the design and implementation of public policies in the forestry sector, as well as in decision-making processes related to conservation, protection, restoration, and the promotion of sustainable forest ecosystem management (Comisión Nacional Forestal [CONAFOR], 2023). In addition, documenting historical changes in land use and vegetation patterns is essential, as it provides a crucial temporal context for ecological studies and informs conservation policy decisions (Duarte et al., 2016; Gordillo-Ruiz & Castillo-Santiago, 2016; Hernández-Cavazos et al., 2023).

Remote sensing technologies enable the accurate quantification and spatial mapping of forest distribution by producing highly reliable vegetation indices for estimating forest cover (Lumbres & Lee, 2014). To enhance the consistency and comparability of change analyses over time, it is essential to standardize land cover classification systems. Such standardization supports the development of consistent maps across different time periods, thereby improving the precision and reliability of assessments related to landscape transformation processes (Duarte et al., 2016).

Assessing land cover changes and current conditions over a specific time period allows for the identification of losses or gains in natural cover and provides key insights into the state of conservation. Land cover and land use change analysis is a widely used tool in studies of deforestation, environmental degradation, biodiversity loss, and alterations to the hydrological cycle (Pontius & Petrova, 2010; Wang et al., 2021). The significance of this type of analysis lies in its ability to quantify transformations and project future scenarios based on historical records (Sánchez-Reyes et al., 2017). To ensure accurate interpretation, it is important to distinguish between land cover and land use. These concepts together offer a more precise understanding of the spatial configuration of a given territory. Land cover refers to the biophysical elements present on the Earth’s surface, such as vegetation, water bodies, or urban areas, whereas land use relates to the human activities that alter or maintain that cover for productive, residential, or recreational purposes (Hernández-Cavazos et al., 2023; Pineda Jaimes et al., 2009; Sandoval-García et al., 2021).

The Nevado de Toluca Flora and Fauna Protection Area (APFFNT) is part of the conservation area system of Estado de México. The region is surrounded by communities undergoing continuous urban expansion, land use changes, and different forms of forest management (CONAFOR, 2023). These environmental and social dynamics have influenced conservation strategies and are particularly reflected in the change to the area’s management category -from its designation as a national park in 1936 to its reclassification as a flora and fauna protection area in 2014 (Olvera & Pichardo, 2017). Therefore, analyzing this area is essential to determine whether the protection and conservation goals have been met within its boundaries. In this context, the objective of the present study was to analyze the transitions and dynamics of forest cover and land use in the APFFNT for the period 1971-2021.

Materials and methods

The study area is located in Estado de México and is one of the 97 Natural Protected Areas (NPAs) of the entity, covering a total area of 53 590.67 ha. Currently, the area is managed by the Comisión Estatal de Parques Naturales y de la Fauna (CEPANAF, 2022).

Cartographic information on land use and vegetation was obtained through the photointerpretation of aerial images, enabling the classification of forest cover and land use based on the land use maps from the INEGI (Instituto Nacional de Estadística y Geografía, 2023) and the classification system proposed by Miranda and Hernández (1963).

Digital files and the final analysis were generated using geographic information systems, following the methodological framework outlined in Figure 1.

The process of gathering and compiling information is often complicated, notably when it requires consistency with existing maps. A critical issue when comparing multiple datasets or studies is the spatial resolution or scale of cartographic materials, which is fundamental for temporal comparison. Occasionally, it is advantageous to produce information from original data sources.

In this study, new cartography was developed using satellite images and aerial photographs of the area, employing remote sensing techniques. Several studies have demonstrated that image classification is an effective tool for generating spatial data, as it provides accurate and reliable information (He et al., 2021; Li et al., 2021; Wang et al., 2021).

For the year 1971, cartography was created through photointerpretation of aerial images obtained from INEGI (2023b), while for 2021, Sentinel-2 images were utilized. The Sentinel images were processed using Erdas Imagine 2015 software with a supervised classification method based on the Maximum Likelihood algorithm (Duarte et al., 2016). Subsequently, land cover and land use were classified according to the system proposed by Miranda and Hernández (1963) and the land use and vegetation map from INEGI (2023b). For the study area, nine classes were generated: agricultural land, water bodies, coniferous forest 1 (Abies religiosa [Kunth] Schltdl. & Cham., Pseudotsuga menziesii [Mirb.] Franco, and Cedrus sp. Trew), coniferous forest 2 (Pinus hartwegii Lindl. and Juniperus deppeana Steud.), bare soil, sand mine areas, grassland, alpine grassland, and human settlements (INEGI, 2023a).

The second map was produced through the photointerpretation of aerial images from 1971 flight missions. These images were corrected and processed using Agisoft software (2023) to generate a composite image of the study area. This image was then digitized in QGIS (2023) based on the previously mentioned classification schemes and exported in TIF format. Spectral signatures were verified using ground-truth sampling points, supplemented with data from the Comisión Nacional para el Conocimiento y Uso de la Biodiversidad (CONABIO, 2023) plant species database. Subsequently, the results were refined using satellite imagery available via Google’s web-based image service (2023), corresponding to the year of study. Finally, the processed image was converted to a TIF file, from which a 1:10,000-scale map of the study area was generated. Before importing the data, a topological cleaning process was performed to eliminate digitization errors or overlaps between land cover classes (Cooper et al., 2021; He et al., 2021; Li et al., 2021).

Following the standardization of the files, an overlay analysis was performed using the TerrSet software (Clark Labs, 2020), specifically employing the Land Change Modeler module to compare two images and assess spatial overlap. Additionally, a CrossTab analysis was conducted using data from the two study periods to generate raw data for a land use change matrix, based on the methodology proposed by Pontius et al. (2004). This matrix produces a cross-tabulation of values, from which formulas were used to identify losses, gains, swaps, total change, and net change, as well as their respective rates of change.

The rate of change (rc, %) is a key indicator used to measure both the magnitude and speed of landscape transformations. Its analysis is essential for understanding the dynamics, proportions, and potential impacts of territorial change (Velázquez et al., 2002). The rc was quantified by analyzing vegetation cover dynamics in relation to space and time. The Food and Agriculture Organization of the United Nations and the United Nations Environment Programme (FAO & PNUMA, 2020) propose the following formula to express the percentage change relative to the area in the initial year:

where,

𝑆1 = area at the initial date (ha)

𝑆2 = area at the final date (ha)

t = number of years between the initial and final dates

Traditionally, the location and quantification of land use changes are carried out through cartographic overlay and the generation of a transition matrix. This process produces maps and data that enable the identification of the magnitude and spatial distribution of land cover change dynamics.

Pontius et al. (2004) proposed a transition matrix-based method that analyzes land cover changes from a general level down to more detailed categories. This approach breaks down observed changes into categories such as gains, losses, swaps, persistence, net change, and total change, thereby facilitating comparison and evaluation of the degree of land cover persistence over time. Based on the identified changes, it is also possible to determine whether the gains and losses follow a random pattern. This is assessed by calculating the difference between the observed and expected percentages, as well as the magnitude ratio. A positive difference indicates that land cover categories at time 1 (t1) experienced greater-than-expected losses under a random process, while a negative difference suggests lower-than-expected losses (Pontius et al., 2004).

The transition matrix numerically represents the changes that occurred during the study period, expressed in relative values corresponding to the percentage of each land cover category and its transitions to other categories. This is achieved by comparing the original area from the initial date with the area observed on the final date, along with their respective transitions, as described by Bocco et al. (2001).

One of the change indices proposed by Pontius et al. (2004) is gain (𝐺𝑖𝑗), which indicates the proportion of the landscape that experienced an increase between the initial time (t1) and the final time (t2). Gain is calculated as the difference between the total area of a category at t2 (𝑃+𝑗) and its persistence (𝑃𝑗𝑗), where persistence refers to the area of each land cover class that remained unchanged over the study period. These persistent values are located along the main diagonal of the transition matrix.

Loss (𝐿𝑖𝑗) represents the proportion of the landscape that decreases between two time periods and is calculated from the difference between the total row of t1 (𝑃𝑗+) and the persistence value (??𝑗𝑗) (Pontius et al., 2004).

Swap (𝑆𝑗) between categories is a concept that simultaneously involves both gain and loss of a land-use category over a given land cover. It occurs when the spatial location of a land-use category changes between two dates while its overall area remains constant. Therefore, for every unit of gain in one category, there is an equivalent amount of loss in another. Swap is calculated as twice the minimum value between the gains and losses (Pontius et al., 2004).

Total change (𝐷𝑇𝑗) is defined as the sum of all proportions of the landscape that changed, whereas net change (𝐷𝑗) refers to the definitive change in the landscape. It represents the difference between total change and swap, expressed in absolute values (Pontius et al., 2004).

Results and discussion

The first phase of the results involved the development of two maps corresponding to the years 1971 and 2021 (Figure 2), with the objective of comparing changes in vegetation cover and generating the land-use change matrix (Table 1). These maps classified nine types of land cover, enabling precise categorization and detailed analysis of tree vegetation, specifically within the APFFNT study area. The maps not only reflect vegetation dynamics but also serve as a foundation for generating the necessary inputs for an in-depth analysis of land use and forest cover in the region.

Table 1. Land Use Change Matrix in the Nevado de Toluca Flora and Fauna Protection Area during the 1971-2021 period

This methodological approach, based on the use of satellite imagery and photointerpretation, is similar to that used by Hernández-Pérez et al. (2022) in their study on land-use change in Veracruz, where maps from different dates were also employed to assess landscape transformation.

After generating the land cover classifications for the study area, transition matrices were constructed for each period with the aim of quantifying exchange trends among categories and characterizing the spatial and temporal patterns of landscape transformation. Each interval showed particular dynamics, such that each land cover type demonstrated specific behavior both in general terms and in the detailed analysis of observed trends. Within this framework, Table 1 summarizes the information related to persistence and transitions of the land covers. The values located along the main diagonal represent the persistence of each category, while the remaining cells reflect conversions and exchange flows between different land uses. Notably, the swap between coniferous forest and agricultural and grassland areas increased from 34 954.03 ha in 1971 to 39 911.74 ha in 2021. These dynamics have also been observed in other geographical contexts, such as the Republic of Congo (Mangaza, 2022) and the Mediterranean region (Coskuner, 2022), where transitions from agricultural use to other land covers were identified; however, such behaviors may vary depending on the natural and anthropogenic factors that influence them.

Another result was the rate of change for each land cover type, which allowed for the identification of land uses exhibiting similar behaviors over time. In other cases, declining trends were observed, as in the case of grasslands, while coniferous forest areas showed an increase in surface area. This is consistent with the findings of Pérez Hernández et al. (2021), who reported similar trends in their study on forest cover change in Oaxaca.

Table 2 clearly shows the growth trend in coniferous forests (types 1 and 2), which exhibit a rate of change of 0.57 % during the period analyzed, equivalent to a gain of 4 957.71 ha. These data suggest that the protection strategies implemented in the area have yielded positive results and have contributed to the significant expansion of this forest cover. An important aspect is identifying the type of land cover from which this gain originates; it is observed that grasslands have experienced a loss in their initial coverage, with a negative rate of change of 1.26 %. The same grassland loss trend was reported by Pineda (2017) in his study on land use and dynamics in Hidalgo.

Table 2. Land cover change rates by period in the Nevado de Toluca Flora and Fauna Protection Area.

During the study period, the pine and táscate coniferous forest has experienced the greatest growth, expanding by more than 3 307 ha, primarily in areas that were previously grasslands. This change has been driven by both reforestation efforts and natural ecological succession. The replacement of grasslands by coniferous forests improves ecosystem recovery and could have positive implications for biodiversity and water balance in the region. Increases have also been recorded in other land covers, such as the coniferous forest of sacred fir, ayarín, and cedar (1 380 ha), as well as in agriculture (253 ha). Urban areas have increased by 0.28 % within the boundaries of the APFFNT, which corresponds to an increase of 148.92 ha. It is important to note that a sand mine has been detected in the protected area, reflecting management deficiencies in the region, impacting approximately 15 ha (Table 2).

Another point of interest is the behavior of urban or construction land covers and sand mine, which together have undergone significant transformation with a rate of change exceeding 30 %, corresponding to an approximate increase of 160 ha. This phenomenon is particularly important since these activities are anthropogenic in origin and are occurring within a conservation area where, ideally, such settlements and activities should not take place. The presence of these transformations underscores the need to evaluate and strengthen management and conservation policies in the area to prevent the expansion of activities incompatible with ecosystem conservation.

To more precisely identify changes in land cover, an analysis was conducted using the change indices proposed by Pontius et al. (2004) through a period matrix that allowed visualization of the most relevant results, such as totals, gains, losses, swaps, total change, and net change (Table 3). The use of this type of analysis is highly recommended by studies such as those by Pineda Jaimes et al. (2009), Sandoval-García et al. (2021), and Hernández-Cavazos et al. (2023), who also applied Pontius’s method in their research to obtain detailed results on land use and cover dynamics. This methodological approach enables a more comprehensive and precise assessment of land cover transformations and facilitates the interpretation of trends and change patterns that may be crucial for management and conservation.

Table 3. Transition matrix for the 1971-2021 period in the Nevado de Toluca Flora and Fauna Protection Area.

The transition matrix in Table 3 reflects the changes that occurred during the period from 1971 to 2021 and highlights the land cover types with the highest percentage of area. Coniferous forest cover (types 1 and 2) shows the greatest increase, accounting for 74.48 % of the total area, primarily composed of pine, táscate, oyamel, ayarín, and cedar species. On the other hand, 16.83 % of the area corresponds to anthropogenic activities such as agriculture, mining, and human settlements. The remaining 8.69 % is represented by other non-invasive land covers that have undergone fewer modifications. This distribution pattern underscores the interaction between natural processes and human interventions throughout the study period.

The gains and losses in land cover reflect significant trends in land use behavior, which can be observed in the results for total change and net change. In particular, it was found that coniferous forests (types 1 and 2) had the greatest changes, accounting for 7.84 % of the total area-a trend consistent with the findings of Cruz-Huerta et al. (2015) in their study on land use modeling in the Chignahuapan-Zacatlán region, Puebla. Grassland, on the other hand, was identified as one of the most dynamic land cover types during the analyzed period, showing a decrease of 4 694 ha, equivalent to 8.76 % of the APFFNT. Considering grassland persistence in relation to its dynamics, the total change corresponds to 6.08 % of the existing 16.89 % land cover in the area. This figure reflects the extent of the dynamics associated with the sum of all land cover types that experienced transformations (Table 3). Moreover, the net area change drops to 9.56 %, as some grassland areas have been lost while others, such as coniferous forests, have expanded-revealing a complex transformation dynamic that characterizes the study area.

One of the most significant land cover changes between 1971 and 2021 was observed in the coniferous forests dominated by P. hartwegii and J. deppeana (táscate), which expanded by 6.67 % of the total area of the APFFNT, equivalent to approximately 3 577 ha. This increase in forested areas is mainly due to the conversion of grassland areas into forests, a phenomenon also noted by Pineda Jaimes et al. (2009) in their study on land use change in Estado de México. The expansion of coniferous forests has been supported by different reforestation programs implemented during the analyzed period, which have significantly contributed to the growth of these forested areas (Figure 3), reflecting efforts to mitigate environmental degradation and promote ecosystem conservation in the region.

The expansion of coniferous forests in the Nevado de Toluca has been driven by forest repopulation programs (Figure 3), although outcomes vary according to the study. The region's forest restoration efforts began with a logging ban from 1947 to 1970, followed by the establishment of PROTINBOS (Protectora e Industrializadora de Bosques) in 1970, which managed the sustainable exploitation of forest resources. However, restoration efforts during the 1970s and 1980s were limited, leading to the dissolution of PROTINBOS and its replacement by PROBOSQUE (Protectora de Bosques) in the 1990s. Despite the implementation of new reforestation programs, the overall impact has been moderate, constrained by urban expansion, unrestricted vehicular access, and agricultural pressure in surrounding areas.

To accurately identify land use changes, a detailed map was developed (Figure 4) to clearly and visually pinpoint the transitions. Notably, a higher incidence of transformations is observed in the northern area of the Nevado de Toluca, as well as in transitional zones between different land cover types, where the boundaries between ecosystems or land uses appear to undergo more pronounced alterations. These changes, which are especially evident in key areas, may reflect dynamic and complex processes such as urban expansion, agricultural activity, or climate change, all of which significantly influence the evolution and transformation of the landscape.

Quantifying land cover changes through a map, complemented by a behavior matrix, organizes the information in a clearer and more structured manner, allowing for more accurate and realistic comparisons. This methodology enables the efficient identification of processes such as permanence, degradation, deforestation, or forest cover recovery, which are essential for understanding the current land use dynamics in the area.

The integration of both approaches provides a more detailed and accurate view of the ongoing changes, enabling a better interpretation of landscape interactions and transformations. Furthermore, this approach allows for the identification of critical areas undergoing significant modifications, thus facilitating informed decision-making for land management and conservation. In this way, not only can the processes of change be better understood, but their potential impacts and future dynamics can also be anticipated.

Regarding land use changes, Figure 5 shows the losses and gains in the study area. A notable reduction in grassland is observed, with minimal gains compared to coniferous forests. While coniferous forests show some loss in area, this is offset by gains in other regions. Agricultural land, on the other hand, shows variation in coverage; however, the total change is minor due to nearly equivalent losses and gains-consistent with the findings of Camacho-Sanabria et al. (2015), who modeled land cover and land use changes in the region.

Studies on land use change are essential for identifying trends and patterns in land cover dynamics. These changes are driven by a combination of physical, social, and, in many cases, economic factors (Brovelli et al., 2020). Understanding the magnitude and dynamics of each land cover type is crucial for assessing their impact on the area. This information is presented in Table 4, which highlights positive outcomes-namely, the favorable expansion of forest cover types, now encompassing approximately 74 % of the area. However, anthropogenic interventions, including agriculture and urbanization, remain a concern, as they now occupy nearly 17 % of the territory.

Table 4. Dynamics and magnitude of environmental impact on vegetation covers in the Nevado de Toluca Flora and Fauna Protection Area.

A notable phenomenon in this area is reforestation, which has led to a change in the spatial organization of the landscape. The main consequence of this phenomenon is the success of the implemented conservation plans, which have positively impacted the area’s recovery. In the past, numerous gaps existed within the forests, but over time, these have been filled, demonstrating the benefits derived from conservation strategies. However, deficiencies in oversight have also been identified, particularly regarding the establishment of mines in the area (15 ha), which requires an analysis of the factors that contributed to this event. Despite this, no significant vegetation losses have been recorded since 1971; rather, forest cover has increased over the past 50 years.

It is important to remember that PNAs are key instruments for the planning and regulation of activities in these zones, both in terms of management and administration. The use of technological tools can be a crucial support for decision-making in specific management strategies, as each zone presents characteristics that require differentiated approaches.

Despite the progress achieved in the forest restoration programs in the Nevado de Toluca, results remain limited due to pressure from urban expansion, agricultural activities, and unregulated access to the park. These factors continue to undermine the effectiveness of the implemented strategies, highlighting the need to adopt more comprehensive and sustainable approaches to ensure the long-term conservation of this important natural area.

Conclusions

The analysis of the Nevado de Toluca Flora and Fauna Protection Area (APFFNT) reflects land use changes that occurred over the 50-year period between 1971 and 2021. Although changes in the area are not extreme, a clear trend toward landscape transformation is evident. An example of this is the northern region of Nevado de Toluca, where there has been a notable loss of grasslands, which have decreased by approximately 8.76 % of the total APFFNT territory, equivalent to about 4 694 ha. This type of cover has been largely replaced by areas of tree vegetation, primarily pine and táscate, whose extent has increased by 6.67 % (3 577 ha). Additionally, urban areas have experienced gradual growth, increasing by 0.28 % of the territory, equivalent to 148.92 ha.