Complex elevational shifts in a tropical lowland moth community following a decade of climate change

Climate change is driving many species towards higher latitudes and higher elevations. However, empirical studies documenting these changes have largely focused on presence/absence based range shifts in temperate regions. Studies in lowland tropical ecosystems that control for detection probabilities are especially lacking.


| INTRODUC TI ON
Climate change is affecting biodiversity across the globe (Araújo & Rahbek, 2006;Garcia, Cabeza, Rahbek, & Araújo, 2014). Distribution shifts towards higher altitudes and latitudes have been observed in many taxa (Hickling, Roy, Hill, Fox, & Thomas, 2006;Pecl et al., 2017). Despite the large number of resurvey and monitoring studies that have compared distributions between historic and current communities (Chen, Hill, Ohlemüller, Roy, & Thomas, 2011;Lenoir & Svenning, 2015), much variation exists in range shift studies exclusively examining shift patterns; variation in methods particularly complicates the attribution and magnitude of range shifts detected across species groups, local climates and habitats (Guo, Lenoir, & Bonebrake, 2018).
Directly calculated distribution limits or range centre shifts often suffer from detectability issues, as imperfect detection and false absences can cause apparent range shifts even when no true shift occurs (Tingley & Beissinger, 2009). Variation in sampling effort across time (documentation of historic sampling effort is usually limited or unavailable) is especially important to consider in comparing distributions in the context of environmental change (Bates et al., 2014;Brown et al., 2016;Tingley & Beissinger, 2009). Furthermore, short-term studies are more likely to track population fluctuation signals, rather than real climate change effects, especially species with high population variability such as invertebrates and annual plants (McCain, Szewczyk, & Bracy Knight, 2016). Under such circumstances, the use of occupancy modelling to correct for imperfect detection when examining range shifts has gained popularity (Moritz et al., 2008;Rowe et al., 2015;Tingley, Koo, Moritz, Rush, & Beissinger, 2012). Occupancy modelling utilizes hierarchical logistic regression models to separate observation processes from occupancy processes to control detection probability variation commonly encountered in ecological sampling (MacKenzie et al., 2002;Tingley & Beissinger, 2009). Partitioning variation caused by observation allows for a more accurate analysis of ecological variation and has been shown to significantly reduce false absence rates compared to traditional direct comparisons (Tingley et al., 2012).
Here, using an annual continuous trapping dataset of geometrid moths, we investigated species range shifts along a lowland elevational gradient in Hong Kong, at the northern edge of the tropics. Moths of the family Geometridae are relatively well studied both biologically and taxonomically (Holloway, 1993;1996;1997;Kendrick, 2002). They are also ecologically important in natural systems (Schowalter, Hargrove, & Crossley, 1986) as well as sensitive to climate change (Chen et al., 2009) and habitat disturbance (Beck, Schulze, Linsenmair, & Fiedler, 2002) for detection probabilities with occupancy models, we tested (a) whether occurrence probabilities across elevation shifted for each species during our study period, (b) whether each species showed shifts of upper or lower elevation boundaries while accounting for false absences, and (c) whether the community of geometrid moths at re-sampled sites showed systematic elevational trends over time.

| Study system
Continuous sampling of geometrid moths in Hong Kong started in 1994 and continued through 2014 by one of the authors using Robinson light traps (Fry & Waring, 1996) (Dudgeon & Corlett, 2004). In total, 556 moth samplings events were opportunistically made between 1994 and 2014, with a maximum of 14 sampling events per year. During each sampling event, light traps were set from dusk to midnight (but sometimes for fewer hours) and species checklists were recorded.

