Hkakabo Razi landscape as one of the last exemplar of large contiguous forests

Deforestation and forest degradation around the world endanger the functioning of ecosystems, climate stability, and conservation of biodiversity. We assessed the spatial and temporal dynamics of forest cover in Myanmar’s Hkakabo Razi Landscape (HRL) to determine its integrity based on forest change and fragmentation patterns from 1989 to 2016. Over 80% of the HRL was covered by natural areas, from which forest was the most prevalent (around 60%). Between 1989 and 2016, forest cover declined at an annual rate of 0.225%. Forest degradation occurred mainly around the larger plains of Putao and Naung Mung, areas with relatively high human activity. Although the rate of forest interior loss was approximately 2 to 3 times larger than the rate of total forest loss, forest interior was prevalent with little fragmentation. Physical and environmental variables were the main predictors of either remaining in the current land-cover class or transitioning to another class, although remaining in the current land cover was more likely than land conversion. The forests of the HRL have experienced low human impact and still constitute large tracts of contiguous forest interior. To ensure the protection of these large tracts of forest, sustainable forest policies and management should be implemented.

Despite human practices and the unprecedented use of natural resources, forests are still widely distributed globally and cover around 30% of the Earth's surface 1,2 . However, ongoing deforestation and forest degradation jeopardize the functioning of biogeochemical and hydrological cycles, climate stability and conservation of biodiversity [3][4][5] . Net loss of forest area occurs largely in the tropics 6,7 and this forest loss continues to impact areas with particularly high conservation value 8,9 .
Tropical forests play a key role in the global carbon cycle and support more than half of the world's biodiversity 10 . Industrial logging, agricultural expansion, fire, mining/resource extraction and urban growth have led to extraordinary loss of tropical forest 11,12 . The amount of forest loss differs between continents 13 , with the highest levels occurring in South America and Asia 2 . In Southeast Asia, Myanmar had the second highest rate of net forest loss between 1990 and 2015, trailing only Indonesia, with a loss rate of 546,000 ha y −1 between 2010 and 2015 1 . Furthermore, this rate of forest loss represented a 25% increase since the 1990s. The driving forces behind the high rates of forest loss in Southeast Asia are logging and the global demand for crops such as oil palm, sugar, and wood fibre 14 .
Despite having the third largest annual forest loss in the world between 2010 and 2015 1 , Myanmar remains one of the most heavily forested countries in Southeast Asia 15,16 . Myanmar is the second largest exporter of Teak (Tectona grandis), a valuable timber species, and much of the rural population continues to depend on forests to supplement their livelihoods 17 . Some forest areas are used for small scale agroforestry and up to 77% of energy demands are covered by traditional energy sources such as fuel wood, charcoal and biomass 18 . Selective logging on government forest reserves has historically been managed under the Myanmar Forest Selection System, which sets harvest quotas to sustain long-term timber yields 19 . In unmanaged forests, though, logging concessions have far less oversight 20 and contribute to the rapid loss of relatively intact forest. Besides wood extraction, agricultural Scientific RepoRtS | (2020) 10:14005 | https://doi.org/10.1038/s41598-020-70917-y www.nature.com/scientificreports/ expansion and infrastructure development are the most common causes of forest loss 21 . Nonetheless, Myanmar has retained much of its original forest cover, stretching across 63% of the country's land 16 .
One of Myanmar's largest remaining area of contiguous intact forest is found in the northern-most part of the country 16 . The Hkakabo Razi Landscape (HRL; also known as the Northern Forest Complex or Northern Mountain Forest Complex) has an exceptionally rich biodiversity and some of the highest concentrations of endemic species in the world 22,23 . New species are continuously being described, including plants [24][25][26][27][28] , mammals 29-31 and birds 32-35 , described to be new for the HRL 36,37 or reconfirmed to occur for the HRL 38 . In 2000, human settlements and paddyfields covered less than 14% of the area 39 , and changes in forest cover from 1990 to 2000 were minimal with an overall deforestation rate of < 1% 39 . Access to the area is limited, with most areas reachable only by walking trails as roads only exist close to towns in the southern part of the HRL 39 .
However, the remaining large forest tracts in Myanmar are under pressure 40 . The country's recent political and economic reforms are attracting investors, leading to far-reaching changes in forestry and other land-based investments sectors 16,[41][42][43] . Therefore, areas previously inaccessible are starting to open up for systematic, largescale resource extraction, timber production, commercial plantations (MOECAF 2011) and extensive tourism 44 as a way to achieve economic growth. Given these rapid developments, assessing the current forest distribution and quantifying the extent of forest loss and degradation in the HRL are required to understand the effects economic and development changes are having on this region, and to develop strategies for sustainable land management.
The objective of this study was to assess the spatial and temporal dynamics of forest cover in the HRL. We mapped land cover, including several ecologically-distinct forest types, using multi-temporal Landsat satellite imagery and quantified forest change and fragmentation patterns across a 27-year interval. We then identified major spatial predictors of forest change. Through this assessment, we aim to determine whether the forests of the HRL are still of high integrity 45 .

