Effect of Deforestation on Land Surface Temperature in the Chiquitania Region, Bolivia

: Neotropical forests offer alternatives to surface cooling and their conservation is an effective solution for mitigating the effects of climate change. Little is known about the importance of tropical dry forests for temperature regulation in Chiquitania, a region with increasing deforestation rates. The impact that deforestation processes are having on the surface temperature in Chiquitania remains an open question. This study evaluated trends in forest cover loss based on land surface temperatures ( ◦ C) in forested and deforested areas in Chiquitania. We hypothesized a positive relationship between higher deforestation and a temperature increase, which would decrease the resilience of highly disturbed Chiquitano forests. We evaluated ten sampling sites (10 × 10 km), including ﬁve in forested areas with some type of protection and the other ﬁve in areas with populated centers and accelerated forest loss. We developed scripts on the Google Earth Engine (GEE) platform using information from the Normalized Difference Vegetation Index (NDVI, MOD13A2) and the daytime and nighttime Land Surface Temperature (LST, MYD11A1) from MODIS products for the period 2001–2021. The statistical signiﬁcance of the trends of the time series averages of the MODIS products was analyzed using a nonparametric Mann–Kendall test and the degree of the relationship between the variables was determined using the Pearson statistic. Our results based on NDVI analysis showed consistent vegetation growth in forested areas across the study period, while the opposite occurred in deforested lands. Regarding surface temperature trends, the results for daytime LST showed a positive increase in the four deforested areas. Comparatively, daytime LST averages in deforested areas were warmer than those in forested areas, with a difference of 3.1 ◦ C. Additionally, correlation analyses showed a signiﬁcant relationship between low NDVI values due to deforestation in three sites and an increase in daytime LST, while for nighttime LST this phenomenon was registered in two deforested areas. Our results suggest a signiﬁcant relationship between the loss of forest cover and the increase in land surface temperature in Chiquitania. This study could be the ﬁrst step in designing and implementing an early climate–forest monitoring system in this region.


Introduction
Tropical dry forests are of great ecological importance but are also highly threatened and have received little attention focusing on the impact of their loss in the context of climate change [1]. Forests offer solutions for surface cooling due to the ability of trees to capture and redistribute the energy of the Sun [2]. Forests can exert a strong influence on local land surface temperature based on biophysical mechanisms as they generally have a daytime and nighttime LST (°C) in forested and deforested areas; (c) to establish the relationship between intact forested areas and areas with deforestation processes and daytime and nighttime LST (°C) in Chiquitania. We hypothesized that the increase in temperature at ground level is related to a reduction in tree cover, which would demonstrate that conservation is key to risk and hazard management. This research will help decision makers formulate risk prevention strategies for natural disasters, adapt to climate change, and establish public policies that can help improve land-use planning, thus avoiding the progressive advance of land-use change.

Study Area
The study area includes the Chiquitania region, an extensive flat to undulating plain located in the Department of Santa Cruz, with mountain ranges of varying amplitude, both vertical and horizontal, formed by transverse faults and water erosion. The Chiquitania region contains a rich cultural heritage made up of indigenous Chiquitano, Guarayo, and Ayoreo peoples, as well as Creoles, indigenous peoples from western Bolivia (Quechua and Aymara), and Mennonite settlers [24].
The climate of this region is warm sub-humid tropical with a rainy period in summer and dry weather in winter, with little average annual thermal variability. The distribution of precipitation determines a strong seasonal rainfall regime, with the rainy season extending from November to March, reaching its maximum in January [24]. The average annual monthly temperature is 25.3 °C, with a difference in the latitudinal gradient; high values between September and March, with a maximum in October; low temperatures between June and August, with lower averages in July [24].
In the Chiquitania region, there are three well-defined vegetation types: forests, welldrained savannas, and savannah wetlands [25]. Forest cover is mainly represented by the Chiquitano Dry Forest (Figure 1), which extends between the Amazonian forests to the north and the Chaco forests to the south [25]. The Chiquitano Forests represent the second largest area of tropical dry forests in the Americas [26] and cover an area close to 24 million hectares, 83% of which are located in Bolivia [27]. The most widely distributed forest type is occupied by large areas of well-drained plains with dominant species such as Acosmium cardenasii, Anadenanthera colubrina, Aspidosperma cylindrocarpon, Aspidosperma tomentosum, and Astronium urundeuva, among others [25].

