UC Merced Frontiers of Biogeography Do sub-groups of butterflies display different elevational distribution patterns in the Eastern Himalaya, India?

Understanding the pattern of biodiversity along environmental gradients helps in identifying diversity hotspot areas that can be prioritized for conservation. While the elevational distribution of several taxa has been studied, responses of the sub-groups within a taxon to elevation and its associated factors are not properly understood. Here we study species richness and butterfly density along an elevation gradient in Sikkim, Eastern Himalaya, India and explore the underlying causes of the patterns. We sampled butterflies using a fixed-width point count method in 16 elevational bands (150–200 m intervals), between a range of 300 and 3300 m a.s.l. We categorized butterflies into various sub-groups based on family, range size, biogeographic affinity, and host-plant specialization. We recorded 3603 individuals and 253 species of butterflies after the completion of 1860 point counts. Overall, species richness in the majority of the sub-groups (except for Riodinidae and Palearctic species) declines with elevation, as does the density of almost all the sub-groups. From a selection of environmental factors, annual actual evapotranspiration has the strongest effect on the species richness pattern of butterflies as well as on the density of the overall butterfly community, especially the Lycaenidae family. The richness and density of butterfly groups display varied responses to the richness and density of trees and shrubs. The conducive climatic conditions and diverse habitats in the lower valleys of the Eastern Himalaya support a high diversity of butterflies (with majority of small range species) and thus warrants conservation attention.