Methods
Study area. Our study area was located in the northernmost tip of Kachin State, Myanmar (27°46′50.55″ N, 97°41′11.27″ E). It is surrounded by Yunnan/China to the east, Arunachal Pradesh/India to the west and the Tibetan Plateau to the north (Fig. 1). The HRL includes the Hkakabo Razi National Park (established in 1998), Hponkan Razi Wildlife Sanctuary (established in 2003), and the currently proposed, but not formally established, southern extension of the Hkakabo Razi National Park. The combined designated and proposed protected areas cover 11,820 km 2 . We added a buffer of 15 km around this landscape and included the Putao plains as part of our study area (Fig. 1). Land-cover mapping. We performed two single-date supervised land-cover classifications for the years 1989 and 2016 to map recent changes in 15 land-cover classes. These classifications were based on two adjacent Landsat 5 TM (1989) and two adjacent Landsat 8 OLI images (2016). The two Landsat images in each year had the same acquisition date and were captured during the dry season to minimize cloud cover interference ( Table 1). All images were pre-processed by converting to Top of Atmosphere (TOA) reflectance and masking out clouds and cloud shadows using the FMASK algorithm 54 . We also used the C correction method to lessen the effects of topographic shadowing through topographic normalization 55 .
Training data for supervised land-cover classification were created by a combination of ground truth data and visual interpretation of fine-resolution imagery freely available through Google Earth. Ground truth data were collected by MSR, MK, and SCR. Visual interpretation of fine-resolution imagery was based on patch size and geometric shape, texture, colour, or vegetation phenology in cases where time series of reference imagery were available. We delineated separate sets of training and validation polygons for each of our 15 landcover classes (Table 1). Training and validation data for forest classes in 2016 were largely identified based on ground truth data and prior knowledge of particular forest areas (field expeditions in December 2013/January 2014, January-March 2016). We created separate training data for forests in three elevation ranges (< 600 m, 600-1,800 m and > 1,800 m) due to expected differences in plant communities and associated spectral characteristics. Training and validation data for secondary (degraded) forest classes represented areas where early dry season canopy cover appeared to be less than 80%. Because intact forest canopies in the HRL are mostly closed 16 , we considered reduced canopy cover to be an indication of degradation. Importantly, the estimated extent of degraded forest, therefore, does not capture forest degradation that does not result in a reduction of canopy cover. Training and validation data for lowland forest were created only in areas less than 600 m in elevation that were visually identified as having flat or gently undulating terrain. The final training dataset consisted of approximately 70 polygons per land-cover class for a total of 1,080 polygons for 1989 and 1,047 for 2016. The separate validation dataset consisted of at least 20 additional polygons per land-cover class.
We also compiled several auxiliary data layers as predictor variables in our land-cover classification, because we expected different plant communities on north-vs. south-facing slopes, ridges vs. valleys, and steep slopes vs. flat flood plains. These included a 30-m digital elevation model (DEM) from NASA's Shuttle Radar Topography Mission 56 , and several topographic measures derived from this DEM: slope, northness 57 , and topographic position index (TPI) 58 . As an auxiliary data source for distinguishing clearcuts from other bare ground (e.g., agriculture), we also calculated the change in Normalized Differenced Vegetation Index (NDVI) 59,60 between 2013 and 2016. This data layer was created in Google Earth Engine 61 , a cloud-based analysis platform, by calculating We conducted a supervised classification using a Random Forest algorithm with the randomForest package 62,63 in the R environment 64 . Random Forest is a non-parametric, machine learning classifier commonly used in remote sensing for classification of satellite or aerial images [65][66][67] . It is flexible concerning distributional assumptions about the training data, generally robust to over-fitting 68 , and thus, potentially better suited than non-parametric approaches for identifying target land-cover classes that include mixtures of spectral information 65,69 . To train the Random Forest classifier, we randomly selected 30 pixels from each training polygon for each land-cover class and used elevation, slope, northness, topographic position, single-band Landsat data (Landsat 8 OLI bands 2-7; Landsat 5 TM bands 2-6), and NDVI change (for 2016 only) as predictor variables. We set the following parameters for the Random Forest classifier: 500 decision trees, 2/3 of the training data sampled with replacement as the training set for each decision tree, and the number of predictor variables randomly sampled at each node split set to the square root of the number of variables. Due to the presence of cloud and cloud shadow gaps in our 2016 land-cover classification (up to 17% of Landsat Scene with clouds; Table S1), particularly at high elevations, we performed a secondary radar-based classification to fill gaps in our primary Landsat-based classification. Predictor variables for this secondary classification were Sentinel-1A VV and VH radar backscatter data and their ratio, collected in the Interferometric Wide (IW) swath acquisition mode. Sentinel-1A data were pre-processed by performing radiometric calibration, speckle noise filtering, Range-Doppler terrain correction, and conversion to decibel scaling using the SNAP 2.0 software 70 . Following gap-filling of the primary 2016 classification, we post-processed both the 1989 and 2016 classified images. We re-assigned forest pixels to the land-cover class corresponding with their actual elevation band if needed. Afterwards, we used a 3 × 3 majority filter to smooth isolated pixels that likely represent classification error.
We validated the accuracy of both single-date land-cover classifications (1989 and 2016) by randomly selecting 10 reference pixels from each of the 20 + validation polygons for each land-cover class in each year and comparing to the final classifications. Based on these reference points, we calculated standard accuracy metrics including an error matrix, producer's accuracies (fraction of reference pixels for a given class that are correctly identified), user's accuracies (fraction of pixels of a given class that are correct in the classified image), and overall accuracy following the recommendations of Olofsson et al. 71 .
Landscape change analysis. As we were primarily interested in information on forest integrity to support conservation activities, we grouped land-cover classes into three main vegetation categories for subsequent analysis: (i) alpine vegetation; (ii) forest, which included forest at different elevation bands and fir/rhododendron; and (iii) shrubland, which comprised secondary forest and shrubs. Grassland/pasture, paddyfield, settlement and clearcut were also merged into "agriculture/developed". We used these land-cover maps for constructing a transition matrix to summarize the area and percentage of land-cover change from 1989 to 2016 based on a pixel-by-pixel comparison 72 . To complement the transition matrix, we plotted the change values in a Sankey diagram (available at https ://githu b.com/csala denes /sanke y, accessed 8 May 2019). Sankey diagrams, traditionally applied in the analysis of flows of energy and materials, have been recently employed for visualizing land-cover change transitions 73 .
In addition, we quantified forest fragmentation by computing the spatial density of forest cover, the Forest Area Density (FAD) 74,75 . FAD is defined as the proportion of the pixels in a surrounding fixed-area neighbourhood that is forest and measures both the area of continuous forest and of patches of forest separated by nonforest lands 76 . Each forest pixel is labelled as forest interior if the associated FAD ≥ 0.9 77 . To account for the scaledependence of fragmentation, we evaluated FAD at five measurement scales defined by neighbourhood sizes equal to 4.41 ha (7 pixels × 7 pixels), 15

Spatial predictors of forest change.
To understand factors that may explain the spatial dynamics of forest change, a multinomial logistic regression was performed 79,80 . We evaluated the effects of physical (slope, elevation, major landform, soil properties), environmental (temperature, precipitation, soil degradation) and geographical (distance from villages, from towns, from roads, from rivers) explanatory variables on the risk of forest transition between 1989 and 2016. Slope was derived from the digital elevation model (DEM) from NASA's Shuttle Radar Topography Mission 56 . Soil properties were obtained from the World Reference Base (WRB), an international system for classification of soils 81 . Soil degradation and landform were based on a map produced by the UNEP Global Assessment of Human-Induced Soil Degradation (GLASOD). Soil degradation was identified based on expert opinion as all areas where human intervention had resulted in a decrease on the capacity of the soil to support life 82 . Mean temperature and annual precipitation for our analysis were based on the Global Climate Data WorldClim 2 83 . Euclidian distance from villages, towns and roads were calculated based on a national database from the MIMU-Myanmar Information Management Unit (https ://themi mu.info/ accessed 3 June 2019) and river data were accessed through the HydroSHEDS database (https ://hydro sheds .cr. usgs.gov/hydro .php accessed 3 June 2019).
We used multinomial logistic regression models to predict transition probabilities between possible landcover categories. In our case, there were three possible discrete outcomes for a pixel belonging to a given Scientific RepoRtS | (2020) 10:14005 | https://doi.org/10.1038/s41598-020-70917-y www.nature.com/scientificreports/ land-cover category (e.g., for forests: forest to forest; forest to shrubland; forest to agriculture/developed). We fit three separate multinomial logistic regression models, one for each land-cover category (forest, shrubland, agriculture/developed), and estimated the changes in transition probabilities due to slope, elevation, landform, soil properties, temperature, precipitation, soil degradation, and distances from villages, towns, roads, and rivers. In the first model, possible outcomes of the dependent variable were forest to forest, forest to shrubland and forest to agriculture/developed (n = 16,058,849). In the second model, possible outcomes were shrubland to shrubland, shrubland to forest and shrubland to agriculture/developed (n = 587,117). Lastly, in the third model possible outcomes were agriculture/developed to agriculture/developed, agriculture/developed to forest, and agriculture/developed to shrubland (n = 630,992). To deal with spatial autocorrelation, each of the models was run using a randomized subsample of 10% of the entire data set. This procedure was repeated 1,000 times and the mean of each coefficient estimate and 95% confidence intervals were calculated. The multinomial models were run in the nnet R package 84 . To aid the interpretation of the estimates of the multinomial logistic regression models (i.e., log of the odds ratio), we calculated the relative risk ratio (RRR) as exp(β) 85 . RRR > 1 indicates that the land-cover transition is more likely and RRR < 1 indicates that the land-cover transition is less likely (i.e., remaining in the current land cover is more probable).

Results
Land-cover mapping. Two land-cover maps, each one for 1989 and 2016, were produced ( Fig. 2) with overall accuracies of 90% and 85%, respectively ( Table 2). The mean per-class user's accuracy was 94% for 1989 and 91% for 2016. Forest at 600-1,800 m was occasionally misclassified in both years but more frequently in 2016. For the 2016 land-cover map, 30% of misclassified forest at 600-1,800 m was due to confusion with secondary forest, shrub/bush and clearcut. Shrub/bush was mistakenly classified in 10% of the cases as secondary forest at 600-1,800 m, and secondary forest at low elevation was confused in 14.5% of the cases with shrub/bush. Most classes were mapped with high accuracy in 1989, but for 2016 four classes were misidentified. Although clearcut was correctly identified in most cases, 34% of clearcut was misidentified. Similarly, secondary forest at low and mid-elevation and shrub/bush were misidentified in comparable magnitude ( Table 2). Based on the error-adjusted estimates, forest covered most of the total area in 1989 (63.6%; Fig. 2, Table 3). Forest at 600-1,800 m was most common (24.5%), followed by forest > 1,800 m (17.5%) and fir/rhododendron (16.1%). Forest at low elevation (< 600 m) covered only 5.6% of the area. The extent of forest degradation, or areas showing canopy damage (i.e. secondary forest and shrub/bush), was low (5.4%). Grassland/pasture, paddyfield, clearcut and settlement (3.9%) occurred mainly at lower elevations in the Putao plains and around Naung Mung village. In 2016, forest also covered most of the area (57.8%). Secondary forest and shrub/bush covered 11.4%, mainly in the north and west of the Mali Raing region, west of the Putao plains, and east of Naung Mung. Clearcut areas covered 3.3%, whereas grassland/pasture, paddyfield and settlement covered almost the same area as in 1989.
Landscape change from 1989 to 2016. The analysis of 27 years of land-cover change indicated that the deforestation rate for the study period was 2.7%. Forest cover decreased from 1,375 thousand ha in 1989 to 1,294 thousand ha in 2016 (annual average rate = 0.225%). The loss of forest was mainly due to an increase in shrub/ bush and secondary forest (5.9%) from 84 thousand ha in 1989 to 127 thousand ha in 2016 (Fig. 3). Agriculture/ developed areas (i.e. grassland/pasture, paddyfield, clearcut and settlement) increased 1.7%, mainly driven by an increase of clearcut areas ( Table 3, Table S2). Forest gain was mainly concentrated near agriculture/developed www.nature.com/scientificreports/ areas, whereas forest loss tended to follow the distribution of shrub/bush and secondary forest in 2016 occurring mainly in the southern extension, the proposed protected area (Fig. 4). The net loss of forest interior area was at least 107 thousand ha, regardless of the scale of analysis, with a maximum loss of 155 thousand ha for the 590.5-ha neighbourhood size ( Table 4). The rate of loss of forest interior area increased monotonically with neighbourhood size (except for the largest one) and was approximately 2 to 3 times larger than the rate of loss of total forest area. Even though deforestation has occurred in the HRL, forest interior was prevalent with little fragmentation (Fig. 5).

Spatial predictors of forest change.
The multinomial logistic regression models showed that physical and environmental variables were the main drivers of either remaining in the current land-cover class or transitioning to another class (Fig. 6, Table S3). Forest was most likely to remain as forest (94% probability) and its persistence was related to slope and landform (i.e. steep land). Forest conversion, although less likely (< 6%), was    www.nature.com/scientificreports/ predicted by soil degradation and temperature. Areas with mild temperatures had a higher risk of forest conversion and becoming agriculture/developed land. Shrubland was most likely to remain shrubland (59%) and its persistence was associated with soil conditions, landform and temperature. Shrubland conversion to agriculture/ developed land, although less likely (8%), was predicted by soil degradation and distance from towns. Shrubland close to towns, in valleys or plain areas, and with degraded soils were more likely to become agriculture/developed land. Agriculture/developed land was most likely to remain stable (68%) and its persistence was associated with soil conditions. Conversion from agriculture/developed land to forest (12%) was predicted by temperature. Low temperature areas were more likely to become forest.

Discussion
Over 80% of the HRL is covered by natural areas, from which forest is the most prevalent land cover. Agriculture/ developed land was low (6%) and occurred mainly at lower elevations in the larger valleys around Putao and Naung Mung. Although the extent of forest degradation represented by shrubs and secondary forest increased from 1989 to 2016, the deforestation rate over the last three decades was low (annual average rate: 0.225%). The annual deforestation rate increased slightly compared with the deforestation rate reported for the region between 1991 and 1999 (annual average rate: 0.01%) 39 . However, the apparent increase in deforestation should be considered with caution. The period, the extent of the study area, the baseline dataset and the analytical methods diverged among studies making a plain comparison difficult. Nonetheless, the increased deforestation rate might indicate that forest loss surged since 2000. Forest loss was mainly located around the Naung Mung plains, along major roads, west of the Putao plains, and the western area of the proposed southern extension. This forest loss was due to the conversion from forest to shrub/bush and secondary forest indicating clearcuts, logging, and extraction of resources. Interestingly, we found a large part of forest loss areas replaced by shrub/bush and secondary forest on high altitude areas. This might be explained by a large wildfire that occurred in 1998 in the Hkakabo Razi National Park between Tazungdam and Tahaundam (own unpublished data Thein Aung). During the wildfire, large tracts of forest burned down. In 2001, fern covered the burned sites, in 2014 shrub/bush vegetation of roughly 5 m height emerged. Only minor  www.nature.com/scientificreports/ Compared to other forested regions in Myanmar, the annual deforestation rate of the HRL is still very low. For example, the deforestation rate at Chattin Wildlife Sanctuary from 1973 to 2005 was 1.86% 86 , an order of magnitude higher than for the HRL. Although, the annual deforestation rate for all forests in the country was between 0.55% between 2002 and 2014 16 and 0.62% from 2005 to 2016 87 , some areas have experienced annual deforestation rates above the country's average. Only 38% of the country's forests can be considered intact with canopy cover > 80% 16 .
Our results show that even though some deforestation occurred in the HRL, forest interior was prevalent with little fragmentation. This suggests that the remoteness of the area and the difficult access have played a role in safeguarding these large tracts of forest. So far, the only access is by air through Putao airport or the Myitkyina-Putao road. However, the recent paving of the Myitkyina-Putao road has reduced the time of transportation to less than two days instead of up to seven days of travel prior to 2014. The latent plans to connect the existing dirt road Putao-Pangnamdim with Yunnan/China would allow easy access to the site and pave the way for subsequent forest losses. Compared to the access from within Myanmar, accessibility from abroad is relatively easy. The lack of efficient in-country accessibility into the HRL and particularly the Hkakabo Razi National Park is regarded as unjust by the local communities (mainly Rawan, Lisu, and Tibetans) and causes major irritations (personal communication, September 2019 with the Rawan Cultural and Literary Organisation Putao). Local Rawan organisations in the HRL learned during the opening up of Myanmar since 2015 that accessibility would allow extraction of natural resources; major irritations come from who has the right to use and extract natural resources in HRL. The lack of efficient in-country accessibility has already resulted into tensions between conservation stakeholders, local organisations representing the Rawan, other local ethnic/religious groups entering into HRL, and government organizations. Addressing and balancing accessibility and extraction of natural resources are central for future management plans and conservation measures 42 to ensure the protection of these large tracts of forest interior and thus the long-term future of the country's forests.
Regression analyses showed that remaining in the current land-cover class was more likely than land conversion. Physical variables (i.e., slope and major landform) were associated with stable forest. As has been shown previously, forests that are difficult to reach have low probability of conversion, as they are located at the top edge of mountains or along steep cliffs 80,[88][89][90] . These results highlight that although major land conversion is occurring in other areas of the country 16 , the HRL has experienced low human impact and may be especially vulnerable to forest changes if sustainable forest policies and management are not put in place. It is important to note, though, that addressing the risk of deforestation does not only encompasses physical, environmental and geographical variables, but also interrelationships among political, institutional, economic, and ethnic and cultural factors should be considered 21,91 .
Based on our assessments, the forests of the HRL have experienced low human impact and still constitute large tract of contiguous forest interior. Therefore, the HRL should be regarded as of high integrity as for 2016 sensu Belle et al. 45 and should have high priority for protection. However, factors threatening the forest are emerging slowly, but steadily and with potential for an exponential increase. Thus, the HRL needs an inclusive forests management plan and effective conservation strategies to support biodiversity, guarantee the sustainable subsistence of local communities, and prevent irreversible forest loss and degradation.

Data availability
The datasets generated during the current study are available from the corresponding author upon request.