Study area

The study area is the FPARSQ, situated at the border of Estado de México and Querétaro, encompassing the municipalities of Amealco, Aculco, and Arroyo Zarco. It is bounded by coordinates 20° 6' 14.70" N - 99° 56' 13.79" W and is situated at an elevation of 2 622 meters (^{Google Earth, 2021}; Figure 1). The nearest locality is San Ildefonso Tultepec, with a total population of 3 632 residents, of which 40 % belong to the Otomí indigenous group: 1 891 women (799 of whom speak an indigenous language) and 1 741 men (692 of whom speak an indigenous language) (^{Instituto Nacional de Estadística y Geografía [INEGI], 2020}).

Source: Compiled by the author based on INEGI data (2016).

Figure 1 Location of the study area in the Forest Protected Area of the southern rivers of Querétaro.  

Cartography, Land Use, and vegetation types

Cartography and land use/vegetation data layers (LuV) from Mexico were used, including data from 1968-1986 (Series I) published in 1997, 2007-2008 (Series IV) published in 2009, and Series VI published in 2017 but with 2016 information. All series were used at a 1:250,000 scale. The data source for all series was the Comisión Nacional para el Conocimiento y uso de la Biodiversidad (^{CONABIO, 2021}).

The study area, encompassing the FPARSQ, was delineated using Google Earth (kml format) and converted to vector shapefile format using ArcGIS v. 10.7.1 (^{ESRI, 2011}). All subsequent processes were carried out using this GIS due to its practicality, consistent use, and broad acceptance in the country.

Estimating changes in land use and vegetation types

Based on the area of each land use (ha and %), extracted from the cartography, three maps for the years 1997, 2009, and 2016 were created. The study considered two periods: 1) 1997-2009 and 2) 2009-2016. The rate of change in land use was calculated using the following formula from the FAO (^{Food and Agriculture Organization of the United Nations, 1996}):

where,

cr = change of rate (%)

a ₁ = area at the initial date

S ₂ = area at the final date

n = Difference in years between the initial and final dates.

For both periods, the corresponding maps were overlapped, and cross-tabulation value matrices were obtained as proposed by ^{Pontius et al. (2004)} (Appendix 1).

Classification of Land Use Changes

Land use and vegetation (LuV) changes were classified into four categories based on the criteria set by ^{Palacio-Prieto et al. (2004)} and the transition matrix (Appendix 2) (^{Pontius et al., 2004}):

a) No Change: Areas where the land use remained the same in the final year as in the initial year, represented along the main diagonal.

b) Deforested or Degraded Forest Loss: The first category refers to changes from tree-based land uses (forests and rainforests) to cultivated areas, induced pastureland, and croplands. Degraded forest loss includes the loss of non-tree vegetation cover (such as other types of vegetation, hydrophytic vegetation, scrubland, and natural grasslands) due to conversion to induced and cultivated grasslands, crops, and other land covers.

c) Replacement Vegetation: Areas that were previously used for agricultural activities (such as croplands or induced and cultivated grasslands) and are now covered by forest, rainforest, scrubland, hydrophytic vegetation, or natural grasslands.

d) Productive Activity: Amount of area in land uses where productive human activities (induced pastureland or agriculture) are observed.

Trends in Land Use and Vegetation Changes

A discrete-time Markov chain is a sequence of random variable X1, X2, …, Xn, where the probability distribution of the next state Xn+1 depends only on the current state Xn y nd not on previous states Xn-1, Xn-2… (Markov property). The chain moves from one state to another, and the probability p _ij of transitioning from states _i to state s _j in one step, known as the transition probability, can be written as pij=Pr(X1=sj|X0=si). Assuming that the transition probabilities depend only on the time interval between the initial year S ₁ and the final year S ₂ of the analysis, the stochastic process is considered homogeneous (^{Linares et al., 2022}). In this study, Markov chain tools were used to analyze the probabilities of land use and vegetation (LuV) change trends using ArcGIS v. 10.7.1 (^{ESRI, 2011}).

Possible Causes of Land Use Changes

Field surveys were conducted to identify potential causes of land use changes in affected areas. Georeferenced control points from these sites were extracted and loaded onto a Garmin Etrex 30 GPS device. The data obtained from the cartography were then verified on-site.

