Long-Term Land-Use/Land-Cover Change Increased the Landscape Heterogeneity of a Fragmented Temperate Forest in Mexico

: The temperate forests of northern Mexico possess a great diversity of unique and endemic species, with the greatest associations of pine-oak in the planet occurring within them. However, the ecosystems in this region had experienced an accelerated fragmentation process in the past decades. This study described and quantiﬁed the landscape fragmentation level of a degraded watershed located in this region. For that, data from the Landsat series from 1990, 2005 and 2017, classiﬁed with the Support Vector Machine method, were used. The landscape structure was analyzed based on six metrics applied at both, the landscape and class levels. Results show considerable gains in surface area for the land use land cover change (LULC) of secondary forest while the Primary Forest (PF) lost 18.1% of its area during 1990–2017. The PF increased its number of patches from 7075 to 12,318, increased its patch density (PD) from 53.51 to 58.46 # of patches/100 ha, and reduced its average patch size from 39.21 to 15.05 ha. This made the PF the most fragmented LULC from the 5 LULCs evaluated. In this study, strong ﬂuctuations in edge density and PD were registered, which indicates the forests of northern Mexico have experienced a reduction in their productivity and have been subjected to a continuous degradation process due to disturbances such as ﬁres, clandestine and non-properly controlled logging, among others.


Introduction
Forests are one of the most widely distributed natural resources on Earth [1]. These are the support of a wide range of environmental services for socioeconomic development, environmental protection and human well-being [2]. However, non-properly controlled timber production, increment of the surface area for agriculture, cattle grazing and the expansion of urban settlements, threaten forest ecosystems in many regions of the world [3]. Habitat fragmentation and forest degradation are some of these major threats to forest ecosystems [4]. These threats cause multiple alterations, such as modification of the biogeochemical cycles, soil degradation, loss of species, microclimate and landscape structure changes, among others [5][6][7][8].
Habitat fragmentation is defined as the division of naturally homogeneous areas into smaller fragments [9]. The division of these areas leads to the formation of isolated patches [10]. Isolated patches limit the landscape connectivity, affecting energy exchange among them and threatening the ecosystems' health [11][12][13][14]. Even though natural ecosystems possess some degree of fragmentation due to natural interactions within them [15], these fragmentation processes were accelerated during recent years in both, space and time, mainly caused by anthropogenic activities [16][17][18][19]. Investigating the pattern and process of change in the forest landscape is essential to understand its fragmentation level. That can serve as a decision tool to mitigate forest degradation and promote forest conservation.
Remote sensing data, together with geographic information systems (GIS), have been widely used to spatiotemporally describe and quantify habitat fragmentation [20]. These tools have served to analyze the pattern dynamics of the forest ecosystem [21]. Remote sensing provides a fast and cost-effective information to quantitatively estimate the biophysical characteristics of the forest landscape [22]. It also monitors the forest landscape over a spatial-temporal scale [23][24][25]. This feature makes it a valuable data source, which has an appropriate accuracy level, to study the landscape structure [26].
In spite the easiness of access to remote sensing data from Landsat, the reliable and accurate production of land cover/land use (LULC) maps derived from satellite images, remains a major challenge [27]. In this regard, the use of automatic learning-based classifiers, such as Support Vector Machines (SVM) [28], represent an alternative and have shown promising results.
Over the past 30 years, several metrics to assess the landscape fragmentation have been proposed [29]. Landscape metrics are commonly used tools to evaluate the ecosystem's structure and composition [30]. They are useful indicators in the study of the ecosystem's health [31]. Previous research has widely used these metrics [32][33][34][35], currently considering them a standard data source in landscape research at different scales [36,37].
In Mexico, temperate forests cover nearly 32 million hectares, possessing a great diversity of unique species of a great biological importance [38]. This region is home of the planet's largest association of pine and oak species [39]. In particular, the 'Sierra Tarahumara' in northern Mexico is characterized by its high biodiversity and its high number of endemic species. It is also recognized by the International Union for the Conservation of Nature as a megacenter of plant diversity (i.e., 4000 plant species) [40]. However, this region has experienced a gradual process of deforestation and is one of the frequently disturbed forest regions of the state of Chihuahua [41]. In addition, research on the fragmentation level of the temperate forests of Chihuahua is lacking.
Thus, this study investigated the LULC dynamics and in its implications on the fragmentation level of a temperate forest ecosystem located in Chihuahua, Mexico. The analysis was performed at the basin level. The period evaluated comprehended 27 years and the LULC classification was based on SVM.