Data Sources
Normalized Difference Vegetation Index. We used the product MOD13A2 Version 6 (https://lpdaac.usgs.gov/products/mod13a2v006, accessed on 12 June 2022), which provides Normalized Difference Vegetation Index (NDVI) values every 16 days at a spatial resolution of 1 km. The NDVI was used to determine trends in values over a time series in forested and deforested areas. The NDVI shows values between −1 and 1 and is determined using the near-infrared (NIR) and red (RED) bands, where values near −1 indicate little photosynthetic activity and, thus, little vegetation growth or reduction, while values near 1 reflect the opposite [28]. MOD13A2 is designed to provide spatially and temporally consistent comparisons of vegetation conditions and was processed from the MODIS L3 daily surface reflectance product (MOD09), corrected for the effects of atmospheric gasses, thin cirrus clouds, and aerosols [29]. In recent decades, MOD13A2 data have been used to quantify vegetation activity and to detect vegetation dynamics in many biological communities [30][31][32].
Land Surface Temperature. Land Surface Temperature (LST) is one of the most important parameters in the physical processes of surface energy and water budgets from local to global scales [8,33,34]. Remotely sensed thermal infrared data provide spatially continuous LST measurements with global coverage to examine the thermal heterogeneity of the Earth's surface and the impact of natural and human-induced changes on surface temperatures [3]. LST data were selected from the daily product MYD11A1 ver. 6 (https://lpdaac.usgs.gov/products/mod11a2v006, accessed on 12 June 2022), obtained during the daytime and nighttime series from a MODIS sensor with a 1 km spatial resolution for the period 2001 to 2020. MODIS land surface temperature and emissivity products provide per-pixel emissivity and temperature values in a sequence from global swath-based products to grid-based products [35]. In recent years, studies using MODIS data were conducted to analyze the temporal and spatial patterns of regional surface temperature; good results were achieved [36]. In addition, daytime and nighttime LST products have already been validated [37]. The original MODIS LST product was obtained in Kelvin ( • K) values, which were then converted to Celsius ( • C).
Terrestrial Coverage Type. We used the MODIS product MCD12Q1.006 (https:// lpdaac.usgs.gov/products/mcd12q1v006, accessed on 1 December 2021), which provides a Land Cover Type (TCT) classification developed by the University of Maryland (Type 2). The TCT consists of a series of coverage maps with annual interventions and a spatial resolution of 500 m [38] available for the period 2001 to 2020. This product was used to define the selection of sampling sites.
Meteorological stations. We used the monthly average data of maximum and minimum temperatures (in • C) provided by the Bolivian Meteorological Service (SENAMHI, https://senamhi.gob.bo, accessed 10 November 2022) for the urban areas of Concepción and San Ignacio for the period 2001-2020.

Sampling Sites and Data Processing
Using a stratified randomized design, ten sampling points were selected, with five in forested areas and five in areas with deforestation. Then, 10 km rectangular buffers were created around these central points to obtain 10 × 10 km quadrants representing an area of 100 km 2 for each sampling site. The criteria for selecting these sites were based on a distance of more than 100 km from each other. The main criterion for forested areas was that the areas were legally protected or were indigenous territories larger than 100 km 2 . We selected five quadrants with intact forests in 2020, with three located in northern Chiquitania and two distributed in southern Chiquitania. The main selection criterion for deforested areas was that they represent important urban areas or show different levels of land-use change between 2001 and 2020 due to productive agricultural or livestock activities, identified based on the product MCD12Q1.006. Deforested areas were represented by quadrants located in two population centers in northern Chiquitania and areas with high deforestation rates in southern Chiquitania ( Figure 2). Table 1  (quadrants) in terms of their names, forest cover, geographic location, and climate. The climate region was identified based on a global map of the aridity index [39].

Sampling Sites and Data Processing
Trends of change. Differences in trends in MODIS product time series averages of daytime and nighttime temperatures (LST) and vegetation cover (NDVI) were analyzed using a nonparametric Mann-Kendall statistical test [43,44]. This test is commonly used as a trend test for the median slope operator and produces Z-score outputs, allowing both the significance and trend direction to be assessed. The trend values were obtained in RStudio using the Trend package [45].
Correlation between variables. To find the significance (p < 0.05) and degree of the relationship between variables, the intersection between NDVI values and daytime and nighttime LST data was determined using GEE and Pearson's correlation coefficient (r). Pearson's coefficient is defined as the association and relationship between two variables and the strength of their relationship. In other words, this coefficient shows how change in one variable causes a change in the other variable, with values ranging between −1 and 1, where a value of −1 indicates a total negative linear relationship, 0 indicates the absence of correlation, and +1 indicates a total positive correlation.
Technical validation. We calculated the monthly averages of the daytime and nighttime LST data and performed an accuracy assessment with the information obtained from the weather stations (maximum and minimum temperatures) based on Pearson's coefficient and the Root Mean Square Error (RMSE).  The analyses were performed using the Google Earth Engine (GEE), a cloud computing platform designed to store and process large datasets (petabyte scale) for decision making [40,41]. We used the 10 polygons to extract the land cover product values MCD12Q1.006 in ArcMap and reclassified the 15 class types into 7 categories: water bodies, forests, shrublands, savannas, grasslands, urban/croplands, and barren. We then calculated annual percentages for each of the categories (2001 to 2020). To assess the level of uncertainty in the resulting classification, 100 random verification points were distributed in ArcMap, grouped into the 7 categories. In GEE, we selected Sentinel-2 satellite imagery to verify the classification. This approach allowed us to obtain a confusion matrix [42] and three types of accuracy estimates-overall accuracy, user accuracy, and producer accuracy-including their 95% confidence intervals. The confidence level obtained from the coverage classification was 90%.
In the GEE platform, the ten polygons (10 × 10 km) were used and a script was developed for analysis purposes using information from NDVI and LST products for the period from 2001 to 2021. Each quadrant was composed of 100 pixels of MODIS products (1 pixel = 1 km 2 ), so the number of samples in quadrants was calculated with a heterogeneity of 50% and a confidence level of 95% through random point scattering in GEE. Following the methodology proposed by Olofsson et al. [42], 86 points were distributed in each quadrant, resulting in a total of 860 sampling points.

Statistical Analysis
Trends of change. Differences in trends in MODIS product time series averages of daytime and nighttime temperatures (LST) and vegetation cover (NDVI) were analyzed using a nonparametric Mann-Kendall statistical test [43,44]. This test is commonly used as a trend test for the median slope operator and produces Z-score outputs, allowing both the significance and trend direction to be assessed. The trend values were obtained in RStudio using the Trend package [45].
Correlation between variables. To find the significance (p < 0.05) and degree of the relationship between variables, the intersection between NDVI values and daytime and nighttime LST data was determined using GEE and Pearson's correlation coefficient (r). Pearson's coefficient is defined as the association and relationship between two variables and the strength of their relationship. In other words, this coefficient shows how change in one variable causes a change in the other variable, with values ranging between −1 and 1, where a value of −1 indicates a total negative linear relationship, 0 indicates the absence of correlation, and +1 indicates a total positive correlation.
Technical validation. We calculated the monthly averages of the daytime and nighttime LST data and performed an accuracy assessment with the information obtained from the weather stations (maximum and minimum temperatures) based on Pearson's coefficient and the Root Mean Square Error (RMSE).

Temporal Dynamics of Land Cover Type
Our results show that the five sites with forest did not feature variations in the percentages over the 20 years of analysis ( Figure 3). However, the urban areas of San Ignacio and Concepción showed a reduction in forest areas and savannas but an increase in grasslands. In the case of El Cerro/California, the forested area was completely reduced, while grasslands and croplands increased. The area of Nueva Esperanza was the most strongly impacted and the area of Santa Ana/Buena Vista reduced its forest cover from 100 to 61% ( Figure 3). centages over the 20 years of analysis ( Figure 3). However, the urban areas of San Ignacio and Concepción showed a reduction in forest areas and savannas but an increase in grasslands. In the case of El Cerro/California, the forested area was completely reduced, while grasslands and croplands increased. The area of Nueva Esperanza was the most strongly impacted and the area of Santa Ana/Buena Vista reduced its forest cover from 100 to 61% ( Figure 3).

Annual Trends of Forest Cover Loss in the Chiquitania Region
The NDVI time series for each study site (Table 2) is shown in Figure 4. There was a small but significant trend (p < 0.05) in the mean values of the index in areas with no forest

Annual Trends of Forest Cover Loss in the Chiquitania Region
The NDVI time series for each study site (Table 2) is shown in Figure 4. There was a small but significant trend (p < 0.05) in the mean values of the index in areas with no forest change (Noel Kempff, Monteverde, San Rafael, and Tucabaca). In sites with increasing deforestation rates during the analysis period, NDVI showed a negative trend. The strongest trend occurred in the El Cerro/California Mennonite colonies (p < 0.001), where all forest cover in 2001 (12% of the quadrant area) was lost by 2021. The next strongest trend was observed at the Santa Ana/Buena Vista site (p < 0.001), where 61% of the original forest area remained by the end of the study period. Significant negative trends were also recorded in San Ignacio, although the area lost 5% of its natural vegetation cover (p < 0.05).    The main variations were associated with deforestation sampling sites ( Figure 5). The areas of the population centers and surroundings of San Ignacio (28.7 ± 0.56 °C) and Concepción (28.7 ± 0.74 °C) showed lower average temperatures compared with the areas of the Santa Ana/Buena Vista ranches (30.5 ± 3.95 °C). Even lower temperatures were observed in the areas of the Mennonite colonies of El Cerro/California (32.3 ± 1.74 °C) and  The main variations were associated with deforestation sampling sites ( Figure 5). The areas of the population centers and surroundings of San Ignacio (28.7 ± 0.56 • C) and Concepción (28.7 ± 0.74 • C) showed lower average temperatures compared with the areas of the Santa Ana/Buena Vista ranches (30.5 ± 3.95 • C). Even lower temperatures were observed in the areas of the Mennonite colonies of El Cerro/California (32.3 ± 1.74 • C) and Nueva Esperanza (33.4 ± 1.34 • C). In the latter, the forest cover was completely lost. Between deforested areas, daytime LST averages showed a difference of up to 4.7 • C. Comparatively, daytime LST averages in deforested areas were warmer than those in forested areas, with a difference of 3.1 • C.
In addition, in the deforested areas, a pattern of temperature increase anomalies in the daytime LST averages was evident for 2002 and 2010-2012 and was even higher for 2019-2020 ( Figure 5), especially in the areas of the Mennonite colonies El Cerro/California and Nueva Esperanza. In forested areas, this pattern was also evident in southern Chiquitania (San Rafael and Tucabaca).

Annual Trends of Daytime LST
Regarding the trends of the increase in daytime LST ( • C), the results showed a positive increase at four of the five sites with the loss of forest cover (Table 3, Figure 6). The highest trend was recorded in the San Ignacio area (Z = 45.10; p > 0.001), followed by the El Cerro/California Mennonite colonies (Z = 8.35; p > 0.001), Santa Ana/Buena Vista ranches (Z = 5.75; p > 0.001), and Concepción town center (Z = 45.10; p > 0.05). For sampling sites with intact forests, no significant changes in diurnal temperature trends were observed (Table 3 and Figure 6).   ranches (Z = 5.75; p > 0.001), and Concepción town center (Z = 45.10; p > 0.05). For sampling sites with intact forests, no significant changes in diurnal temperature trends were observed (Table 3 and Figure 6).

Relationship between Intact Forest and Areas with Deforestation Processes and Daytime and Nighttime LST ( • C)
The results also indicated that there was a statistically significant medium to high correlation between NDVI and daytime LST, especially for the Santa Ana/Buena Vista ranches (r = 0.76), San Ignacio (r = 0.60), and El Cerro/California Mennonite colonies (r = 0.60) (Table 5, Figure 9). As the NDVI levels decreased, the temperature increased. For nighttime LST, a significant relationship was found between the loss of forest cover and the increase in nighttime temperature in the Santa Ana/Buena Vista (r = 0.59) and El Cerro/California (r = 0.58) areas (Table 5, Figure 10). These results suggest that deforestation increased the land surface temperature in the Chiquitania region.  Figure 8). The rest of the sites, despite indicating an increase and decrease, were not statistically significant.      Table 6 shows the validation statistics (performance metrics) of the monthly average daytime LST and nighttime LST in relation to air temperature data from weather stations in the urban areas of Concepción and San Ignacio. The results of these correlations were also significant and very high for the two study sites. However, the RMSE values were lower for Concepción, with 5.4 °C for daytime LST and 5.5 °C for nighttime LST.  Table 6 shows the validation statistics (performance metrics) of the monthly average daytime LST and nighttime LST in relation to air temperature data from weather stations in the urban areas of Concepción and San Ignacio. The results of these correlations were also significant and very high for the two study sites. However, the RMSE values were lower for Concepción, with 5.4 • C for daytime LST and 5.5 • C for nighttime LST. Table 6. Pearson correlation (r) and Root Mean Square Error (RMSE) between the monthly averages for LST and the air temperatures of the meteorological stations (Concepción and San Ignacio). The highest correlations are highlighted in bold.

Maximum
Temperature

Discussion
Deforestation in Bolivia is growing at an alarming rate. Recent estimates indicate that the country lost a total of seven million hectares by 2021, 86% of which will be concentrated in the department of Santa Cruz [19]. In addition, the advance of deforestation has been associated in recent years with areas destined for forest use and exploitation [20], as well as legally protected areas [20]. The expansion of soybean and cattle ranching has been the main direct cause of deforestation [19]. Our analyses sought to identify changes in forest cover in sites located in protected areas and indigenous territories using the MODIS product MCD12Q1.006, as well as changes in areas known to be the fastest growing areas of deforestation in Chiquitania. While the area of Nueva Esperanza presented the highest concentration of cropland, we found an evident reduction of natural cover in urban areas (Concepción and San Ignacio), as well as in the areas of El Cerro/California and Santa Ana/Buena Vista. If this deforestation trend continues, it is estimated that forest could disappear completely by 2050 [20].
Assessments of vegetation cover status, changes, and processes are important components of global change research programs and are topics of considerable societal importance [46]. Spectral vegetation indices (e.g., NDVI) are among the most widely used satellite data products and provide metrics for climate, hydrological, and biogeochemical studies, land cover change detection, phenology, and natural resource management [47]. In addition, these indices are important for the analysis of vegetation growth and decrease trends [48]. Our results were mainly based on a comparison of vegetation trends in forested and deforested areas. Using NDVI, we detected continuous vegetation growth in intact forests (Noel Kempff, Monteverde, and San Rafael), which showed no evidence of forest degradation processes in these areas. The explanation for these results is that high NDVI values are mainly related to the density of green leaves in a given area, making them a good indicator of vegetation cover and vitality [49]. Our assessments considered a long time series (2001-2020) to ensure that the results would not be influenced by seasonal and interannual variability in the forests of the Chiquitania region [50]. In contrast, we found sites with negative and statistically significant trends in NDVI values in three of the five sampling sites (El Cerro/California Mennonite colonies, Santa Ana/Buena Vista ranches, and San Ignacio), where a decrease in forest cover was recorded at different sites in Chiquitania. This result was expected, especially considering the current scenario of accelerated land-use change processes in the region [19]. In addition, there was evidence of a negative trend marked by decreases of the NDVI in the Tucabaca forests. These results are very close to those of San Rafael, a forested site with a positive trend, which suggests disturbances due to illegal logging; however, ecological disturbance processes (e.g., changes in stomatal conductance) that cause physiological weakening and tree mortality could not be ruled out. However, the results for Tucabaca should be interpreted with caution and field evaluations should be carried out to identify the factors causing this degradation.
Land Surface Temperature (LST) monitoring using remote sensors offers several advantages over a wide observation range, including easy access and strong spatial continuity.
Multiple studies have been conducted globally in recent years to determine LST [3,8,51,52], which is a key parameter for obtaining a more direct measure of surface conditions when analyzing surface-climate interactions [33,52]. Due to these qualities, we were able to determine, over a long period of time (2001-2020), the diurnal and nocturnal temperatures at the local level in Chiquitania. Of the 10 study sites, the highest diurnal LST levels were found in Nueva Esperanza, a completely deforested area. In addition, in the Cerro/California Mennonite colony area, which lost its remaining forest cover, a notable increase in diurnal LST was identified from 2007 onwards, which coincides with the decline in the trend of photosynthetic activity of the vegetation presented in Figure 3. In addition, we were able to compare the differences in LST between forested and non-forested areas. Our results showed that the average values of daytime LST between forested sampling sites varied by up to 2.7 • C and those between deforested sites varied by up to 4.7 • C. However, when comparing forested and deforested sites, the difference in the averages was 3.1 • C. Differences in the diurnal temperature range are known to be smaller for forested areas and larger for non-forested areas [53], because LST increases in deforested areas may be exacerbated by changes in roughness length, which impede energy dissipation via sensible or latent heat fluxes [54]. In contrast, nocturnal heating is caused by the release of stored thermal energy during the day [3]. However, comparisons of nighttime LST averages for the ten sampling sites did not show notable differences and, for this reason, nighttime LST could not be considered a good indicator of the differences between forested and deforested areas in Chiquitania. Nevertheless, our data help highlight the large differences between daytime and nighttime LST averages, which for the forest was 7.1 • C and for the deforested areas was 10.2 • C. This information is valuable, because it highlights the importance of forests in regulating temperatures in the Chiquitania region.
It is known that anomalies of maximum LST capture spatial patterns that are associated with droughts and heat waves on the land surface [8]. In our study, both forested areas and those with deforestation processes in Chiquitania showed interannual anomalies of temperature increases in the daytime LST averages for the years 2002, 2010-2012, and 2019-2020, with the latter increases being higher than the others. In addition, in terms of nighttime LST, similar patterns were recorded for 2002 and 2015. In the Amazon area, generalized anomalies in LST values and large-scale directional changes towards higher temperatures were found in 2005 and 2010 [8], where large and severe droughts were recorded [55]. However, by 2020 different severe to extreme meteorological megadroughts events were recorded at the continental level. One of these events occurred in southern Brazil and Paraguay [56,57] and was also observed in Bolivia. In addition, on a local scale, drought events in the Chiquitania region were identified based on different climate stations, showing significant trends of increased frequency and intensity [21]. Further studies are needed to understand the relationship between heat waves and droughts, which could improve our understanding of the impacts of land-use change on the local climate.
Recent studies indicate that in some lowland areas of Bolivia, temperatures tend to increase mainly in urban and agricultural/livestock areas [19]. In our research on Chiquitania, it was shown that there is a tendency for an increase in daytime LST values, which is associated with deforestation processes, independent of the differences between the climate regions (humid and dry sub-humid areas) of the sampling sites. The trend of an increase in diurnal LST was statistically significant in four of the five areas with deforestation, mainly in San Ignacio, which registered the highest value in the Mann-Kendall test. However, in the sampling sites with forests, there were no statistically significant trends in daytime LST values, because tropical forests can exert a strong influence on local LST. As such, forests generally have a lower surface albedo due to absorbing more shortwave radiation during the day and higher evapotranspiration compared with areas with bare vegetated surfaces [3,58]. Nevertheless, nighttime LST values corresponded to another scenario, with a significant increase observed for the two forested areas (Noel Kempff and Bajo Paraguá) and two urban areas (San Ignacio and Concepción). Although an increasing trend in areas with deforestation was expected, new questions emerged regarding nighttime LST trends in forests. A possible explanation for this phenomenon is the presence of a large number of granitic rocky outcrops (inselbergs) in the Noel Kempff and Bajo Paraguá areas [54], which concentrate the heat during the day and radiate heat to the surrounding forest during the night [59,60]. However, further research is needed in this area.
Globally, deforestation in tropical regions is causing strong warming between 0.38 ± 0.02 [52] and 2.4 ± 0.10 • C [3] (LST values). Our results showed a statistically significant correlation between NDVI and LST at three sites for diurnal values and two sites for nocturnal values. This result is in addition to the trend of 0.1 • C per decade observed for the air temperature increase rate [61]. This result is concerning because the Chiquitania region is home to hundreds of thousands of people who carry out productive livestock and agricultural activities. The increase in LST values in the Chiquitania region could lead to a series of effects that have not been evaluated, including migration, economic production, and human capital to impacts on biodiversity [62]. Given that forest loss is at the heart of both global and local warming, initiatives to reduce deforestation must remain a priority [63]. Climate change adaptation strategies based on maintaining forest integrity are required to mitigate the increase in LST in the Chiquitania region [64], where the conservation of existing forests and restoration of those lost to deforestation should be prioritized [65].
Protected areas continue to be the main defense against forest cover loss and the best strategy for maintaining the ecological integrity of forests [66]. Bolivia has approximately 130 protected areas, of which approximately one-third are located in the Department of Santa Cruz [67]. Most of these protected sites contain ecosystems of high conservational value and, in many cases, intact forests [68]. However, there is direct and indirect pressure around these areas owing to the rapid expansion of the agricultural and livestock frontier [19,69]. In the Chiquitania region, this expansion has caused a loss of forest connectivity [19], a reduction of critical ecosystems [70,71], and a loss of habitats of key and priority species for conservation [72]. If this scenario continues in these protected areas, up to half of forest area harbors could be lost in the next 30 years [68].
Reducing current deforestation trends in Chiquitania is a high priority. It is necessary to redefine norms and public policies regarding forest conservation, aimed at territorial management at multiple scales. Undoubtedly, revising the design criteria for land-use plans, enforcing the economic and social functions of private and communal property, and promoting the protection and restoration of conservation easements will improve opportunities for thermal regulation and the maintenance of ecological functionality at the landscape scale in the Chiquitania region. It is also necessary to establish a monitoring system to evaluate changes in natural vegetation cover, especially temperature and precipitation, at regional and local scales. With these changes, there will be greater opportunities to increase socio-ecological resilience to the impacts of climate change.

Conclusions
This is the first study conducted in the Chiquitania region of Bolivia over a long time series to identify the relationship between land surface temperature increases and forest cover loss through remote sensing products. NDVI remains an essential tool for monitoring the health of the field and determining the loss or maintenance of forest cover, especially in the present scenario where deforestation is advancing rapidly. In addition, LST was shown to be a key parameter that can help estimate temperature variations at different scales of space and time. Our findings demonstrate the importance of using combined time series data of spectral indices and LST to fully understand temperature trends and their drivers in the Bolivian lowlands.
These results indicate that the forests of the Chiquitania region play a fundamental role in the regulation of temperature. Our study highlights the importance of protected areas in the Chiquitania region for the conservation of forests, which, in turn, affects the regulation of land surface temperatures. The evaluated forest sites represent national (Noel Kempff), departmental (Tucabaca), and municipal (Bajo Paraguá and San Rafael) protected areas, as well as an indigenous territory (Monteverde), demonstrating that the legal allocation of land tenure and territorial management focusing on the conservation and sustainable use of wild resources are important for maintaining the integrity and conservation of natural ecosystems. However, this scenario may change in the face of land-use change trends. In this sense, conservation strategies and actions are very important for safeguarding the maintenance of forest integrity, biodiversity, and environmental functions. For this reason, it is necessary to continue working on management and capacity building to protect these sites.