To identify socioeconomic factors as potential causes of these changes, a survey was designed, including closed-ended questions about productive activities related to forestry and other aspects of the population's perceptions (Appendix 3). The sample size for the application of the survey was determined with the following formula (^{León et al., 2016}):

where,

n = sample size sought

N = population size

Z = statistical parameter depending on the confidence level (p = 0.05)

e = maximum accepted estimation error

p = probability of occurrence of the studied event (success)

q (1-p) = probability that the event studied will not happen.

According to the Population and Housing Census of ^{INEGI (2020)}, the municipality of Amealco de Bonfil contains 46 localities, including San Ildefonso Tultepec, which has 10 neighborhoods totaling 864 households. Nearby localities affected by the sites include San Ildefonso Centro, Cuisillo, Xajay (Texquedo, Llano Largo), and Bothe, which collectively comprise 125 households. Based on the calculated sample size, 94 surveys were applied.

Additionally, data from the Comisión Nacional Forestal (^{CONAFOR, 2020}) were reviewed, and a literature review of forest fires in the FPARSQ from 2009 to 2016 was conducted.

Land use, changes and trends in the period 1997-2009

At the beginning of the studied period, in 1997, the following land uses and vegetation types were identified: oak forest, oak-pine forest, pine-oak forest, agricultural use, secondary shrubby vegetation of pine-oak forest, induced pastureland, secondary shrubby vegetation of oak forest, and induced pastureland with erosion. It should be noted that oak forests covered 77 % of the studied area. Table 1 shows the results for this period.

Table 2 shows the LuV change transition matrix. According to ^{Farfán et al. (2016)}, the matrix is an asymmetrical arrangement table: the horizontal axis represents the land use categories from the first date, and the vertical axis represents the categories from the second date. The diagonal of the matrix shows the area of each land use category that remained unchanged during the considered period. On the other hand, the rest of the cells show the area that experienced a change to some other type of category. In the second-to-last column of the matrix, the total areas of all categories from time 1 are shown. The second-to-last row presents the total sum for the land uses from time 2. An additional column has been added to the matrix to represent the losses of each LuV category. Additionally, an extra row has been added to show the gains.

Tables 1 and 2 highlight the change in pine-oak forest, which lost 472 ha, of which, according to the transition matrix, the majority (402 ha) changed to secondary shrubby oak forest vegetation (SSv/QF). Additionally, the oak forest was reduced by 12 %, equivalent to a loss of 958 ha, with transitions to induced pastureland (480 ha) and secondary vegetation (326 ha). Regarding land uses that recorded substantial increases in occupied area, SSv/QF stands out with 9 % (742 ha) and induced pastureland with 6 % (513 ha).

For the period 1997-2009, Table 3 shows the trends in LuV using a probability matrix. According to the results, the LuV with the highest probability of remaining are oak-pine forest, oak forest, induced pastureland and agriculture with 90, 80, 60, 60 and 60 %, respectively. The secondary vegetation of oak forest has the lowest probability of remaining (40 %). On the other hand, pine-oak forest has a high probability (90 %) of migrating to shrub secondary vegetation of oak forest.

Land uses, changes and trend for the period 2009-2016

Table 4 shows the dynamics of LuV for the period 2009-2016. Among the changes observed, the loss of 394 ha of induced pastureland, representing 5 %, stands out. Also, the oak forest had losses that were mainly occupied by secondary vegetation (375 ha) and induced pastureland (40 ha), so in net terms the loss was 68 ha.

Table 5 shows that in the 1997-2009 period, 76 % of the area remained unchanged; however, forest loss was about 631 ha, representing 8 % of the total area. In the 2009-2016 period, the area that remained unchanged stands out, which increased slightly by 5 % compared to the previous period, from 6 085 ha to 6 516 ha. The forest loss is 432 ha less than the loss observed in the period 1997-2009. Also, there was a decrease of 117 ha in replacement vegetation compared to the previous period.

Table 6 shows the probability matrix for the period 2009-2016. The land uses with the highest probability of remaining over time are pine forest, oak-pine forest and oak forest with 100, 100 and 90 %, respectively. Agriculture recorded 50 % probability of remaining and the other 50 % could change to oak-pine forest (20 %) or oak forest (20 %) and pine forest (10 %).