Introduction
There has been an upsurge in studies assessing biodiversity patterns across broad spatial scales, explaining the underlying processes, and exploring any conservation implications (Stevens 1992, Sánchez-Rodríguez and Baz 1995, Rahbek 2005, Acharya et al. 2011a, Wu et al. 2013a, Li and Feng 2015, Rana et al. 2019, Supriya et al. 2019. These studies may serve as a baseline for understanding the response of biological assemblages to climate change (Hodkinson 2005). Additionally, environmental gradient studies Keywords: Biodiversity pattern, biogeographic affinity, butterflies, elevational gradient, environmental factors, Himalaya, spatial factors.

This paper is part of an Elevational Gradients and Mountain Biodiversity Special Issue
Frontiers of Biogeography 2021, 13.3, e49643 © the authors, CC-BY 4.0 license 2 help to identify diversity hotspots that need to be prioritized for conservation (Hunter and Yonzon 1993, Bhardwaj et al. 2012, Chettri 2015. The majority of studies focusing on the latitudinal gradient have consistently found a decline in species from the equator to the poles (Pianka 1966, Stevens 1989, Gaston and Blackburn 2003. Elevational gradients are thought to proxy latitude in their distribution pattern of biological biomes and are therefore used to elucidate patterns and processes that influence biodiversity of various taxa around the world (Rahbek 1995, Gaston and Blackburn 2000, Colwell et al. 2004, McCain and Grytnes 2010, Sanders and Rahbek 2012. Patterns along elevational gradients are less general, however, due to complex biophysical processes that are involved in shaping the fine-scale spatial patterns on mountains. Commonly observed patterns of species richness along an elevational gradient are: (i) monotonic decline with increasing elevation, (ii) mid-elevation peak, (iii) low-elevation plateau and then linear decline, and (iv) low-elevation plateau with mid-elevation peak (McCain 2004, Rahbek 2005, Weins et al. 2007).
Factors such as contemporary climate, habitat heterogeneity, evolutionary events, and area or space have been proposed to explain diversity patterns along elevational gradients (Wiens et al. 2007, McCain andGrytnes 2010). The climatic factors include temperature, precipitation, humidity, and cloud cover (Rosenzweig 1992, Sánchez-Rodríguez and Baz 1995, Despland et al. 2012. These factors regulate the productivity and water-energy dynamics of mountainous ecosystems (O'Brien 2006, Wu et al. 2013b, Hu et al. 2017, Vetaas et al. 2019. The dynamic interaction of energy and water that controls liquid water availability may explain the richness patterns of various taxa (Wu et al.2013b, Hu et al. 2017, Vetaas et al. 2019. Similarly, energy drives the net primary productivity that facilitates species richness (Wright 1983, Rosenzweig 1992, Currie et al. 2004. The habitat heterogeneity hypothesis predicts that complex habitats provide more niches and diverse resources, thus increasing species richness (Bazzaz 1975, Kerr 2001, Levanoni et al. 2011. The species− area relationship predicts that along a montane gradient, environmental zones cover a wider area at the base and thus can sustain more species than zones such as mountain tops, which cover smaller areas (Rosenzweig 1992, Sánchez-Rodríguez andBaz 1995). The disparity in patterns between various organisms and associated factors makes it difficult to develop a universal model for explaining the variation of biodiversity (Supriya et al. 2019, Vetaas et al. 2019). Therefore, taxon or regionally specific studies need to be conducted so as to develop more specialized models.
The Himalaya forms the world's largest mountain complex and offers a unique system for understanding the elevational patterns of species richness. The Himalayan region is part of the globally significant biodiversity hotspot (Critical Ecosystem Partnership Fund 2020) and forms a transition zone between the Oriental and Palearctic biogeographical realms (Mani 1974, Holt et al. 2013). The region is well-suited for diversity studies along elevational gradients due to a sharp but continuous transition in vegetation and climatic conditions within a small geographical range. Much work has been carried out in the Himalaya on the distribution patterns of various taxa (Vetaas et al. 2019 and references therein), including plants (Grytnes and Vetaas 2002, Bhattarai et al. 2004, Oommen and Shanker 2005, Acharya et al. 2011b, Sharma et al. 2019, fishes (Li et al. 2009), amphibians (Fu et al. 2006, Chettri and, reptiles (Chettri et al. 2010), birds (Acharya et al. 2011a, Wu et al. 2013a, and mammals (Wu et al. 2013b, Hu et al. 2017. Pollinators, including butterflies, are facing a threat of extinction worldwide due to global climate change and anthropogenic activities, resulting in a pollination deficit which might trigger food shortages in the future (Allen-Wardell et al. 1998, Vanbergen 2013. Prior information on their distribution and influential variables forms the basis for predicting their responses to climate change and human disturbances. Such information would also be crucial in identifying hotspots that should be designated to conserve viable populations of taxa in the face of recent threats (Hunter and Yonzon 1993, Whittaker et al. 2005, Bhardwaj et al. 2012, Chettri and Acharya 2020. Biogeographical studies are particularly important in mountain regions such as in the Himalaya where the effects of climate change and human disturbance are more pronounced than in other parts of the world (Singh et al. 2011). In the Himalayan region, biogeographic studies on butterflies are scarce (Bhardwaj et al. 2012, Acharya and Vijayan 2015, Chettri 2015, Vetaas et al. 2019 and little is known either quantitatively or qualitatively about butterfly richness at local and regional scales (Rahbek 2005, McCain andGrytnes 2010).
Butterflies can be grouped into meaningful ecological sub-groups (Oommen and Shanker 2005) and can be studied at both the whole and sub-group level. Range size, biogeographic affinity, taxonomic categories, and feeding guilds have been consistently used in macroecological studies to categorize a wide variety of taxa such as plants, moths, mammals, birds, and amphibians (Oommen and Shanker 2005, Beck and Chey 2006, Fu et al. 2006, Wu et al. 2013b, Hu et al. 2017, Maicher et al. 2018, Zhou et al. 2019, Chettri and Acharya 2020, Subedi et al. 2020. Grouping taxa into smaller subsets may unveil differences in diversity patterns and their responses to abiotic and biotic factors. For example, the distribution of largerange species is likely to be affected by geographical constraints, whereas small-range species are more influenced by environmental factors (Fu et al. 2006, Wu et al. 2013b, Hu et al. 2017. Additionally, small-range or endemic species are often rare when compared to large-range species and might therefore be more vulnerable to extinction (Elsberry et al. 2018). Similarly, species having a tropical biogeographic affinity are more narrowly distributed while temperate species show a wider distributional range across elevation due to their higher environmental tolerance Feng 2015, Zhou et al. 2019 are also thought to influence elevational range size; however, no concrete evidence exists for such an assumption (Brehm et al. 2007).
Here, we aim to address the following questions. What is the pattern of richness and density of butterflies along the elevation gradient of Rangeet Valley, Eastern Himalaya? What are the factors associated with these patterns? How do the different sub-groups (categorized according to their family, elevational range size, biogeographic affinity, and larval host-plant specificity) respond to elevation? Are the underlying mechanisms similar in all the sub-groups? What are the conservation implications of the data?

Study area
The study was conducted in Rangeet Valley situated in the south-west district of Sikkim, Eastern Himalaya, India (Fig. 1). The elevation of the valley ranges from 300 m to 8586 m (the height of Mount Khangchendzonga) above sea level (a.s.l.). Rangeet River, an important tributary of the River Teesta, flows through the valley. The river originates from the Rathong glacier as Rathong Chu River at around 4674 m a.s.l. in west Sikkim and merges with the River Teesta near Melli at around 300 m a.s.l. The region experiences a broad range of climatic conditions from hot tropical at low elevations to cold alpine at high elevations. The temperature shows a linear decrease with increasing elevation (lapse rate = -0.50°C per 100 m) while precipitation is highest at mid-elevations (see Supplementary Figure S1). The rapid transition of climatic conditions along the elevational gradient influences the type of vegetation growth from tropical forest at lower elevations to alpine at higher elevations.

Butterfly sampling
We used fixed points along selected transects to sample butterflies following the methods of Acharya and Vijayan (2015). Mainstream methods such as Pollard's walk and transect count (Pollard 1977, Isaac et al. 2011 were not feasible because of the steep topography. Our method has been used in a considerable number of previous studies (Chettri et al. 2018, Dewan et al. 2019a, Sharma et al. 2020 and is recognized as an appropriate technique to sample We covered all the months to ensure that the majority of the species were encountered and recorded during sampling. We conducted a total of 1860 point counts during this study and the detail of sampling effort is provided in Table 1. Sampling was conducted on clear sunny days between 10:00 hrs and 13:00 hrs to ensure optimal weather conditions at least within one season. Additionally, we conducted the point counts in alternative order along the transect (i.e., starting from 1 st point in first sampling but with the last point in the next sampling, and so on) in order to avoid any time bias with respect to any particular point (although variation in weather conditions between this short time frame, i.e., 10:00 h to 13:00 h would be minimal). Butterflies were identified during sampling using the illustrated guide-books of Haribal (1992) and Kehimkar (2016). Butterflies that could not be instantly identified in the field were photographed and later identified by referring to guide-books and the ifoundbutterflies 1 website.

Species grouping
Butterflies observed in the study areas were broadly grouped according to their respective families, range size, biogeographic affinity, and hostplant specialization. Families represented include Nymphalidae, Hesperiidae, Lycaenidae, Papilionidae, Pieridae, and Riodinidae. In terms of the range-size category, butterflies with an elevational range greater than half (1500 m) of the total elevational range covered in the study (3000 m) were considered as large-range species and the rest as short-range species, following the approach adopted by Wu et al. (2013b). Species observed in only one elevation band were assigned a 100 m range (±50 m of the point elevation), assuming the species to be present within this range (Stevens 1992). Holloway (1964Holloway ( , 1974 ascertained the distribution ranges of selected Indian butterflies and identified the centres of their diversity. Following a simplified version of Holloway's method (Kunte 2011), we grouped the observed butterflies according to their affinities to respective biogeographic realms into: (a) global (having a centre of diversity in at least two regions), (b) Oriental (affinity to hot, humid, evergreen forest habitats), and (c) Palearctic (affinity to colder and temperate regions). A few species with an affinity to the African region, and others that did not show affinity to any biogeographic realm, were excluded from the group-based analyses. Based on larval host-plant specialization, we grouped butterflies into monophagous (larva feeding on plants in only one genus), oligophagous (larva feeding on plants in a single family, but more than one genus), and polyphagous species (larva feeding on plants in more than one family or order) (Zhang 2019). The data on host plants were obtained from various sources (Haribal 1992, Kehimkar 2008 supplemented by field observations. Data on host plants of 72 species were deficient and larvae of two species were carnivores, hence, analysing the trends of these species was not viable.

Area
We downloaded Digital Elevation Model (DEM) imagery (covering the Sikkim Himalayan region) generated from the Cartosat-1 satellite, built and operated by the Indian Space Research Organization (ISRO). The Cartosat-1 data are freely available on Bhuvan 2 , an online Indian geospatial platform. The DEM raster image was first classified into the new classes that correspond to the respective elevation bands (classified according to this study) using the reclassify tool in the spatial analyst toolbox in ArcGIS 10.4. We then calculated the area of each reclassified elevation band using the zonal geometry tool in ArcGIS 10.4. Zonal geometry calculates the geometry measures (e.g., area) for each band in a dataset.

Vegetation
Trees and shrubs were surveyed once during the study period along all the 16 transects established for sampling butterflies. Quadrats of 10x10 m were established at each point, adjacent to the butterfly sampling point, and two sub-quadrats (5x5 m) were laid diagonally within each 10x10 m quadrat for surveying shrubs. Plants with DBH (diameter at breast height) ≥20 cm were considered as trees. We estimated species richness and density of trees and shrubs in the quadrats and then pooled the data for each elevation band.
The Normalized Difference Vegetation Index (NDVI) was used as a surrogate for above-ground productivity (Nieto et al. 2015). We used three years (2016-2018) of Landsat8 imagery data (available at 30 m resolution) for the Sikkim Himalayan region available from the USGS 3 website. We first averaged three years of red and near-infrared imagery data and then calculated NDVI from these averaged outputs using the formula

Climatic variables
All climatic data were downloaded from the CHELSA (Climatologies at high resolution for the Earth's land surface areas) dataset (Karger et al. 2017a, b). Among 19 bio-climatic variables, mean annual precipitation (MAP) and mean annual temperature (MAT) were selected because they are important in affecting the distribution of butterfly biodiversity. Besides, other bio-climatic variables are mostly derived from these two variables. The CHELSA dataset has a resolution of 30 arc seconds (1 km 2 grid). The selected variables were then categorized into equal elevation bands of 200 m, except for the lowest band (<500 m) as land below 300 m does not occur in the study area. Temperature and precipitation values for consecutive elevation bands within the whole Sikkim Himalaya area were obtained by averaging the grid values falling into each band using ArcGIS 10.4. We also calculated annual evapotranspiration for all the transects using standard equations such as provided by Turc (1954) for actual evapotranspiration (AET;cf. Kluge et al. 2006) and Holdridge et al. (1971) for potential evapotranspiration (PET). AET is a function of water availability and temperature and, hence, has been used as a measure of water-energy balance (Bini et al. 2004).

Data Analysis
Observed species richness is the total number of species observed in all seasons per elevational band during the study. Since it is practically impossible to detect all species present, simple observed species richness may not always be a reliable estimate of richness. Hence, non-parametric estimators were also used to estimate richness (Colwell 2013). Chao1 and Jackknife1 estimators were selected owing to their high precision in estimating richness (Hortal et al. 2006). Species accumulation curves were generated using these estimators in order to assess the completeness of the sampling effort. To reduce the bias of unequal sampling effort, we estimated sample-based rarefied richness. Here, species richness was rarefied to the lowest number of counts conducted for any site (110 point counts). A preliminary analysis showed that Jackknife1 predicted a higher number of species in most of the sites. Hence, we used Jackknife1 estimated richness as the measure of species richness.
We also recorded abundance (total number of individuals) of butterflies in each elevation band. To account for variability in abundance due to unequal sampling, we converted the abundance records into density of butterflies in each elevational band. Density is the number of individual butterflies recorded per where D = butterfly density (numbers ha -1 ), n = total number of butterflies observed in all counts within the specific radius, r (specific radius is the average radial distance of butterflies from the observer), and C = total number of counts conducted, following the approach used for birds (see Reynolds et al. 1980). From the overall pooled data, species richness and density of all the sub-groups in each elevational band were estimated.
To assess the relationship between elevation and observed species richness, estimated richness, rarefied species richness, total density, and species richness and density of the sub-groups, we drew scatter plots using the ggplot2 package (Wickham 2016) in R (R Core Team 2019). Richness estimates and density of most sub-groups show a linear trend with elevation: hence, we used ordinary least squares regression to test the significance of the relationship. We also analysed the relationships between species richness and density with various predictor variables. Since MAT (r = 0.998, p <0.01), MAP (r = -0.874, p <0.01), and PET (r = 1, p <0.01) are highly correlated to AET, we do not consider MAT, MAP, and PET for further analyses (see Currie et al. 2020). Generalized linear modelling (GLM) with a log link function assuming a Poisson distribution of error was used to explore the relationship of selected explanatory variables and total species richness or density of butterflies in the various sub-groups (cf. above). A total of 128 GLMs were generated using the package glmulti in R (Calcagno and de Mazancourt 2010). From all the models generated, the best fitting GLM is the one with the lowest corrected Akaike information criterion (AIC c ) value. Models with a ΔAIC c <2 from the model with the lowest AIC are considered equally likely (Burnham and Anderson 2002). Hence, we used a model averaging approach to compare all the likely models and estimated the relative importance of each of the predictor variables from these models (Johnson and Omland 2004) using the package MuMin in R ( Barton and Barton 2013). GLM takes into account the deviance explained in each of the models because the lowest AIC c has the minimum residual deviance.

Species richness and density of butterflies along the elevational gradient
A total of 3603 individual butterflies representing 253 species and six families were recorded during the study ( The observed species richness of butterflies shows a declining trend with increase in elevation ( Table 2, Fig. 2). For each elevation band, estimated richness (Jackknife1) predicts a slightly higher number of species, suggesting more species could be counted with further sampling. The species accumulation curve, however, predicts that the rate of addition of species would be uniformly low, thus suggesting that the sampling effort was almost complete (Supplementary Figure S2). Estimated species richness (Jackknife1) also shows a declining trend but with a slight hump at around 500 m elevation (R 2 = 0.868, p < 0.01). Similarly, rarefied species richness shows a declining trend with increasing elevation (R 2 = 0.883, p < 0.01) ( Table 2).
The species richness of the majority of the subgroups declines with increasing elevation and fits well to a linear declining model ( Table 2, Fig. 2). When assessed, as butterflies grouped into higher taxonomic levels, the results are similar to total species richness for families such as Nymphalidae (R 2 = 0.806, p < 0.01), Papilionidae (R 2 = 0.833, p < 0.01), Hesperiidae (R 2 = 0.684, p<0.01), Lycaenidae (R 2 = 0.890, p < 0.01), Pieridae (R 2 = 0.768, p < 0.01) with the exception of Riodinidae, which shows no definite trend. The species richness of the small-range butterflies shows a distinct linear decline with elevation (R 2 = 0.836, p < 0.01), while the large-range species have two distinct peaks (at 500 m and 1700 m), making a poor fit to the linear regression model (R 2 = 0.649, p = 0.01). Global (R 2 = 0.844, p < 0.01) and Oriental (R 2 = 0.836, p < 0.01) species mirror the overall species richness pattern and decline linearly with an increase in elevation but Palearctic species do not show any definite pattern. Oligophagous (R 2 = 0.909, p < 0.01), monophagous (R 2 = 0.583, p = 0.01), and polyphagous (R 2 = 0.786, p < 0.01) butterflies show a declining trend with elevation confirming that butterflies of the study area decline linearly with the elevation.

Determinants of butterfly species richness and density
The overall richness of butterflies is explained by two sets of best candidate models that have the lowest AIC c (Supplementary Table S2    -0.050 0.020 0.320 -2.560 0.022 * Coefficient of regression, standard error (Std. Error), R 2 representing the proportion of variance of regression, and t-value along with overall significance of the regression are presented. *Significant at p <0.05, ** significant at p <0.01. Negative relationships are indicated by a minus (-) sign. inference of the two most likely models suggest that AET, followed by tree species richness and density are the best explanatory variables of the major variation in overall butterfly richness along the elevation gradient (Table 4). Other sets of models explain the richness patterns of butterflies of the different sub-groups. AET significantly influences the species richness pattern of most sub-groups of butterflies except for Riodinidae, large-range, and Palearctic species. Tree species richness strongly affects the species richness patterns of Nymphalidae, Hesperiidae, small-range, and monophagous butterflies, whereas tree density is an important determinant of oligophagous butterflies. The species richness of the Riodinidae family, largerange, and Palearctic butterflies shows no significant relationship with spatial or any environmental variables.
For the density of butterflies, average model sets suggest that AET followed by shrub density are the most significant predictors of total butterfly density along the elevation gradient ( Table 5, Supplementary  Table S3). Amongst the different sub-groups of butterflies, AET is a significant variable for the density of Lycaenidae only. Habitat variables such as shrub density  significantly affect the density of Pieridae, smallrange, oligophagous, and polyphagous butterflies. Similarly, tree density significantly affects the density of Hesperiidae butterflies. However, species richness of shrubs is found to have a negative influence on the density of certain sub-groups of butterflies such as Papilionide, Pieridae, and Palearctic species.

Species richness and density along the elevation gradient
This study examines the pattern of species richness and density of butterfly communities along an elevation gradient in Rangeet Valley in Sikkim, Eastern Himalaya, India. The species richness and density of butterflies generally decline with an increase in elevation. A midelevation peak is the most common pattern of species richness in mountain ecosystems for the majority of taxa (Rahbek 1995, McCain and Grytnes 2010, although for butterflies, a monotonic decline has been frequently reported from the Himalaya and elsewhere (Sánchez-Rodríguez and Baz 1995, Kumar et al. 2009, Bhardwaj et al. 2012, Leingärtner et al. 2014, Acharya and Vijayan 2015, Chettri 2015. A monotonic decline in species richness with increasing elevation might therefore be the general pattern for butterflies.
While the species richness and density of the different sub-groups often mirror the overall richness and density patterns, we find a few exceptions to this general pattern. We also find varied responses of the sub-groups of butterflies to spatial, environmental, and biotic variables. Differences between the subgroups indicate that the trends strongly depend on the sub-groups or species considered (Wu et al. 2013b). The variety of trends and responses to explanatory variables may be attributed to differences in physiological adaptation, ecological requirements, and the evolutionary history of the species groups (Wu et al. 2013b, Zang 2019. The richness and density of five butterfly families, namely Nymphalidae, Papilionidae, Pieridae, Hesperiidae, and Lycaenidae, follow a declining trend. Hesperiidae and Papilionidae are mostly restricted to an elevation below 2000 m, probably due to physiological requirements for their energetic lifestyle. The species richness and density of small-range butterflies decrease linearly with an increase in elevation, whereas large-range species do not show a clear linear decline. Several studies have shown that small-range species are likely to be affected by environmental variables while large-range species (having wider environmental tolerances) with humpshaped diversity distribution might also be influenced by geographic constraints (Jetz and Rahbek 2002, Colwell et al. 2004, Brehm et al. 2007. Larger ranges are more likely to overlap in the middle of the domain causing a mid-elevation peak in richness (Colwell and Hurtt 1994). This perhaps explains the distinct midelevation peak in richness as well as density of the large-range species found in our study.
The Lepidoptera of the Eastern Himalaya are mostly dominated by Oriental species (mostly Indo-Chinese and Malayan forms), with less representation of Global and Palearctic elements (Holloway 1974, Mani 1974. The Oriental biotas are mostly represented by species adapted to the tropical hot/humid climate, whereas Palearctic elements are considered to be representative of the colder temperate region (Holloway 1974). The differences in niches of the Palearctic and Oriental biota can be observed in the Himalayan butterflies. The mixing of faunal elements having different biogeographic affinities provides direct evidence that historical events such as continental drift, Himalayan uplift, and colonization were important in shaping the current distribution of butterflies (Miehe et al. 2015).
While we find distinct variation in richness and density patterns between many sub-groups, there are no differences in the trends between butterflies categorized according to their feeding specificity as species richness and density always decrease with increasing elevation. The elevational niche-breadth hypothesis predicts that the diet breadth of herbivores increases with increasing elevation (Rasmann et al. 2014); hence, it will be pertinent to assume that a higher number of species will be polyphagous at higher elevations while lower elevations will have more specialist species. The deviation in our results from this hypothesis may be due to: (i) non-availability of sufficient information on larval host plants for Himalayan butterflies, and (ii) the exclusion of a large spatial extent of alpine area (>4000 m) in our study due to logistical reasons (harsh climatic conditions, steep gradient, accessibility, etc.). Species in stressful habitats (such as alpine area in our study) are more likely to have different life-history strategies compared to their lowland counterparts. Evidence for the niche-breadth hypothesis is mixed and varies according to region. Pellissier et al. (2012) show that in temperate climates, diet breadth of butterflies decreases with elevation, while Rodríguez-Castañeda et al. (2010) find the opposite pattern in the tropics. Novotny et al. (2005) find no significant difference in moth diet-breadth with elevation in the tropics. More research is necessary to understand how species with different dietary requirements are segregated along environmental gradients (for example, elevation).

Determinants of species richness and density along the elevational gradient
Among all the variables, annual AET is the most important variable affecting the overall species richness patterns and total density of butterflies along an elevation gradient in the Eastern Himalayan landscape. Annual AET has been found to strongly influence butterflies (Acharya and Vijayan 2015) and trees (Acharya et al. 2011b, Rana et al. 2019) along elevational gradients in the Eastern Himalaya. AET is reported to decline with elevation (Trabucco and Zomer 2010), resulting in the decline in species richness of butterflies. AET is known to function in two ways -(1) directly affecting the physiology of organisms via temperature/light stress and water availability (water-energy balance or water-energy dynamics) and (2) by influencing the productivity of the ecosystem (Rosenzweig 1995, Waide et al. 1999, Hawkins and Porter 2003, Whittaker and Heegaard 2003. Waterenergy dynamics has subsequently been demonstrated to be a better explanation than net primary productivity in explaining species richness patterns of various taxa, including butterflies (Vetaas et al. 2019). Since butterflies are ectotherms, thermal energy is crucial to their basic physiology and their feeding behaviour by influencing water availability in all forms (nectars, mud puddles, fruit juices; Fleishman et al. 2005, Kehimkar 2008). It can, therefore, be concluded that AET affects species richness both indirectly by influencing primary productivity and, most importantly, directly due to the physiological requirements of the butterflies. Largerange and Palearctic species, being widely distributed, are less affected by the AET gradient due to their higher level of environmental tolerance. Also, the density of most of the sub-groups seems less affected by AET and more by habitat variables, indicating that resource abundance is necessary to maintain the population of the species (Curtis et al. 2015).
Resource availability and habitat condition are also considered a strong determinant of species richness and density (Ribas et al. 2003, McCain andGrytnes 2010). We find a strong relationship between habitat variables (tree species richness, tree density, shrub density) and the species richness pattern and density of butterflies. The ambient climatic condition  (Schulze et al. 2004, Vu 2009). Monophagous butterflies, in particular, show a strong relationship with tree richness, indicating that their distribution is mostly affected by host-plant distribution. Moreover, it is evident that higher plant diversity at lower elevations results in a more heterogeneous habitat, resulting in an increase in butterfly diversity. Habitat heterogeneity also influences species richness because complex habitats provide more diverse resources, thus increasing species diversity (Bazzaz 1975). An increase in area has often been linked to an increase in species richness (Rosenzweig 1992). At the regional or global scale, the extinction rate decreases due to more populations in larger areas, and speciation increases due to the potential for the formation of barriers. At the local scale, larger areas support more diverse habitats, allowing more species and individuals to thrive. Thus, along an elevational gradient, the species-area relationship may matter along with these two scales (Rosenzweig 1995, McCain 2007. In this study, we do not find a statistically significant relationship between area as a potential driver and species richness or density of the different butterfly sub-groups. In contrast to most mountain areas, where area decreases with increasing elevation, in the Sikkim Himalayan region areal extent has two distinct peaks (500 and 1500 m, Supplementary Figure S1). Studies in other parts of the Himalayan region show similar results where the relationship between area and the species distribution pattern along an elevation gradient are not significant (Hu et al. 2017).

Conclusions and Conservation Implications
Butterfly species richness and density decline with an increase in elevation, with the highest values below 500 m in Rangeet Valley in Sikkim, Eastern Himalaya. The trends in species richness and density and associated biotic, abiotic, and spatial factors vary with the sub-groups considered. This indicates that sub-groups within a taxon may respond differently to climatic changes and anthropogenic pressures. We find that the trends in species richness and density are mainly explained by climatic factors and habitat variables. Reports on the range shift of butterflies due to global climatic changes are on the rise (González-Megías 2008, Foristera et al. 2010, Braby andHsu 2019). Small-range species, Oriental species, the majority of the butterfly families, and polyphagous species are more likely to be affected by changes in temperature and precipitation gradients caused by climate change. Additionally, monophagous species, due to their exclusive dependency on habitat variability, are also threatened by habitat loss which will exacerbate the effect of climate change (Fonseca 2009). Such differences in resilience and vulnerability are mainly due to the variation in life-history associated with a particular group. Long-term studies are needed to document the life-history traits of the various butterfly sub-groups and to understand their responses to the energy-elevation gradient. Such studies will provide significant insights to inform better directed conservation policies for different groups of butterflies.
Climatic conditions and diverse habitats in the lower sub-tropical valleys of Rangeet support a high diversity of butterflies along with majority of smallrange species, and thus requires high conservation attention. The lowland forest in the Himalaya and elsewhere is under immense anthropogenic pressure, leading to the extinction of species (Pandit et al. 2007). In the Sikkim Himalaya, 31% of the total geographical area is within a protected area network, but most of the protected areas are above 1500 m a.s.l. (Forest, Environment and Wildlife Management Department 2019). The forest in lowland areas below 500 m in Sikkim covers only 40 km 2 (Forest Survey of India 2017). The low-elevation landscape is mostly dominated by agricultural lands, industry, dams, towns, and road networks. Apart from natural forests, the traditional agroecosystem has been shown to be important for the conservation of butterflies and odonates (Dewan et al. 2019b, Sharma et al. 2020. Hence, specific policies are required to safeguard these lowland landscapes (including traditional agroecosystems), which are vital for the conservation of butterflies in the long run.
Frontiers of Biogeography 2021, 13.3, e49643 © the authors, CC-BY 4.0 license 14 Figure S2. Species accumulation curves of butterflies observed at different elevations Table S1. Details of butterfly species recorded during the study Table S2. Best candidate generalized linear models describing the relationship between species richness of butterflies and selected predictor variables Table S3. Best candidate generalized linear models describing the relationship between density of butterflies and selected predictor variables