| Local climate and habitat change
To quantify local climate change during the study period, we downloaded monthly averaged daily mean temperature, monthly averaged daily maximum temperature, monthly averaged daily minimum temperature and accumulated monthly rainfall from number of dry days (daily rainfall = 0 mm) and annual number of flood days (daily rainfall>100 mm) for Hong Kong from the same source. We were only able to find climate data from 1997 to 2014.
To quantify the changes of these climatic variables over the study period, we applied simple linear regressions with year as an explanatory variable and the yearly climatic variables as response.
Frost days were very rare during the sampling period, so we did not analyse these data statistically, but given the important ecological influence of frost on Lepidoptera (Wang et al., 2016), we qualitatively compared frost days over time.
We used a high-resolution (30 m by 30 m) forest cover change dataset as an indicator of land use change over the last decade (Hansen et al., 2013). We generated forest cover loss and gain pro-

| occupancy models and occurrence probability change
Using occupancy modelling in a Bayesian framework, we estimated and tested for changes in occurrence probability from 2000 to 2014 across elevation. Because sampling was highly uneven at sites across years, we undertook a series of data partitioning steps prior to analysis. We thus subset the moth dataset to two time  Additionally, because only a small number of sites (six out of 24) were visited in both time periods, we applied an "unpaired site" model parameterization as described by Tingley and Beissinger (2009), which looks for changes in occupancy across a defining gradient (here, elevation), rather than explicit colonization-extinction at paired sites. We only included sites visited more than once, resulting in 16 historical sites (92 sampling events) and 14 current sites (98 sampling events). Despite the fact that the geographic distances between sampling sites were greater in the historic sampling, the elevational distances are similar (Two-tailed T test, historic 177.4 m, current 167.8 m, p = 0.6). Detailed site information is shown in Figure 2 and Supporting information Table S1.
We only analysed species that occupied at least 10% of both historic and current sites, as occupancy models cannot reliably be fit in data-sparse situations. Consequently, we focused the analysis only on species for which we were confident that we could successfully fit an occupancy model, resulting in 123 species in total (Supporting information Table S3).
As all moths were collected by the same person with identical trap design, for each species, we assumed a constant detectability between the two sampling periods. While peak flying time for adult moths in Hong Kong is either at the beginning of the wet season (March to May), or at the end of the wet season (September to November), or both (Kendrick, 2002), sampling was concentrated during these seasons for both historic and current periods. Additionally, there was no seasonal bias in sampling within elevational bins across the two time periods (Supporting information Figure S1).
For occupancy parameterization, we included both sampling period and elevation in our models, resulting in four competing occupancy models: sampling period, sampling period +elevation, elevation and constant. We did not include quadratic terms because the elevation gradient covered is such that it was unlikely that a "bell-shaped" curve could be fitted. We built occupancy models in a Bayesian framework using the package "R2WinBUGS" in R 3.3.1 (Gelman, Sturtz, Ligges, Gorjanc, & Kerman, 2013;R Core Team, 2014). For each model (4 models for each of 123 species =492 models in total), we used three chains and sampled 10,000 iterations with a thinning rate of 8. For each model, we calculated Bayesian information criterion (BIC) scores (Ntzoufras, 2011). As all BIC scores of competing models for each species were close (differences usually within 10 units), we calculated BIC weights for each occupancy model to produce occupancy curves along elevational bands by model averaging the candidate set (Hooten & Hobbs, 2015). We then calculated 95% confidence intervals based on unconditional variance of averaged models according to Buckland, Burnham, and Augustin (1997), and assumed that non-overlapping confidence intervals within elevational bands identified significant changes in occurrence probability.