Study Area
The study area is located in the southwest region of the state of Chihuahua, in the 'Sierra Tarahumara', covering a surface area of ≈ 497,159 ha. It is located within the coordinates 108 • 0 W, −28 • Figure 1). The relief of the study area consists of an extensive system of deep canyons, producing a mix of temperate and tropical ecosystems, with semi-cold and semiwarm type climates, respectively. This ecosystem comprises sites where pure pine, pure oak or mixed pine-oak forests predominate. The 'Sierra Tarahumara' presents relicts of ancient forests coming from the temperate forests of North America, which have been centers of oaks diversification. The mean annual precipitation is between 1052 mm to 2348 mm. The mean annual temperature fluctuates from 10 • C to 11.7 • C [42].
Productive activities in the region include forestry, extensive cattle ranching and irrigated agriculture [43]. With the rapid economic development experienced in recent years, the 'Sierra Tarahumara' has been used for extensive agriculture and cattle grazing. That has increased the deforested area, which has impacted the soil, making it vulnerable to water erosion. This condition has become critical in some mountain areas, where the rocks have emerged and where the soil has been considered completely lost [44].
The reduction in natural vegetation has brought several ecological problems, such as the decrease of forest species, and the increase of floods' rate and duration with a weakened capability of the forest to adjust runoff and reduce floods. Thus, it is critical to investigate the pattern and dynamics of the forest landscape structure to contribute to a better decision-making and to improve the ecological management in this region. Productive activities in the region include forestry, extensive cattle ranching and irrigated agriculture [43]. With the rapid economic development experienced in recent years, the 'Sierra Tarahumara' has been used for extensive agriculture and cattle grazing. That has increased the deforested area, which has impacted the soil, making it vulnerable to water erosion. This condition has become critical in some mountain areas, where the rocks have emerged and where the soil has been considered completely lost [44].
The reduction in natural vegetation has brought several ecological problems, such as the decrease of forest species, and the increase of floods' rate and duration with a weakened capability of the forest to adjust runoff and reduce floods. Thus, it is critical to investigate the pattern and dynamics of the forest landscape structure to contribute to a better decision-making and to improve the ecological management in this region.

Framework
Firstly, the study focused on the standard processing of the Landsat images. Secondly, the LULC classification of the study area was carried out by using the Support Vector Machine classifier (R Foundation for Statistical Computing, Vienna, Austria). Thirdly, the forest landscape analysis was based on the result of the supervised classification, resulting in a set of entities or patches with LULC attributes (Figure 2).

Framework
Firstly, the study focused on the standard processing of the Landsat images. Secondly, the LULC classification of the study area was carried out by using the Support Vector Machine classifier (R Foundation for Statistical Computing, Vienna, Austria). Thirdly, the forest landscape analysis was based on the result of the supervised classification, resulting in a set of entities or patches with LULC attributes (Figure 2).

Data Collection
Three satellite images from the Landsat sensor were used (Path 33, Row 41). Two of them were from the Thematic Mapper (TM) 5 and one was from the Operational Land Imager (OLI) 8. The spatial resolution of the satellite images was 30 m × 30 m. The images are from the years 1990, 2005 and 2017; they were acquired cloud free, and all of them correspond to the month of March to reduce variability among the images and avoid con-

Data Collection
Three satellite images from the Landsat sensor were used (Path 33, Row 41). Two of them were from the Thematic Mapper (TM) 5 and one was from the Operational Land Imager (OLI) 8. The spatial resolution of the satellite images was 30 m × 30 m. The images are from the years 1990, 2005 and 2017; they were acquired cloud free, and all of them correspond to the month of March to reduce variability among the images and avoid confusion at the time of classification. All the images were downloaded from the United States Geological Survey image server [45] and their characteristics are shown in Table 1. Table 1. Characteristics of the images used in the study.

Sensor
Capture

