Study areas
Our study was conducted in three areas of conservation interest in central, western, and southern Myanmar (Fig. 1). Site 1 (Latitude: 17.1013 – 18.1960, Longitude: 95.7043 – 96.4787) is located in the central part of Myanmar in the foothills of Bago Yoma mountain ranges. Historical unsustainable teak extraction in this site created a highly disturbed forest mosaic that is increasingly being invaded by other human land uses, including the construction of hydroelectric reservoirs, settlements, as well as commercial teak, sugarcane, rice, and rubber plantations. Site 2 (Latitude: 16.0554 – 17.0842, Longitude: 94.1860 – 94.6838) is a mountainous area along the west coast of the Ayeyarwaddy state that stretches from north to south creating an elongated forest with hard boundaries on east and west. Rice plantations dominate the matrix between forest patches in this site, where rubber and peppercorn agricultural use is also prevalent. Site 3 (Latitude: 10.7141 – 12.0981, Longitude: 98.3356 – 99.4626) is part of the larger Dawna Tanintharyi Landscape which extends from mountain ridges along the border with Thailand to the coastal plain. Land use at site 3 is primarily composed of oil palm and betel nut plantations, surrounded by lowland deciduous forests. Threats of human encroachment, road development, and agricultural expansion into the remaining forest are rising in the area.
The study areas are strongly seasonal, with rainfall records demonstrating the extended dry season occurs between early December and late March, and the wet season between early June and mid-September [30]. During the dry season, human-elephant conflicts (HECs) peak in relation to the harvest of rice, sugarcane, and other agricultural products [29, 31]. Rainfall is significantly higher at site 3, resulting in markedly different forest composition. Forests at site 3 are predominantly lowland evergreen forests, while at site 1 and 2, they are mostly mixed deciduous forests with strong seasonal leaf-fall patterns.
Elephant capture for GPS collaring
All capture and animal sedation were performed by veterinarians from Myanmar Timber Enterprise (MTE). MTE is the Myanmar government agency responsible for the management of logging elephants, and their staff have extensive experience in veterinary care of captive and wild Asian elephants, including sedation. Individuals were independently captured during the collaring period, and no collared elephants were found in the same social unit. All capture and handling procedures followed or exceeded the guidelines provided by the American Society of Mammalogists [32]. Elephants were immobilized using Etrophine and Xylaxine for sedation and Naltrexone for reversal. The immobilization and collaring process took approximately 30 minutes per individual on average and was carried out early in the morning or late in the afternoon when air temperature was relatively lower (<35° C). All the collars are set to record a GPS fix every hour. Due to high collar failure and poaching soon after collar deployments [19], telemetry datasets were often patchy and covered relatively short periods. Consequently, we only included individuals with 1) > 60 days of tracking data and/or 2) that had an established range (based on a semi-variogram analysis of range stability described below). To assess whether animals established ranges during the tracking periods, we used methods described by [33] in their continuous-time movement modeling package (ctmm) in R. When the semi-variogram function for the relocation data of an elephant approached an asymptote, we classified that dataset as capturing an established range [33], which occurred within 60 days each season for the elephants in this study.
For dry season range analysis, we analyzed data from eight individuals from site 1 (4 females: 4 males), six from the site 2 (1 female: 5 males), and eight from site 3 (2 females: 6 males) – totaling 22 different individuals. We performed data analysis on data collected from December 2016 through March 2020. Therefore, our analysis covers four dry seasons. There were four individuals whose tracking periods covered two dry seasons. To avoid problems with pseudo-replication when developing our regression model set, we excluded the year with fewer data points for each of these four individuals such that each only supplied one season to the analysis.
For the wet season ranges estimation, we utilized data from five individuals from site 1 (1 female: 4 males), three from the site 2 (3 males), and two from site 3 (2 males). There were two individuals who had data spanning two wet seasons, allowing estimation of 12 wet-season ranges. Because of the relatively small sample size, we did not run a regression analysis on wet season data.
For annual home range estimation, we included individuals with > 365 days tracked, which amounted to 8 individuals: five individuals from site 1 (1 female: 4 males), one from the site 2 (1 male), and two from site 3 (2 males).
Range estimation
We employed a ctmm framework [33] to estimate seasonal (dry and wet) and annual range sizes among individuals. We compared the fit to our data of independent and identically distributed (IID), Ornstein-Uhlenbeck (OU), and Ornstein-Uhlenbeck Foraging (OUF) movement models using an autocorrelation estimation method. We picked the best fitting model and applied it to fit the autocorrelated density estimator (AKDE) function to estimate range size. We calculated 95% and 50% AKDE percentile level ranges for all individuals. We assumed 50% AKDE level as core areas within the respective ranges where animals spent 50 percent of their time. To enable comparison with other studies, we calculated and presented 95 percentile Minimum Convex Polygon (MCP) ranges.
Predictor variables and candidate models
To assess which landscape conditions were related to dry season range size, we applied gamma regression models with estimated AKDE range sizes as a response variable based on the assumption that our response variable can only be a non-zero positive number. Our covariate dataset of landscape properties was developed by classifying Landsat 8 imageries to develop land cover maps for each of our study areas (Chan et al. unpublished data). We used the ‘landscapemetrics’ package in R to derive different measures for characterizing landscape metrics from our land cover map [11]. To describe the landscape of each individual range, we calculated several different shape, area, edge, and aggregation metrics for water, agriculture, and natural vegetation classes (Table 1). In addition, we quantified landscape-level metrics, including Shannon’s diversity index, relative patch richness, and relative patch density (Table 1). We computed 48 landscape metrics in total.
To simplify these 48 metrics for our regression analysis, we relied on exploratory factor analysis with oblique minimal rotation of principal factor axes to reduce the data dimensionality. This approach relaxes the assumption of normality [34] and allowed us to identify the variables that best characterized variations in our landscapes (Table 1). Specifically, we included the highest positive and negative loading variables from the first five principal factor axes to reduce the metrics to the primary explanatory variables while explaining sufficient variance in the data. Afterwards, we compared single variable models among metrics belonging to the same land cover class. We kept the variables if the AIC corrected for small sample size (AICc) score was within 8 of the top model and excluded variables that did not meet the criteria in our candidate model set. AICc is the metric used to rank the models in your candidate model set in such a way that the most parsimonious model will have the lowest AICc value among the model set. This allowed us to eliminate variables with relatively low explanatory power. We also assessed the Spearman’s rank correlation coefficients between all the variables before including them in the final candidate model sets (all the variables included in the model set were less than 0.6).
From the retained variables (Table 1), we then developed different biologically meaningful combinations of agriculture and natural vegetation indices in the model set for both the 95% and 50% AKDE level for dry season ranges. We included a model with a quadratic term for percentage of agriculture to determine whether we could assess the threshold relationship between agriculture and range size. We also assessed the effect of sex, site, and year by adding these covariates to our best performing model and ranked the models using AICc for both 95% and 50% AKDE top models. We investigated these effects further in 50% AKDE range sizes analysis by dropping uninformative parameters in our model sets based on model weights and parameter estimates and presented the most parsimonious and biologically meaningful model since the effect of site/region came out stronger in our model set [35].
In addition to our range size models, we developed a candidate model set to assess the correlation between landscape metrics and average daily distance moved by the elephants. We calculated average daily distance moved by calculating the sum of hourly distance moved (straight-line distance between the two consecutive points) and dividing by the total number of days tracked for the particular individuals. We did not include days where fix success rate was below 80% in our daily distance traveled calculation. For this particular candidate model set, we tested several hypotheses using the most informative variables from the 50% AKDE dry season analysis. We tested whether sex, site/region, and/or two agriculture metrics (percentage of agriculture presence and perimeter-area ratio of agriculture patches within the range) influenced average daily distance moved by elephants by fitting gamma regression model as described above. We set female and study site 2 as a reference category for sex and region categorical predictor variables, respectively, in the model.
We used AICc to rank models in the candidate model set [36]. We selected the model with the lowest AICc as the best/top model in respective candidate model set. To account for variation in range sizes driven by sampling differences, we included the number of days tracked as an additional variable in the top model. We retained the number of days tracked variable if it was included in a model within 2 AICc scores of the top model. All variables were standardized to a mean of 0 and a standard deviation of 1 before fitting the model for easier interpretation of the results and standardize the effect size of all covariates. All analyses were conducted in R version 3.6.3 using ‘ggplot2’ (version 3.3.0), ‘dplyr’ (version 0.8.5), ‘ctmm’ (version 0.5.9), ‘landscapemetrics’ (version 1.4.3), and ‘AICcmodavg’ (version 2.2.2) [37–41].