| Elevational limit shifts
In addition to comparing estimated occurrence probabilities for each species, we tested explicitly for shifts in range limits by calculating the cumulative probability of false absence where raw occupancy data suggest a range shift has occurred (Moritz et al., 2008;Tingley & Beissinger, 2009). This approach is robust in minimizing errors of reporting range shifts when no actual range shift has occurred (Tingley et al., 2012). To calculate the probability of false absence (P fa ) for an apparent range shift, we used the detection probability (p i ) for each species from the best occupancy model (lowest BIC value) and the formula as reported by Tingley and Beissinger (2009). We considered a P fa <0.05 as a significant range shift; that is, species were likely to have truly changed their occupancy status rather than an apparent shift arising from false absences.
As all the other metrics we compared so far in this study were drawn from species that occur in both time periods, we made a direct comparison of species checklists to check for potential changes in colonization or local extinction. During the entire study period (1994 to 2014), the moth sampling data were not only limited to our study site (Tai Mo Shan elevational gradient), but also included 5,358 species occurrence records spread across 182 sites in Hong Kong.
Using that list, we compared species of different time periods and listed 15 candidate species that only appeared in either the current or historic sampling period. We verified this list with the same P fa approach described above, but this time we ad hoc calculated detection probability for each species by dividing its total number of detections by the total number of visits to occupied sites, as we could not fit occupancy models for these species due to the haphazard nature of sampling. We calculated the P fa of sites with visits where we failed to record each species and also used visits from the geographically nearest sites for unpaired sites. Additionally, we compared those 15 species to records from iNaturalist, a citizen science program with wide participation in Hong Kong (https://www.inaturalist.org/projects/hong-kong-moths), to see whether any validated records (records with a photo associated that could be verified as to the species identification) exist within our study period but were not recorded in our own surveys.

| Community level elevation changes
We calculated "community elevation scores" for each shared site between time periods (n = 6: 70 m, 175 m, 210 m, 291 m, 307 m, 468 m) similar to Chen et al., (2009). This community metric treats the elevational centre of each species in the historic sampling period as a fixed "trait." By averaging these traits of species occurring at one site together, we produced the "community elevation score" for that site. A lower community elevation score in the current period indicates a higher proportion of species with lower elevational centres occurring in that site, which suggests a community-wide shift to higher elevations. Paired-Wilcoxon tests were performed to determine whether community changes across all sites were significant.
To calculate the elevational centre for each moth species, we first calculated the relative occurrence probability for each species by dividing occurrence in each site by the number of sampling times in that site. Then, we standardized these occurrence probabilities to one, and used them as weights to average the elevations of sites for each species.

| RE SULTS
Environmental data suggest that the climate has changed in Hong Kong over the sampling period, but the natural environment has mostly remained constant. For climate data, we detected no significant change in annual mean temperature, mean daily minimum temperature or annual mean precipitation from 1997 to 2014 (Supporting information Figures S2 and S3) When checked with P fa , 32 out of 123 moth species showed lower or upper range limit shifts (Supporting information  (Figure 3 and Supporting information Table   S3).
Of the 15 candidate species that only occurred in one time period, the P fa approach found that Alex palparia and Herochroma viridaria (Supporting information Table S2) were likely truly absent from historic recordings and were newly established in the current sampling period. Records from iNaturalist of these two species in the current sampling period (2010-2014) outside of the study site further confirm the recent colonization of Hong Kong by these two species. Based on iNaturalist records and our own surveys, Calletaera dentata may also be newly established-however, we cannot rule out the possibility of false absence due to limited historic sampling effort in our surveys relative to its detectability (Supporting information   Table S2). For the other 12 species, due to limited sampling intensity in our dataset, the P fa method could not rule out the possibility of false absence. However, iNaturalist records indicated that they were actually present elsewhere in Hong Kong across the whole sampling period. In this case, we do not know whether they went locally extinct or are new arrivals to the Tai Mo Shan area, but we know they were and are present in Hong Kong.
Community-wide, we found that the average elevational centres of species present at re-surveyed sites decreased over time at 5 out of 6 sites, indicating widespread upslope shifts (Wilcoxon signedrank test, p = 0.06). For these five sites, current community elevation scores were 3 to 27 m lower than historic scores, representing more

| D ISCUSS I ON
This study represents one of the first examinations of climate change effects along an elevational gradient in the lowland tropics ( Figure 1). Although moths are generally nocturnal, day-time temperatures can still affect resting individuals and larvae (Xing et al., 2018). With minimal habitat change but significant warming over the past decade (as reflected by a 0.5°C increase in maximum temperatures), we detected complex changes of geometrid moth species at both species and community levels. Increased averaged daily maximum temperatures, and increased numbers of extreme hot and cold days, result in higher fluctuations in daily temperatures. These fluctuations potentially are more challenging for organisms than pure increases of daily mean temperatures (Paaijmans et al., 2013;Vasseur et al., 2014). This might explain why, despite the relatively short resurvey timespan (15-year study duration), we detected significant changes in the moth communities, particularly in range limits and at community levels, despite using a highly conservative inferential approach. These changes not only indicate high sensitivity of moths to environmental change, but also provide possible early warning signs for larger climate change impacts in the future (e.g., 48 of 123 moth species showed temporal trends in occupancy). In other tropical montane sites, similar heterogeneous shifts have been reported in resurveys of birds (Forero-Medina et al., 2011;Freeman & Freeman, 2014), plants (Feeley et al., 2011;Jump et al., 2012), moths (Chen et al., 2009;Chen, Hill, Shiu, et al., 2011) and reptiles (Raxworthy et al., 2008). In addition, studies across different taxa in other continents and cooler climates have also found similar heterogeneous shift patterns (Forister et al., 2010;Rowe et al., 2015;Tingley et al., 2012 which could change without detection by the satellite image forest cover change criteria we used here) are known to lag warming impacts due to relatively low dispersal abilities (Rehm, 2014), while predatory birds might react differently and move faster or more slowly (Forero-Medina et al., 2011;Rehm, 2014). These interactions could indirectly influence moth presence to a greater extent than direct temperature or precipitation changes and produce complex shift patterns.
Additionally, the context of this study system might also affect the interpretation of these results. Community elevation scores results may only represent patterns of changed site occupancy but not "true" elevational shifts, as it only includes six sites. In addition, rare species are common in the tropics (Price, Diniz, Morais, & Marques, 1995) and are potentially more sensitive to climate change impacts (Iverson, Schwartz, & Prasad, 2004). However, such species also have low detection probabilities, which therefore require greater sampling effort to detect changes and shifts. Different species might also react to sampling methods differently and produce underestimated or overestimated changes in occupancy or abundance (e.g., Eumelea ludovicata is not easily attracted to light or light traps but is frequently observed in Hong Kong).
Finally, these results lead to recommendations for long-term monitoring programs in the tropics for the detection and management of climate-driven species redistribution (Basset et al., 2017;Bonebrake et al., 2018). First, sampling effort should be balanced between sites and sampling methods should be diversified to cover broad communities. Second, repeated sampling is crucial to effectively calculate detection probability and record abundance changes. Third, citizen science data (in this case iNaturalist records) can be used as a reliable source to verify or support findings from standardized monitoring efforts. Where possible, more population-level studies are also needed to examine species redistributions and the processes underlying these changes (Brown et al., 2016). And for those species that do show trends of range shifts and abundance changes, we need to prioritize and support the monitoring of their populations and distributions, in case intervention is required if further climate change endangers their persistence. Meanwhile, species with different life history traits will respond differentially to the same F I G U R E 4 Comparison of moth community elevation scores between the two time periods in Tai Mo Shan area, Hong Kong SAR. Except for one site at 70 m, all sites show that the community elevations score has decreased over time, defined by an increased proportion of species with low-elevation centres, suggesting an upslope shift of low-elevation distributed species change in climate (Garnas, 2018). Collecting additional natural history information on such species is important and will allow more sophisticated trait-based analyses which can uncover mechanistic processes behind range shifts (Dulle et al., 2016).

CO N FLI C T O F I NTE R E S T.
Authors declare no conflict of interest.

ACK N OWLED G EM ENTS
This work was funded by the Research Grants Council GRF (HKU 760213).

DATA ACCE SS I B I LIT Y
Data are available in the Supporting information (Figures S1-S5 and Tables S1-S4).