Data Processing
The scenes were radiometrically corrected. The conversion from digital numbers (DN's) to reflectance values was carried out with the Dark Object Subtract 1 method, which is based on the image's properties and is the most widely used method for change detection [46]. This method is performed to correct the reflectance differences caused by the solar distance and the solar angle [47]. That serves to avoid any discrepancies because the data comes from multi sensors and dates [48,49]. Equations (1) and (2) were used to carry out this conversion.
where: L λ denotes the spectral radiance, QCAL is DN, QCALmin and QCALmax are the minimum and maximum quantized calibrated pixel value, respectively, Lmin λ denotes the spectral radiance scales to QCALmin while Lmax λ represents the spectral radiance scales to QCALmax, p λ denotes the reflectance, d is the distance from the earth to the sun, ESUN λ is the mean solar exoatmospheric irradiance and θ S is the solar zenith angle. The radiometric correction was carried out in the QGis v.3.10 software (QGis Development Team, open source) with the SemiAutomatic Classification plugin developed by Congedo [50]. Landsat 8 images were radiometrically converted by applying Equation (3).
where: p λ is the top of the atmosphere (TOA) planetary reflectance, with correction for solar angle and θ SE is the local sun elevation angle.

Support Vector Machine for Image Classification
The classifier based on machine learning known as Support Vector Machine (SVM) was applied to the 1990, 2005 and 2017 image stacks to obtain the LULC information. The SVM was run using the caret package in the R statistical software (R Foundation for Statistical Computing, Vienna, Austria) ver. 3.4.2 [51]. The SVM is a supervised non-parametric statistical technique. The classification objective of the SVM is to determine a hyperplane to optimally separate two classes. An optimal hyperplane is determined by using training data and its separation capacity is verified with validation data [52]. Therefore, the primary goal of the SVM is to find the optimal separation plane among all possible hyperplanes. The SVM approach uses the principle of structural risk minimization [53]. The three images were classified into five classes of LULC, which were: Areas with No Apparent Vegetation (AWV), areas where vegetation is not present or where vegetation is not perceived by the satellite images; Secondary Forest (SF), forest with a canopy partially covering the ground; Human Settlements (HS), areas where population exist, where buildings and infrastructures for human service are detected; Water Bodies (WB), surface covered temporarily or permanently by water; and Primary Forest (PF), forest with a canopy completely covering the ground.

Classification Accuracy
The K APPA Index was used to assess the accuracy of the LULC classification. The K APPA Index is a discrete multivariate technique used to validate LULC classes through a matrix [54]. It is commonly used to establish the degree of similarity between the simulated mapped values and the real ones. A K APPA value equal to one indicates 100% similarity between mapped and real values. Conversely, a value of zero suggests 0% similarity. The K APPA Index is represented by Equation (4).
where: K APPA is the K APPA Index, k is the total number of rows in the matrix, X ii is the number of observations in row i and column i (along the diagonal), X i+ and X +i are the total marginal for row i and column i, respectively, and N is the total number of observations.

Conventional Estimation of Areas of Land Use/Land Cover
After the classification of LULC, a conventional quantification of areas was carried out. In addition, rates of changes were calculated with Equation (5), which was adapted from [55]: where: r = rate of change (% y −1 ), S 1 = area at time 1, S 2 = area at time 2 and n = the number of years studied (time 1-time 2).

Landscape Metrics
Forest landscape fragmentation has a close relationship with the patch size, isolation, shape, connectivity, and edge [56]. Based on its complex nature, forest landscape fragmentation cannot be measured using a single metric. Thus, a group of landscape metrics is usually employed to describe the forest landscape fragmentation.
This study analyzed the landscape based on metrics applied at both, the landscape (L) and class (C) levels. The metrics show numerical information about the landscape composition, its configuration, and dimension [57,58]. A metric at the landscape level quantifies several map categories or classes, i.e., the set of classes making up the LULC of a region. Meanwhile, a metric at the class level quantifies the specific spatial characteristics of the patches making up one class or several types of LULC in a region. With that, several configuration metrics can be formulated either in terms of individual patches or in terms of the entire class or landscape, depending on the scope of the study.
Despite the existence of countless effects of the forest fragmentation, many of these effects are due to a reduction in the area and connectivity of the forest landscape. The forest fragmentation is also due to an increase in the proportion of habitats influenced by the edges. For example, the perimeter-area fractal dimension is a measure of the complexity of the patch shape, which can be calculated for each patch and then averaged either at the class or landscape level. The percent of landscape describes the patch size and quantifies the edge created by these patches, which are useful indicators in landscape fragmentation analyses. The analysis of these metrics plays a key role in the process of forest loss and fragmentation: forest loss and fragmentation generally involve the disaggregation of continuous patches of forest into more dispersed and/or separate (i.e., subdivided) habitats and more isolated patches [59].
For each LULC cover map, the metrics were calculated in a raster format with the FRAGSTATS software, ver. 4.2© (UMASS, Amherst, MA, USA; https://www.umass.edu/ landeco/research/fragstats/fragstats.html (accessed on 15 June 2021). Such metrics are described in Table 2. Metric selection was based on the studies by Del Castillo et al. [60] and Tian et al. [61]. The analysis of forest fragmentation measures the configuration of a particular type of patch and the connectivity among the patches based on their neighborhood.

Land Use/Land Cover Classification of 1990, 2005 and 2017
Five types of LULC were obtained for the study area, as they are the most representative for the region (Figure 3). They included areas with no apparent vegetation, secondary forest, human settlements, water bodies and primary forest. According to the classification method applied for the dates evaluated in this study (1990, 2005 and 2017), the two LULC with the greatest surface area were secondary forest and primary forest. According to the K APPA Index, the accuracy was higher for the classifications corresponding to 2005 and 2017, both with overall coefficients of 0.85, than that of 1990, which obtained a coefficient of 0.84 (Table 3).

Patterns of Land Use/Land Cover Change
Results from the analysis of land use/land cover change (LULCC) showed considerable gains for the secondary forest (Table 4). Meanwhile, the surface area of the primary forest decreased from 55.8% in 1990 to 37.7% in 2017. The areas with no apparent vegetation increased from 4.11% to 4.87% during the period 1990-2017. The categories of human settlements and water bodies showed an increase from 0.03 and 0.01% in 1990 to 0.10% and 0.03% in 2017. In general, the primary forest was the only land use showing a decreasing trend. The LULCs of secondary forest, human settlements and water bodies showed the highest change rates, with 0.87%, 0.98% and 1.07%, respectively, during 1990-2005 and 3.37%, 10.98% and 5.80%, respectively, during 2005-2017.

Patterns of Land Use/Land Cover Change
Results from the analysis of land use/land cover change (LULCC) showed considerable gains for the secondary forest (Table 4)

Metrics at the Class Level
The metrics for the LULC evaluated, at the class level, are shown in Table 6. From 1990 to 2017, the metric of percentage of landscape, corresponding to the portion (%) of the total area occupied by each class, showed an increment for the classes of areas with no apparent vegetation, human settlements and water bodies. In contrast, the class of primary forest showed a considerable decrement in 2017. Regarding the metric of patch density, it showed an increment for the LULC of areas with no apparent vegetation, where patches increased during 1990-2017. The LULCs of areas with no apparent vegetation and human settlements showed and increment in their number of patches from 3095 to 5577 and from 6 a 43, respectively. This increasing trend also occurred for patch density and percentage of landscape. Conversely, mean patch size and the Euclidean nearest-neighbor distance showed a negative trend. Regarding the LULC of water bodies, its patch density, percentage of landscape, mean patch size, number of patches and patch cohesion index showed a positive trend. In the case of the secondary forest, it number of patches got reduced from 13,035 to 8761. This same trend occurred for percentage of landscape and Euclidean nearest-neighbor distance. Conversely, patch size, patch density and patch cohesion index showed an increment during 1990-2017 for this class.

Land-Use/Land-Cover Classification
The K APPA Index was greater than 0.85 for the LULC use maps. That is considered as a high accuracy by Manandhar et al. [63] and good enough to perform an analysis of the study area, since they exceeded the minimum value stipulated by Anderson et al. [64] for LULC use maps derived from satellite images. This same threshold has been adopted in previous studies [65][66][67].
Regarding the changes in LULC, the results of this study agree with previous ones. For example, Costa et al. [59] analyzed forest cover changes during the periods of 1976-1996 and 1996-2014. They found that, in the second period, there was a reduction in deforestation, where a regeneration processes occurred, increasing the patch size and reducing the number of patches. This may positively have affected biodiversity, considering that larger and more compact patches may offer better habitat quality for forest-dependent species [68]. Gounaridis et al. [69] used Landsat images to analyze 28 years of forest cover change where they found reductions of 2-3% in the surface area. This shows an accelerated change in the forest structure during the same period. Rosa et al. [70] evaluated forest cover changes during 27 years from Landsat images. The forest cover reported in one of their study areas at the end of the studied period was 31%. In this study, the primary forest covered 37%.

Metrics at the Landscape and Class Levels
In this study, landscape quantitative metrics measuring LULC patterns and changes were applied to infer about the fragmentation of the landscape. The analysis of the landscape metrics was performed at both, the landscape and the class levels. The results from such analysis provide an in-depth understanding of key trends in the LULC structure.
At the landscape level, the overall increase in the number of patches, as well as in the patch density, indicated an increment in the fragmentation level during the evaluated period [37]. The reduction in the mean patch size also indicated forest fragmentation, as pointed out in previous research [71]. Patch density and mean patch size are complementary, since an increase trend in patch density and a decrease trend in mean patch size confirm the landscape has undergone a process of fragmentation [72]. These results are consistent with those reported by Fuller [73], who used Landsat images in Virginia, USA during the period 1973-1999 to analyze forest fragmentation, and those reported by Vázquez et al. [74], who analyzed fragmentation in forests of Durango, Mexico during 1974-2011. Results from an additional metric, indicating that a fragmentation process occurred in the study area, was edge density, which registered an increasing trend comparable to the fragmentation rates in a study shown by Cottam et al. [75]. The increase in edge density, together with the patch density, confirms the landscape has been experiencing a process of fragmentation [76][77][78][79][80][81]. Regarding the metric of perimeter-area fractal dimen-sion, its value fluctuated over time, indicating that landscape changes are based on large alterations in both, quantity and shape [38]. These fluctuations in the perimeter-area fractal dimension may be associated to anthropogenic interventions, which modify the shape of patches, making them more irregular, due to the effects associated with deforestation caused by agriculture and the establishment of induced grasslands [77]. Heterogeneity, which was assessed through the Shannon diversity index, did not change notoriously during the 27 years studied, starting with values between 0.83, increasing to 0.86 and then decreasing to 0.84. However, these fluctuations may be caused by forest management. Decisions regarding collective forest management are often subjected to arbitrary changes due to the transition of local chiefs, in general, resulting in lack of sustainable forest management [78].
At the class level, the primary forest increased its number of patches, patch density and Euclidean nearest-neighbor distance. Conversely, it reduced its mean patch size and percentage of landscape. The increase in patch density and number of patches explains the fragmentation process in this class [79]. If the metric of percentage of landscape decreases with time while patch density increases, a fragmentation process is occurring. This same behavior between percentage of the landscape and patch density was reported by Zhou and Wang [80].
The increase on the metric of Euclidean nearest-neighbor distance indicates that large patches of primary forest are in a process of isolation, as mentioned by Tolessa et al. [81]. The secondary forest showed a reduction in its number of patches and an increment in patch density. The LULC of areas with no apparent vegetation, human settlements and water bodies registered an increment in the number of patches, patch density, percentage of landscape and a reduction in mean patch size. That indicates the creation of new patches and a fragmentation process in the class [73]. Regarding the Euclidean nearest-neighbor distance, it showed a reduction for all the LULCs, indicating a reduction in the patch distance of these classes [82].
The patch cohesion index showed very similar values for all the classes, only with small variations. In previous studies, connectivity was measured using these same metric [83,84]. The areas with no apparent vegetation were the only class showing a decrease in the values of patch cohesion index, indicating this class had the lowest physical connection among the patches of the same class.
The trend towards the dominance of the secondary forest is evidence of a transition process, in which the primary forest becomes secondary forest in most of the cases. However, the primary forest also experienced transitions to areas with no apparent vegetation and human settlements. This is an indicator of forest fragmentation, where continuous and long fragments of primary forest became more compacted, isolated while new small fragments of areas with no apparent vegetation and secondary forest appeared. The main driving forces of this transition, which led to a landscape modification, are development of rural areas [85], as well as inappropriate land use practices and policies [86,87]. In addition, forest fires are sometimes poorly prevented and controlled, causing huge damages to forests. Fires are generally caused by excess of biomass, which indicates appropriate forest management is lacking [88]. In addition, a recurrent drought of more than a decade in the study region has played a role as a driving force, producing appropriate conditions for forest fires, pests and diseases. Deforestation for agricultural activities, has also been an important driving force. The aforementioned has caused a change in the structure of the temperate forest in the region, transitioning its growth mainly towards a secondary forest.
Finally, this study uniquely employs methods for the analysis of landscape fragmentation through classifying LULCs based on a machine learning-based technique (i.e., SVM). The combination of this classification method and the fragmentation analysis was useful to infer about the landscape structure. Moreover, this study is the first one of its kind for the study area.

Conclusions
In this study, a novel approach for the analysis of the temperate forest fragmentation in the 'Sierra Tarahumara' was performed by classifying LULCs with the machine learningbased technique SVM, resulting in LULCs maps of a high accuracy. The study area experienced a fragmentation process, notoriously impacting the primary forests. The patches of primary forest got reduced and more isolated while other classes, such as the secondary forest, increased in their covered area.
Information from the metrics at the landscape and class levels could be applied in forest local management planning, as a tool for conservation and protection of natural resources. This type of information can be used to improve policies, promote sustainable forest management, public awareness and scientific-based decision making.