Statistical analysis
To determine if the range shift (centroid distance: distance between the current and future centroids) is associated with the amount of current suitable habitat and elevational distribution, I ran two separate GLMMs in R using the ‘glmmTMB’ package (Magnusson et al., 2017). To calculate whether the range shift (centroid distance) is associated with the current suitability distribution along the elevational gradient, I added the centroid distance as the response variable and the interaction between mean elevation (current) and migratory status as the explanatory variable. I used the interaction of migratory status with the response variable to account for how migratory species can adjust to the new habitat faster than the non-migratory species and assess if Bangladeshi migratory butterflies follow a similar trend. To calculate whether the range shift is associated with the amount of suitable habitat (current), I ran a similar model considering the range shift as the response variable and the interactive effect of the amount of suitable habitat (current) and species migratory status as the explanatory variable.
In both cases, I used the ‘nbinom2’ family to account for overdispersion and butterfly families as a random factor. Finally, if any interaction was significant, I conducted post hoc pairwise comparison tests using the ‘emmeans’ R package (Lenth et al., 2019).
To assess the vulnerability status of the threatened and non-threatened species, I ran two separate GLMMs (for centroid distance and elevation shift) in R using the ‘glmmTMB’ package (Magnusson et al., 2017). To test whether threatened species experienced more range shift or elevation shift, I added centroid distance and elevation shift periodically as the response variables and the threat status as the explanatory variable. I used the ‘nbinom2’ family to account for overdispersion and butterfly families as a random factor. Given that there were negative values in the elevation shift (species for which the mean elevation declined), before fitting the model, I transformed the values by adding (1+ maximum negative values) for each species.