Landscape and Species Traits Co-Drive Roadkills of Bats in a Subtropical Island

: The expansion of roads has threatened wildlife populations by driving casualties due to vehicle collisions. However, the ecological drivers of wildlife roadkills are not yet fully explored. We investigated the strength of landscape features and ecomorphological traits in determining spatial patterns of bat roadkills in Taiwan. In total, 661 roadkills that belonged to 20 bat species were acquired by citizen scientists between 2011 and 2019. The number and species richness of victim bats declined with increasing elevations with varying species compositions. Elevation and artiﬁcial light had signiﬁcantly negative effects on the occurrence of roadkill, whereas protected area and its interaction with elevation had positive effects. Ordination analyses showed that roadkills were driven by different ecomorphological traits and landscape features. At low elevations, road casualties were associated with an aerial hawking hunting strategy. At higher elevations, roadkills were associated with higher elevational distribution. Roadkills of non-cave bats were associated with brighter environments, suggesting that bats might be exposed to higher risk when hunting insects near artiﬁcial light. Our ﬁndings suggest that management agencies shall consider both species traits and landscape features when planning impact assessments and mitigation practices of roadkills for bats and probably other wildlife, particularly when long environmental gradients are covered.


Introduction
The development of modern transportation is considered one of the major threats to global biodiversity in the Anthropocene. The increased use of automobiles and expansion of roadway infrastructures can impact terrestrial ecosystems via fragmentation and degradation of natural habitat, pollutions (light, noise, and air), disruption of energy flow, nutrient cycles, and water regeneration process, accidental introduction of alien species, as well as wildlife collision mortality [1,2]. Among all these threats, collision with vehicles might be the most immediate road-associated process that drives the direct loss of wildlife individuals, which could lead to declines and degradation of vulnerable populations [3][4][5] and damage efforts of existing conservation works. For example, road mortality of 124 vertebrate species, including 17 IUCN (International Union for Conservation of Nature and Natural Resources) Endangered and four Vulnerable species, is reported from a confined road section (47 km in length) of the Dong Phayayen-Khao Yai World Heritage Site in Thailand [6]. In extreme cases, subsequential loss of key ecosystem functions may occur if collisions suppress the keystone species due to population declines and road avoidance [7]. Therefore, understanding the influences and underlying drivers of wildlife roadkill is critical, especially since the expansion of roadway infrastructure and traffic volume increases continuously every year [8].
With over 1400 species [9], bats are one of the most species-rich and ecologically diverse mammalian orders. Bats are usually considered as good indicators of environmental changes due to their strong associations with habitats and high sensitivity to disturbances [10], including mortality due to vehicle collision [11]. Although bat roadkills are only reported from a few countries [6,[12][13][14][15][16], the fact that bats and viceless can be found globally suggests that it is an overlooked threat to global bat populations. Incidence rate and spatial patterns of bat roadkills have been demonstrated to be associated with roadside microhabitats, landscape features, for example, distance to critical resources, and roadway designs [6,14,[17][18][19]. Other factors, including volume and speed of traffic [19], sex, age [13,18], and seasonality of reproductive phenology, migration, and hibernation behaviors [6,14,18], have also been found to be related to bat roadkills in certain areas. Interspecific variation in roadkill rate within a bat assemblage is likely to be driven by differences in behavioral responses to road infrastructure, flight mode, and habitat selection mediated by the constraint of species traits. A meta-analysis based on 12 empirical studies by [13] concluded that low-flying bats are more vulnerable than high-flying bats. Other cases show that trophic level [15], abundance [17], and tolerance to human disturbance [19] are also important traits in mediating the roadkill of bats. In addition, evidence shows that bats adjusted their activity patterns near roads based on their sensitivity to adjacent habitat types [20], traffic noses [21], and streetlights [22][23][24]. Although not documented, such behavioral adjustments could directly affect the likelihood that bats collide with vehicles. In contrast to the success of environment variables and species traits in explaining the spatial patterns of bat roadkills at the local scale, such factors have not been found to explain spatial patterns at the regional and global scale. Vulnerable species and key environmental and trait determinants vary between studies, locations, and geographic regions [11,18]. Studies on Brazilian bats show that frugivores and omnivores (family Phyllostomidae) have more roadkill incidence than insectivore and carnivore species [15,19]. However, a significantly higher road mortality rate of insectivores compared to frugivores (family Pteropodidae) was observed in Thailand, contradicting the findings in the Neotropics. A similar case but less complicated scenario includes examples from Poland. Lesiński (2007) [18] reported bat roadkills were more common for the ground gleaners and rare for the high aerial hawkers. In contrast, Lesiński et al. (2011) [25] found that Nyctalus noctula, a common high aerial hawker in Poland, contributed almost near one-third of the total bat roadkills in the central region of the country, representing the most vulnerable species. The discrepancy in the findings among studies suggests that environmental features or trait diversity alone is not enough to explain the regional variations between areas with distinct bat fauna, which is the case when comparing patterns of bat roadkills between Brazil and Thailand, as mentioned earlier.
Taiwan Island houses a diverse bat fauna that contains 33 species of six families [26]. With more than 200 mountain peaks attaining heights over 3000 m on the island, the terrain, climatic conditions, vegetation diversity, and intensity of anthropogenic disturbances vary significantly along the elevational gradient in Taiwan [27][28][29]. Unsurprisingly, the island's bat species also show strong distribution patterns across the elevational gradient. Bat fauna of the higher elevation areas is characterized by the presence of some Eurasian taxa, such as genera Plecotus and Barbastella, whereas bat fauna of the lower elevation areas is characterized by the presence of Indochinese taxa, including Pteropodidae, Hipposideridae, Kerivoula, and Scotophilus [30,31]. The densities of road and vehicles in Taiwan were 118 km/km 2 and 204 unit/km 2 in 2017, respectively, which are considered high compared to global standards [32,33], and hundreds of incidences of bat roadkill have been reported throughout the island [34]. However, the information about roadkill of individual bat species and the ecological causes have not been described and evaluated. The considerable variations in abiotic, biotic factors, and the diverse bat fauna that occur along the elevational gradient present us with a unique opportunity to understand how the spatial patterns of roadkills, environmental variables, and species traits are associated with each other.
In the present study, we revisited the same dataset used by Chyn et al. (2019) [34] with recent updates to assess the associations between bat roadkill, landscape features, and ecomorphological traits in Taiwan. Our objectives were three-fold. First, we investigated the spatial and temporal patterns of bat roadkill incidence. Then, we identified the influences of different landscape features in driving the spatial pattern of bat roadkill occurrence. Lastly, we tested the responses of individual ecomorphological traits to the influential landscape features associated with the occurrence of bat roadkill in space. We predicted the following: (1) the occurrence of bat roadkills would covary with species-specific elevational distributions. More specifically, we expected Indochinese species to have a higher chance of road collisions at lower elevation, whereas Eurasian species will have a higher chance at higher elevation; (2) consequently, assemblage structures of bat roadkill fauna would also show elevational variations. Such variations would be mediated by trait responses that covaried with the elevational variations in environmental features. Specifically, we expected that forest bat species that were usually characterized by a high-pitched echolocation call and clutter foraging strategy would have higher chances of road collisions at higher elevation where forest cover is generally high. On the other hand, species adapted to a disturbed landscape that were usually characterized by low-pitched calls and an open space foraging strategy would have higher chances at lowland where forest cover is low.

Materials and Methods
The study area is located on Taiwan Island, which is centered on the Tropic of Cancer in East Asia (21 • 55 -25 • 20 N, 119 • 30 -122 • 00 E). The island is 36,193 km 2 in area. In the central and eastern parts of the island, the terrain is mountainous and has an altitude of up to almost 4000 m asl. The western part is typically composed of plains and is hilly, which are usually less than 100 m in elevation. The island has a typical marine subtropical climate with an average annual precipitation of over 2500 mm. Typhoons often visit the island at least twice per year during July to September. Cold monsoon weather originates from the northeast to the island during the winter (November to early March). The temperature on the island exhibits a strong elevational pattern, which decreases significantly from the low to high elevation. The average temperature of the coldest and hottest months ranges from 16 to 20 • C and 28 to 29 • C, respectively, in the lowland, and it ranges from 0 to 6.5 • C and 8 to 15 • C, respectively, at the elevation over 3000 m. During the winter, the mountains peaks are often covered with snow in January and February [35].
In the lowlands and foothills of the mountains, particularly in northern and western Taiwan, the landscape is highly populated and dominated by agricultural lands, secondary forest fragments, and urban habitats, including high density of transportation infrastructure [29,36,37]. The mountainous areas are generally forested, except where the altitude is over 3200 m asl and is dominated by open habitats, such as short bamboo stands and bare land. The forest types vary along the elevational gradient. At low elevations (usually <1400 m asl), the natural forest is typically dominated by evergreen board-leaved trees of Fagaceae, Lauraceae, and Theaceae, and it is often mixed with large bamboo stands. At the mid-elevation level (c. 1400-2600 m asl), the natural vegetation is composed of evergreen and deciduous broad-leaved forest and mixed cloud forest dominated by trees of Fagaceae, Lauraceae, and Theaceae. The major forest type at the higher elevation level of 2500-3400 m is alpine forest dominated by coniferous trees of Cupressaceae and Pinaceae. At elevations above 3200 m, the vegetation is gradually replaced by conifer woodland and scrubland, including Fargesia bamboo stands, with increasing altitude and, hence, a decreasing dominance of canopy trees [28]. Old artificial forest of Japanese Cedar (Cryptomeria japonica), various types of orchards, or tea and vegetable plantations are also common at elevations up to 2600 m.
Bat roadkill data were extracted from the Taiwan Roadkill Observation Network (hereafter TaiRON [38]), which is a database designed by the Endemic Species Research Institute to collect roadkill incidence data of vertebrates and terrestrial crabs since 2011 [34]. For the present study, we restricted the records to those collected by citizen scientists, of which most records were attached to images of animals taken in the field and allowed us to confirm the initial identification by the observers. In some cases, animal corpses in good condition were taken as evidence by the observers and later preserved as specimens at the zoological museum of the Endemic Species Research Institute in Jiji Township, Taiwan. Putative taxa of the bats in photos were first identified by the collectors in the field and then by volunteers online in TaiRON's Facebook group. The results were further revised by examining both the images and specimens by the first author. Samples were identified to the finest taxonomic level by morphological and craniodental traits following the existing keys in the literature [26,31].
For the species level, we summarized the elevational patterns of the occurrences of roadkills of individual bat species and compared the results with the elevational range of the known distribution for each species. To understand how bat roadkill "assemblages" varied along elevation, samples were assigned to different elevation subregions following a classification suggested by Chang (2015) [30] for bats in Taiwan. The classification subdivides Taiwanese bat fauna into five categories based on species elevation limits, which are 0-100 m, 100-500 m, 500-1500 m, 1500-2500 m, and >2500 m. Due to the small number of sampling grids at high elevation, we combined the highest two subregions in Chang (2015) [30] into one category, >1500 m. The new classification granted a more balanced sampling among the four elevational subregions (10,943, 8954, 9404, and 7670 grids for the subregions I, II, III, and IV, respectively) (also see Supplementary File S1 for the statistical justification). We first calculated the number of grids with bat roadkill detected and species richness of victim bats among elevation subregions. Then, we tested if the proportions of grid numbers and species richness differed from the expectations based on the ratios of the overall roadkill sample grids and known species richness across subregions using a Chi-square test of equality of proportions. The beta diversity of the bat roadkill assemblage between neighboring subregions was estimated to further evaluate the degree of intersubregion variations in species compositions. We partitioned the overall beta diversity into species turnover (β SIM , replacement of species) and nestedness (β SNE , loss of species without replacement) components to account for the contributions of independent and dependent variables respectively, on richness differences between subregions [39,40].
To assess the associations between bat roadkill and landscape features, we excluded the roadkill data collected prior to the TaiRON citizen science standard procedure officially launched in 2013 to minimize sampling bias. All records were mapped to the grid system by ArcGIS software version 10.6, which subdivided the land area of Taiwan into 40,254 1 km-by-1 km girds according to the georeferenced data attached to the evidence images. Grids without any record of vertebrate roadkill were excluded from the analyses to avoid false absence due to low detectability or lack of survey effort. Due to the difficulty of distinguishing cryptic species of corpses morphologically in most cases, a hybrid of species and morphospecies was used for the analyses of trait responses to landscape features in association with bat roadkills. Two morphospecies are used for the following two species groups, including "Myotis frater complex" for M. frater and M. soror, and "other Pipistrellus" for all known Pipistrellus species, except P. abramus. This grouping arrangement only potentially reduced a small proportion of the victim "species" richness (up to 3) but gave opportunities to increase samples of the rare species in the sample due to the constraints of identification.
Twenty landscape features were used to test the associations between bat roadkills and environmental variables. Three topographic variables were selected, including the mean elevation, the mean slope, and distance to the nearest freshwater body. Ten landcover variables were selected, namely areas of farm field, meadow, forest, farm wetland, bushland, wetland, urban area, waterbody, bare land, and buildings. Habitat heterogeneity, which was estimated as the Shannon's Diversity Index based on all habitat types, except the area composed of buildings, is included. The area composed of buildings is not included because it is acquired separately from the other habitat categories. Six anthropogenic landscape measures, including night light intensity, total length of roadway, number of human populations, occurrence of protected area, and distance to the nearest protected area and a fragmentation index are used. Due the lack of accessible traffic data, the number and density of the human population were used as approximates of the number and density of vehicles. Data of all variables were acquired from open sources [41,42] except for the following eight variables: habitat heterogeneity, night artificial light intensity, human population size, human population density, occurrence and extent of protected area, distance to the nearest protected area, and the fragmentation index, which were generated in this study. All variables were originally digitalized at the same spatial resolution of 1 kmby-1 km as the roadkill data, or, if not, rescaled to the same grid system by extrapolation and interpolation methods using the weight mean algorithm. Detailed definitions of the landscape variables are given in Appendix A Table A1, and the measurements are provided in the online Supplement Table S1.
Ten ecomorphological traits were selected to test the associations between traits, landscape features, and roadkills, seven of which are quantitative measures, namely forearm length, body mass, wingspan, frequency of maximum energy in echolocation calls, minimum of elevational distribution, maximum elevational distribution, and range of elevational distribution. For the first four quantitative measures, median values were estimated by the maximum and minimum values across literatures, as most available information are given in range or a single value. The remaining three traits are qualitative, namely echolocation call type, roost type use, and hunting mode. For the call type, we primarily followed the classification suggested by Jones and Teeling (2006) [43] but with an additional call type later described from Coelops frithii, which is one of the collision victim species in this study, by Ho et al. (2013) [44]. We followed Denzinger and Schnitzler (2013) [45] to assign species into five different hunting guilds according to foraging strategy, namely open space aerial-hawking, edge and gap aerial hawking, water trawling, and gleaning. Lastly, the use of cave, rock crevice, tree cavity, or foliage roosts is presented as presence-absence separately. Despite diverse roost types described in the literature, roost use category assignment was based on the shared roost attributes, such as physical structure and surface material, regardless of whether they are natural or artificial manmade structures. For example, the use of bat boxes was assigned to the category of tree cavity based on the confined space surrounded by woody materials. Trait data were summarized from existing information in the literature (see detailed citations information in the online Supplementary Table S2), the authors' unpublished data from empirical studies on wild individuals and museum collections, previous field observations, and elevation data from the present project. To reduce measurement bias by different researchers, e.g., from differences in the minimum units among measurement kits used, all quantitative measures were rounded to the first decimal place. Trait data of each morphospecies were combined with the measurements of all species within each group. Definitions of the selected traits are listed in Table 1.  [31,44] All statistics were carried out using R version 3.5.3 under RStudio environment (desktop version 1.1.456). Chi-square tests of equality of proportion were deployed using the prop.test function in the default stats package. Occurrence-based beta diversity measures were estimated with the Sorensen Index using the beta.pair function in the betapart package v.1.5.2. [57] for the following three elevation subregion pairs: subregions I vs. II, II vs. III, and III vs. IV. We used generalized linear mixed models (GLMM) with log link to test the significance of the effects of 20 landscape features (fixed factors) on the occurrence of bat roadkills in each grid with controls of the influences of event year (random factor) and month (fixed factor) using glmer function in the lme4 v. 1.1.23 package [58]. Quantitative measures were rescaled to z-score to avoid failure in estimating variances due to extreme values in some variables. We first ran the model with all 20 variables and the subsequential models based on selected variables with statistical significance in the previous model and their interactions. Then, we performed Akaike information criterion (AIC) and Bayesian information criterion (BIC) to choose the best predicted subset among the selected models based on the scores of both fit criteria using Ictab function in the bbmle v. 1.0.20 package [59]. To avoid inaccuracy in parameter estimates due to multicollinearity, only the selected landscape features of the best fitting GLMM model were used in the analyses of trait responses. RLQ and fourth corner analyses were used to assess trait-based responses to landscape features. Both methods are ordination-based correlations of three tables, namely L matrix (species by grids), R matrix (landscape features by grids), and Q matrix (ecomorphological traits by species). The RLQ analysis measures the overall correlations between the three matrices. The relationship between the traits and landscape features can be tested by a global permutation test. The fourth corner analysis measures the strength of each individual trait-landscape feature relationship and tests the significance using a permutation procedure [60]. The RLQ and fourth corner analyses were performed using the rql and fourthcorner functions, respectively, in ade4 v. 1.7.25 package [61]. The significance of the global structure of RLQ analysis was tested using the randtest function. Permutations tests of the outputs of both methods were run with 14,999 iterations using the model 6 procedure to avoid inflated type I error [62]. In addition, the fdr method was used for estimating the adjusted p-values of the permutation test for the fourth corner model to avoid false discoveries [63]. In the present study, both RLQ and fourth corner analyses were built based on occurrence-based community data.

General Patterns of Bat Roadkill Incidence
A total of 661 bat roadkill incidences were collected by the TaiRON from 2011 to March 2019. Of these, 62.2%, 4.5%, and 15.6% of the samples were identified to species, genus, and family levels, respectively ( Table 2). Twenty species belonging to 14 genera of four families were recognized, 12 of which are endemic to the island at either species or subspecies levels. Two species, namely Coelops frithii formosanus and Myotis formosus flavus, were listed as Nationally Vulnerable in Taiwan and Near Threatened by the IUCN Red List. The remaining 18 species were either not considered as threatened both nationally and globally, or lack of sufficient data for further assessment for their conservation status ( Table 2). With 98, 75, and 68 records respectively, Pipistrellus arbramus, Miniopterus fuliginosus, and Hipposideors armiger terasensis were the top three roadkill victim species, together contributing 58.2% of the samples with species identified. Scotophilus kuhlii, Murina puta, Rhinolophus monoceros, and Eptesicus pachyomus horikawai were also common in the samples, contributing to 6-10% of the total samples. The remaining 12 species were rare in the samples, which contained less than 10 records of each species and together accounted for only around 9.2% of the total samples with species identified (Table 2). Bat roadkills (n = 598) also varied between months during the year 2013-2018. June to August, which was the summer season, represented the major peak of bat roadkills in all surveyed years, except in 2013, in which more bat roadkills were reported between August and October. September-October and March-April represent two minor peaks of bat roadkills within a year, and the remaining months had few records (Figure 1). Bat roadkill incidences were found throughout the island in general at a wide elevation range from the lowland to 3152 m (   [64].

Elevational Patterns of Occurrences of Bat Roadkill at Species Level
There were seven morphospecies and 17 morphospecies out of the 19 bat morphospecies found with roadkill records at elevations lower than 100 m and 500 m, respectively (

Elevational Patterns of Occurrences of Bat Roadkill at Species Level
There were seven morphospecies and 17 morphospecies out of the 19 bat morphospecies found with roadkill records at elevations lower than 100 m and 500 m, respectively ( Figure 3). Roadkills of seven and four species were recorded at elevations over 1000 m and 2000 m, respectively. Submyotodon latirostris is the only species recorded at over 3000 m and the only species with an elevation lower limit over 1000 m. Four species, namely H. armiger, R. monoceros, M. puta, and M. fuliginosus, showed a wide elevational range of roadkill distribution between 1500 and 2000 m. The remaining 15 species showed narrower elevational ranges (<1000 m) of roadkill distribution, although 13 of these species had a much wider elevational range in the species distribution (1578-3165 m) (Figure 3).

Elevational Patterns of Bat Roadkill at Assemblage Level
The elevation subregion I had the highest number of grids with bat roadkill recorded (n = 250), following by the subregion II (n = 144), subregion III (n = 39), and subregion IV (n = 9). The proportions of grid number of bat roadkill to the grid number of overall vertebrate roadkills of each elevational subregion differed significantly among the four subregions (X 2 = 248.1, d.f. = 3, p < 0.001, test for equality of proportions). With 16 species, the elevation subregion I had the highest species richness of bat roadkill victims. Species richness of the subsequential subregions decreased with increasing elevation. The subregions II, III, and IV had 15, 9, and 6 species recorded. The proportions of victim species richness to the known species richness of each elevational subregion (17,17,16, and 17 species for the subregions I, II, III, and IV, respectively) differed significantly among the four subregions (X 2 = 18.1, d.f. = 3, p < 0.001, test for equality of proportions). The overall measures of beta diversity increased with increasing elevation (Sorensen measures increased from 0.16 to 0.47 from the subregion pair I vs. II to the pair III vs. IV), suggesting increasing differentiations in species compositions between neighboring subregions at higher elevations ( Figure 4). Nestedness contributed higher proportions than species turnover to the beta diversity of the two mid-elevation subregions (β SIM : 33.3%, β SNE : 66.7%). In contrast, species turnover processes contributed higher proportions than nestedness to the beta diversity of the pairs of the subregions I vs. II (β SIM : 81.3%, β SNE : 18.8%) and the subregions III vs. IV β SIM : 70.2%, β SNE : 27.8%).

Elevational Patterns of Occurrences of Bat Roadkill at Species Level
There were seven morphospecies and 17 morphospecies out of the 19 bat morphospecies found with roadkill records at elevations lower than 100 m and 500 m, respectively ( Figure 3). Roadkills of seven and four species were recorded at elevations over 1000 m and 2000 m, respectively. Submyotodon latirostris is the only species recorded at over 3000 m and the only species with an elevation lower limit over 1000 m. Four species, namely H. armiger, R. monoceros, M. puta, and M. fuliginosus, showed a wide elevational range of roadkill distribution between 1500 and 2000 m. The remaining 15 species showed narrower elevational ranges (<1000 m) of roadkill distribution, although 13 of these species had a much wider elevational range in the species distribution (1578-3165 m) (Figure 3).

Elevational Patterns of Bat Roadkill at Assemblage Level
The elevation subregion I had the highest number of grids with bat roadkill recorded (n = 250), following by the subregion II (n = 144), subregion III (n = 39), and subregion IV (n = 9). The proportions of grid number of bat roadkill to the grid number of overall vertebrate roadkills of each elevational subregion differed significantly among the four subregions (X 2 = 248.1, d.f. = 3, p < 0.001, test for equality of proportions). With 16 species, the elevation subregion I had the highest species richness of bat roadkill victims. Species richness of the subsequential subregions decreased with increasing elevation. The subregions II, III, and IV had 15, 9, and 6 species recorded. The proportions of victim species richness to the known species richness of each elevational subregion (17,17,16, and 17 species for the subregions I, II, III, and IV, respectively) differed significantly among the four subregions (X 2 = 18.1, d.f. = 3, p < 0.001, test for equality of proportions). The overall measures of beta diversity increased with increasing elevation (Sorensen measures increased from 0.16 to 0.47 from the subregion pair I vs. II to the pair III vs. IV), suggesting increasing differentiations in species compositions between neighboring subregions at higher elevations ( Figure 4). Nestedness contributed higher proportions than species turnover to the beta diversity of the two mid-elevation subregions (βSIM: 33.3%, βSNE: 66.7%). In contrast, species turnover processes contributed higher proportions than nestedness to the beta diversity of the pairs of the subregions I vs. II (βSIM: 81.3%, βSNE: 18.8%) and the subregions III vs. IV βSIM: 70.2%, βSNE: 27.8%).

Associations between the Occurrence of Bat Roadkills and Landscape Features
Of the 20 landscape features, 12 variables, namely the mean elevation (ELE), artificial light intensity (ALight), presence of protected area (PA), and all types of landcover, except the area of building, were significant in the full GLMM model (Table 3). ELE and ALight had negative effects on the occurrence of bat roadkills, whereas PA and the nine landcover types had a positive effect. A second GLMM model based on the 12 variables suggests that the effects of all landcover variables were not significant. Therefore, they were excluded from the subsequential models. Then, we tested the effects of the remaining three variables and the interactions between variables on bat roadkills. The results of the six

Associations between the Occurrence of Bat Roadkills and Landscape Features
Of the 20 landscape features, 12 variables, namely the mean elevation (ELE), artificial light intensity (ALight), presence of protected area (PA), and all types of landcover, except the area of building, were significant in the full GLMM model (Table 3). ELE and ALight had negative effects on the occurrence of bat roadkills, whereas PA and the nine landcover types had a positive effect. A second GLMM model based on the 12 variables suggests that the effects of all landcover variables were not significant. Therefore, they were excluded from the subsequential models. Then, we tested the effects of the remaining three variables and the interactions between variables on bat roadkills. The results of the six models based on subsets were similar. The three selected landscape features had similar effects on bat roadkills to the grand model, and their effects were all significant. The interaction between PA and elevation had a significant positive effect, suggesting that the influence of protected area on bat roadkills increases with increasing elevation or vice versa. The interaction between ALight and PA had a significant negative effect, suggesting that the influence of protected area on bat roadkills decreases with increasing night light intensity or vice versa. No significant effects were found in the interaction between ELE and ALight as well as the interaction among three variables (Table 3). Combining the results of AIC and BIC criteria, the GLMM model contains ELE, ALight, PA, and the interaction between ELE and PA was the best fitting model ( Table 3). Table 3. Evaluation results of seven generalized linear mixed models (GLMM) models with selected landscape features using Akaike information criterion (AIC) and Bayesian information criterion (BIC) criteria. In model formula, a cross sign denotes that an interaction of the prior and posterior variables is considered, whereas a plus sign denotes that interaction is not considered. Variable abbreviations-ELE, ALight, PA, BU, and BL stand for the mean elevation, artificial light intensity, presence of protected area, bushland, and bare land. In AIC and BIC scores, numbers in bold indicate the two smallest values among all results in either criterion.

Trait-Based Responses to Landscape Features in Terms of the Occurrence of Bat Roadkill
The overall association between the ecomorphological traits and the three GLMMselected landscape features in the RLQ analysis was significant (permutation test: p = 0.0018). The first two axes of the RLQ analysis explained approximately 92.6% and 6.3% of the overall variation, respectively. Among the 75 possible trait-environment associations, four were significant ( Figure 5). Bats that exhibited open space aerial hunting modes were negatively correlated with the mean elevation, whereas the lower limit of elevational distribution is positively correlated with the mean elevation. The use of cave roost was negatively correlated with night light intensity and vice versa (fdr adjusted p < 0.05, permutation test).

PA
0.768 ± 0.318 0.016 1 The temporal variables (year and month) and the dependent variable (occurrence of bat roadkill) are not presented. The general formula for all models is Bat roadkill occurrence = "selected landscape features" + Month + (1|Year). 2 The results of the two temporal variables and variables with no significance detected are not shown. 3 The seven landcover types include farm field (FF), meadow (MD), forest (FO), farm wetland (FW), bushland (BU), wetland (WT), urban area (UB), waterbody (WB), and bare land (BL).

Trait-Based Responses to Landscape Features in Terms of the Occurrence of Bat Roadkill
The overall association between the ecomorphological traits and the three GLMMselected landscape features in the RLQ analysis was significant (permutation test: p = 0.0018). The first two axes of the RLQ analysis explained approximately 92.6% and 6.3% of the overall variation, respectively. Among the 75 possible trait-environment associations, four were significant ( Figure 5). Bats that exhibited open space aerial hunting modes were negatively correlated with the mean elevation, whereas the lower limit of elevational distribution is positively correlated with the mean elevation. The use of cave roost was negatively correlated with night light intensity and vice versa (fdr adjusted p < 0.05, permutation test).

Discussion
In this study, we found significant elevational variations of bat roadkills from our long-term data in Taiwan. The overall abundance and species richness of collision victims declined significantly from low to high elevations, whereas evenness increased with increas-ing elevations. On the other hand, the beta diversity increased with the increasing elevation, suggesting increasing changes of species compositions of bat roadkills. The variations in species compositions between the mid-elevation areas with the lowland and high elevation were contributed more by species replacement than nestedness processes, whereas the variations within the mid-elevation areas show reversed patterns. At lower elevations, most of the roadkill victim species are Eurasian genera in the subfamily Vespertilioninae and family Miniopteridae. In contrast, at higher elevation, many of the road victims are species in the Indochinese families Hipposideridae, Rhinolophidae, and two vespertilionid subfamilies Kerivoulinae and Murininae. The results contradict our first prediction even though the elevational variations in incidences of roadkills were partially mediated by the constraints of species distribution-the lower elevation limit in this case rather than the upper limit as we predicted. Nevertheless, our findings support the local-assemblagestructure hypothesis that has been proposed in explaining the geographic variations of bat roadkills in Poland [18]. To our knowledge, this is the first study that demonstrates that elevation can affect spatial patterns of bat roadkills, revealing an overlooked but potentially important element in previous studies on wildlife road ecology.
Past studies have demonstrated the effectiveness of landscape features [6,13,14] and species traits [13,15,18,21] in explaining roadkill patterns of bats, respectively. Our data clearly show trait-environment associations in mediating bat roadkills in Taiwan. Beyond elevation and distribution limit, we also found that open-space bats are more likely killed at lower elevations outside protected areas. These species, which include Pipistrellus abramus, Miniopterus fuliginosus, Scotophilus kuhlii, and Eptesicus pachyomus, are commonly detected in disturbed lowland habitats, such as agricultural land, suburban, or urban habitats during acoustic surveys in the lowland [65], where the density of road and vehicles is generally high. It is noteworthy that all these species are common and widely distributed in East Asia or Southeast Asia [66], but only S. kuhlii had roadkill records reported from Thailand ( [6], in the Supplementary File S1). Further investigations of whether these bat species are also roadkilled in other regions, within the species distributional ranges, are needed to have a better understanding of the impacts of vehicle collision to their global populations.
Our second prediction hypothesized that there were roadkill-associated trait-mediated responses to different environmental features along elevation. We found that bat roadkills at higher elevation were associated with the distributional limit of individual species. Common victims include species of vespertilionid subfamilies Murininae, Kerivoulinae, Myotinae, genera Hipposideros, and Rhinolophus. Most of them are common and abundant during trapping sessions in mountain areas [51,67], foraging at understory and ground levels [68], and are seldom recorded from the lowland [30]. The decrease of open space aerial hawkers in the roadkills at higher elevation is likely due to the restriction of the elevational upper limit of species distribution in this case. For example, P. abramus, the most common victim species in this study, is the most abundant bat in highly disturbed landscapes in the lowland and hills [65,69] but is generally less common in the mountain areas (e.g., [67]), despite its distribution up to 1000 m [70]. Other open space bats, e.g., E. pachymous and Nyctalus plancyi, are also reported as common during acoustic surveys in mountain areas with elevations over 2000 m [51]. In contrast, these two species are rarely recorded during understory-level trappings in mountain areas (e.g., [51]). It is possible that these bats fly higher in mountain areas to avoid difficult terrains [71] or to follow insects at a higher elevation in the sky [72,73] as observed in other open space species.
We found a positive effect of protected area on the occurrence of bat roadkills, and the influence increased with increasing elevation. Although not significant, two traits, namely frequency of maximum energy and open-space hunting strategy, showed trends in associating bat roadkills with the presence of protected area (0.05 < fdr adjusted p < 0.01, permutation tests, this study). The effects of protected area and its interaction with elevation have general associations with non-open space bats using constant frequency (Rhinolophus and Hipposideros) or broadband echolocation calls with high frequency (Murininae and Kerivoula). The higher casualty occurrence rate in the protected areas might be associated with the higher-quality habitats within, e.g., old-growth forest. Unexpectedly, we did not find any significant associations between bat roadkills and habitat types. The lack of bat roadkill-habitat relationship in this study might be a result of failing to test habitat effects at appropriated spatial scales. Few studies show that the effects of environment variables on wildlife roadkills are scale dependent [6,21]. The effects of protected area on bat roadkills could be explained by tourist traffic as many of the protected areas in Taiwan are popular tourism sites and connected with some major roadways in the country [74]. Therefore, higher wildlife roadkills at areas with high visitation rates and traffic volumes are expected. In the present study, we use the number and density of human population as approximate measures of traffic, which are unlikely representative for tourism activity. Unfortunately, there are no open-source data of traffic at the national scale available for further analyses.
Artificial light is identified as an emerging anthropogenic impact on bats [75] and could be associated with bat roadkills in several ways. Bats with negative phototactic behavior may avoid areas with strong artificial light during the night [24,76], including bat-friendly underpasses, overpasses, and gantries [11,22]. The light avoidance behavior can force bats to fly across roads in darker areas to reach critical resources [77,78], which may increase the risk of bat-vehicle collision. This might explain why bat roadkills show a negative correlation between the cave roosting behavior and night light. Cave specialists, e.g., Rhinolophus bats, are known to be sensitive to light conditions near roost entrances [79]. Coelops frithii and Myotis formosus, the two most vulnerable species in our study, are known to roost in caves and show light-avoidance behaviors near the roost entrance and during foraging [80,81]. The anti-phototactic behavior may bring these bats closer to potential risks of collision with vehicles when they forage in areas or roost in caves close to roads in dark environments. On the other hand, bats that exhibit positive phototaxis might be attracted by swarming insects surrounding lights [82,83]. For example, P. arbramus, a non-cave roosting bat and the most frequently recorded victim species in this study, is often observed hunting insects near streetlights (JCCH personal observation). The food resourcerelated positive phototaxis could also bring this non-cave bat closer to greater risks of collision with vehicles, which may explain why we detected a correlation between nonuse of cave and night lights. Future studies shall focus on the elucidations of how taxa-specific responses of collision victim bats, particularly species with distinct trait-environmentroadkill relationships, to different models and settings of artificial light. Such knowledge is critical to the implantation of bat-friendly light into future road designs to minimize the impact of roadkills to bat populations and their associated ecosystem functions [76,77,82].
There are several limits of the present study. First, the roadkill data presented in this study were collected occasionally by volunteers and citizen scientists. The unknown sampling efforts, lack of a standard survey protocol, and unbalanced sampling constrain the usage of the incidence data and the interpretation of the spatial variation of bat roadkill patterns. Lack of unbiased estimation of roadkill incidence would also lead to the missed opportunity to evaluate how road mortality may impact populations of different bat species, particularly for those with high conservation concerns. It is also suggested that analyses based on different spatial scales could lead to different results when interpreting the influences of environmental features to roadkill patterns (e.g., [6]). Further investigation into the associations between bat roadkills and landscape features at different spatial scales at the hotspot identified in this study is recommended. Moreover, several important determinants of wildlife roadkill mortality, i.e., road designs [6,14,[17][18][19], vehicle speed and volume [19], and species traits, including abundance [18] and wing morphology measures [15,84], are not included due to lack of accessible data. Such information would be useful for a more comprehensive understanding of the potential underlying processes of bat roadkill on the island. It is noteworthy that information of important traits identified in this study, i.e., echolocation call features, foraging mode, roost ecology, and elevational distribution, for bats are generally deficient in the global mammal trait databases [85,86]. Therefore, filling the knowledge gaps of trait data is crucial for a better understanding of the impacts of roadkill and other road-related disturbances to global bat conservation.

Conclusions
In summary, the present study demonstrates that combinations of species traits and environmental variables can effectively explain the spatial patterns of bat roadkills along the elevational gradient in Taiwan. As conservation biologists and management agencies have started to pay more attention to the impacts of development of road systems to biodiversity in the recent years, it is reasonable to expect that data on wildlife roadkills would be generated quickly from different areas in the near future. The needs of the "fourth corner", the trait-environment associations, in understanding causes and processes of bat roadkills might be more relevant when elucidating data collected from large spatial extents, along long environmental gradients, or more functionally diverse taxa. The large-scaled sampling designs may produce results of mixed effects from divergent trait-environment relationships by including mega-species diversity. Caution shall be taken to avoid potential false conclusions of assemblage-level reactions if cause-effect analyses are drawn based on either the environment or the trait dimension only. We recommend that management and conservation agencies shall consider both species traits and environmental factors when planning impact assessments of vehicle collisions for wildlife populations and designing mitigation practices, i.e., bat-friendly infrastructures, especially when long environmental gradients are covered.
Supplementary Materials: The following will be available online at https://www.mdpi.com/14 24-2818/13/3/117/s1, File S1: Comparisons of roadkill samples of bats and non-bat vertebrates across elevational subregions in Taiwan from 2013 to 2018; Table S1: Measurements of 20 landscape features in 1 km-by-1 km resolution of Taiwan; Table S2: Measurements of ten species traits of 19 bat morphospecies of the roadkill victims in Taiwan.

Data Availability Statement:
The vertebrate roadkill data in this study is archived in Taiwan Roadkill Observation Network (https://roadkill.tw/en (accessed on 19 January 2021)) and are available on request from the corresponding author. The data are not publicly available because it is in a preparation for a dataset publication to Global Biodiversity Information Facility (https://www.gbif.org/ (accessed on 19 January 2021)). The environment data are available in the Supplement Table S1 and https://data.depositar.io/en/dataset/taiwan_environmental_dataset (accessed on 19 January 2021). The trait data are available in the Supplement Table S2.

Acknowledgments:
The authors would like to thank Tuanjit Srithongchuay and Guo-Jing Weng for the statistics support, Ying-Yi Ho for sharing call recordings of Murina bicolor, National Museum of Natural Science and the Mammal Division of the Endemic Species Research Institute for the permissions to access their collections of bats.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results. LAND COVER [41] Farm field, FF Continuous m 2 Area of rainfed cropland [41] Meadow, MD Continuous m 2 Area of herbaceous cover [41] Forest, FO Continuous m 2 Area of forested land [41] Farm wetland, FW Continuous m 2 Area of irrigated and post-flooding farmlands [41] Bushland, BU Continuous m 2 Area of shrubland [41] Wetland, WL Continuous m 2 Area of flooded land with fresh, saline, and brackish water [41] Urban Area, UB Continuous m 2 Area of artificial surfaces [41] Waterbody, WB Continuous m 2 Area of freshwater and saltwater bodies [41] Bare Land, BL Continuous m 2 Area of land without any landcover [41] Building, BD Continuous m 2 Area of buildings regardless usage [42] Habitat heterogeneity, H Habitat Index value Heterogeneity of land cover composition estimated as Shannon's Diversity Index. BD is not included in the estimation due its alternative source from other land cover variables.  Effective mesh density (seff)-a measure of to which movement within a landscape is interrupted by transportation infrastructure and geographic barriers. See more details in [91] This study