In general, during the 20 years analyzed, three important results can be distinguished in the FPARSQ: 1) the land uses and vegetation (LuV) with tree cover that lost the most surface area in net terms were oak forest with 1 026 ha (16 %) and pine-oak forest, whose 472 ha disappeared from the FPARSQ (100 %); 2) secondary shrub vegetation/oak forest increased from 21 ha to 1 088 ha, representing an increase of over 5 000 %; 3) regarding LuV associated with productive activities, induced pastureland grew by 119 ha (98 %) and agriculture increased by 128 ha (6 %).

Studies related to natural protected areas (NPA) provide different kinds of results. For example, ^{Rioja-Nieto et al. (2015)} found no clear evidence supporting the conservation objectives of the PNA. ^{Blackman et al. (2015)} concluded that PNA do not have clear protection effects on the clearings occurring in them. On the other hand, ^{Morales-Hernández et al. (2016)}, in their study of the Bahía de Banderas PNA, found that there is a process of change and replacement of land cover and uses mainly associated with agriculture at the expense of forests and rainforests.

Mexico has lost large areas of its ecosystems in recent decades; by 1976, the original total coverage in Mexico had been reduced by 38 %, and by 1993 it only covered 54 % (^{Challenger et al., 2009}). A significant portion of these areas belong to rural communities; their land use changes affect biodiversity and reduce vital ecosystem services (^{De Alba et al., 2020}). Under these circumstances, the identification and characterization of vegetation units have been of great importance, and GIS have been widely used in their analysis (^{Flores et al., 2020}). GIS generates data using remote sensors that allow for the integration of forest cover information.

Causes of Land Use Change

Figure 2 shows the analyzed dynamics from 1997 to 2016 and the identified sites with some impact on the cover, which were visited to corroborate the land use change.

Source: Compiled by the authors based on data from ^{INEGI (1997}, ²⁰¹⁶⁾.

Figure 2 Dynamics of change (left) and main affected sites (right) in the Forest Protected Area of the southern rivers of Querétaro. 

A preliminary study was conducted using a survey applied to a sample of 94 households from five localities in the region, aimed at identifying people's perceptions regarding the dynamics of change of LuV. This survey, applied on November 2021, indicated that 67 % stay at home, 11 % are involved in agriculture, 10 % are engaged in artisanal activities, and the rest in other activities. Regarding the energy used for cooking, the responses indicate that 32 % use liquefied petroleum gas (LPG), 30 % use firewood, and 38 % use both LPG and firewood; respondents stated that the firewood is extracted from nearby forests.

Regarding the residents' perception of the possible causes of the impacts, 90 % identify illegal logging as the main problem, which is associated with charcoal production and the extraction of firewood for domestic use. Charcoal production carries a high risk of causing fires. CONAFOR reports that between 2009 and 2016, a fire affected nearly 30 hectares of vegetation in the study area.

The conditions and trends of biodiversity change, as indicated by ^{Morales-Hernández et al. (2016)}, are rooted in social, economic, and political factors known as root causes. These factors often lead to changes in land cover and land use, overexploitation of organisms, introduction of invasive exotic species, anthropogenic climate change, and the addition of polluting products. Consequently, regarding the causes of vegetation changes, it is common for frequently reported factors to include anthropogenic activities and the high content of woody material that is vulnerable to forest fires. Another contributing factor is the expansion of agriculture and livestock farming, which leads to the extension of infrastructure (^{Ellis et al., 2017}; ^{Monjardín-Armenta et al., 2017}).

GIS methodologies have frequently been used in the study of natural protected areas (NPA) in the country. For example, ^{Neri-Suárez (2018)} evaluated conservation gaps (known as gap analysis) in vegetation types relative to PNAs in Guerrero, Mexico. Additionally, ^{Chico et al. (2015)} assessed land use changes using orthophotos from 1994, 1999, and 2007 in the Insurgente Miguel Hidalgo y Costilla National Park, also known as La Marquesa, in Mexico City. ^{Rioja-Nieto et al. (2015)} evaluated the effects of management in El Palmar State Reserve on the Yucatán Peninsula using remote sensing techniques based on hierarchical analysis of habitat characteristics. A similar case, but applying geospatial modeling of agroecological and forestry variables, was conducted by ^{Villegas-Martínez (2021)} in the Otomí-Mexica Park in Estado de México.