An assessment of the niche centroid hypothesis: Pteropus lylei (Chiroptera)

. Recent ecological theories propose that species reach their highest abundance and genetic diversity in the center of their ecological niche and decline toward the edges. We assessed whether Lyle's ﬂ ying fox, Pteropus lylei , abundance and genetic diversity were correlated with niche centroid distance using an ecological niche model as a proxy for fundamental niche ( N F ). Alternatively, we assessed whether P. lylei abundance and genetic diversity were correlated with ﬁ ne-scale environmental factors as a proxy of the species ’ realized niche ( N R ). We examined relationships between abundance and environmental factors at coarse and ﬁ ne scales as proxies of N F and N R , respectively. For coarse scale, ecological niche of P. lylei was modeled using all available occurrence records in Thailand, Cambodia, and Vietnam, coupled with climatic data. We conducted ﬁ eld measurements of P. lylei abundance and used genetic structure data across a large portion of the species ’ range. We measured Euclidean distances between abundance and genetic data and the niche centroid in environmental dimensions. Spearman ’ s correlation was estimated between abundance and genetic diversity vs. distance to the niche centroid. Complementarily, for the ﬁ ne-scale test, we measured multiple regression models between abundance and genetic diversity versus the normalized difference vegetation index (NDVI), local temperature, percent area of waterbodies, human density, and number of Buddhist temples. We failed to detect relationships between abundance and genetic diversity with the distance to the niche centroid in the coarse-scale model. When using the ﬁ ne-scale, landscape-level data, we found negative correlation between genetic diversity and number of temples. The data available were unable to support niche centroid hypothesis for the current distribution and abundances of P. lylei . We note that our failure to ﬁ nd an association does not support nor reject the niche centroid hypothesis. Instead, our capacity to test the niche centroid hypothesis may be limited by our ability to use empirical data to accurately reconstruct N F from ﬁ eld observations only. Future research may require physiology-based experimental approaches to explore relationships between species abundances and the niche structure.


INTRODUCTION
Understanding geographical patterns of species abundances is one of the most important challenges in ecology. Knowledge of the geography of abundances can help elucidate the finescale causes of species range limits, gene flow within populations, and population dynamics. At a coarse scale, the geography of abundances can help explain macroecological patterns such as species abundance distributions, species-area relationships, similarity-distance relationships, and fraction of clumped species . The abundant center hypothesis states that population abundance will peak at the center of a species' geographic range and decline toward the periphery (Hengeveld andHaeck 1982, Brown 1984). An alternative hypothesis is the niche centroid hypothesis formalized recently by Osorio-Olvera et al. (2019) and an extension of Hutchinson's ecological niche concept (Hutchinson 1957). The niche centroid hypothesis proposes a negative association between species' abundance and the distance to the niche centroid (Borregaard and Rahbek 2010, Mart ınez-Meyer et al. 2013, Manthey et al. 2015, Weber et al. 2017, Dallas et al. 2017. The niche centroid hypothesis proposes that species abundances should be measured in environmental dimensions rather than in geographic space. Based on the niche centroid hypothesis, higher production of crops, higher pathogen transmission, and higher genetic diversity are expected in the core areas of a species' ecological niche (Holt 2009, Mart ınez-Meyer et al. 2013, Lira-Noriega and Manthey 2014, Yañez-Arenas et al. 2014, Manthey et al. 2015, Ram ırez-Gil et al. 2018. For example, Lira-Noriega and Manthey (2014) found inverse relationships between the genetic diversity of a species and the distance to their niche centroid. This finding aligns with patterns of species' genetic diversity that are generally consistent with the niche centroid hypothesis (Dixon et al. 2013). Thus, the niche centroid hypothesis is strongly dependent on the species' ecological niche (Brown 1984). The niche, however, must be carefully defined and reconstructed to assure a correct study design and interpretation of the patterns detected.
Lyle's flying fox Pteropus lylei (Andersen 1908; Fig. 1) is a member of the fruit bat family (Chiroptera: Pteropodidae) and has a distribution that encompasses Thailand, Cambodia, and Vietnam (Bumrungsri et al. 2008). It is a primarily nocturnal species that emerges from roosts at twilight, traveling up to 23 km a night to forage on fruit, nuts, nectar, and flowers in native and cultivated forest plantation, farmland, residential areas, and gardens (Weber et al. 2015, Aziz et al. 2016, Choden et al. 2019. The species is listed as vulnerable on the IUCN Red List due to substantial population declines and in Appendix II of CITES, with trade controlled to avoid overutilization of the species (IUCN 2015, CITES 2015. Threats to P. lylei vary across the species' range. However, empirical evidences suggest that hunting and agricultural intensification that leads to loss of foraging areas and roosting sites are largely responsible for P. lylei population declines through the reduction of foraging availability and roosting sites (Chaiyes et al. 2017).
Abundance is a crucial indicator of wildlife population stability, reflecting factors such as reproductive success, longevity, and a populations' susceptibility to extinction (Morrison et al. 2006). Understanding the factors affecting P. lylei abundance is important to better design conservation strategies for the species in the long term. Pteropus lylei is geographically restricted to South-East Asia (Bumrungsri et al. 2008) and is easily detectable. Indeed, the species is commonly found in urban, rural, and natural landscapes (Weber et al. 2015, Chaiyes et al. 2017, Choden et al. 2019, making it an ideal candidate species to study ecological theories that require high species detectability. The goal of this study was to assess whether the abundance and genetic diversity of P. lylei are associated with environmental features based on the estimated fundamental (N F ) and realized (N R ) ecological niches of the species. Environmental dimensions at a coarse scale are commonly used as a proxy for a species N F -that is, the environmental conditions where the species can survive and maintain populations in the long term without the need for immigration . A species N F is typically reconstructed based on abiotic variables. Here, N F was estimated using climatic data and ecological niche modeling methods (see Peterson et al. 2011). However, P. lylei has a narrow distribution in tropical areas. Thus, abundance data for this species are restricted to areas where climate information may not be particularly heterogeneous (Collen et al. 2008), and this may limit the information provided to estimate N F . In this vein, some studies have failed to support the niche centroid hypothesis based on ecological niche models and species abundances (Dallas et al. 2017, Santini et al. 2019. Other approaches argue that genetic structure and environmental features at the landscape level are the real drivers of species abundance (Gilman 2006, Lester et al. 2007. Therefore, to address this distinction, we also assessed the fine-scale environmental conditions in area of P. lylei occurrence. More specifically, we used landscape-level data as a proxy for the species realized niche N R -that is, the environmental conditions where the species actually occurs, and considers abiotic variables and biotic interactions . Landscape variables were selected based on our understanding of the species roosting habits and known distribution (Duengkae et al. 2019) and included a vegetation index as a proxy for vegetation phenology and foraging resources, Buddhist temples as a proxy for disturbance, or lack thereof, waterbodies as a proxy for water resources, and human density as a proxy for potential competitors or predators. Our findings may help to elucidate the value of the niche centroid hypothesis approach in designing wildlife conservation and management plans.

Occurrence data (N F model)
We estimated a proxy of P. lylei's N F using multivariate climatic space and occurrence data. We obtained P. lylei occurrence data from the following digital public databases: Global Biodiversity Information Facility (GBIF.org 2018), VertNet beta (VertNet 2018), and Mammal Species of the World (Smithsonian National Museum of Natural History 2018). We reviewed the geographic coordinates of each occurrence based on the descriptions of localities in the original databases and only used occurrences from the period 1980-2018. Additionally, we systematically searched for P. lylei's occurrences in scientific journals and theses, using the bibliographic databases Scopus (www.scopus.com) and ISI Web of Science (www.webofknowledge.com), with different combinations of the terms Pteropus lylei and Lyle's flying fox in the search field for the period 2000-2018. We verified the species identity and temporal and spatial precision of all occurrences from these data sources, discarding dubious records (e.g., records outside the known range and based on a questionable identification of the species).
Occurrence records obtained from digital databases and publications were inspected manually using Google Earth to confirm the geolocation corresponded to the metadata of roosting sites. Forest cover was used as a landmark of roosting sites. In total, 36 roosting sites were obtained and used in the ecological niche model calibration (see Model calibration area). The bibliographic details for these records are available in Appendix S1: Table S1.

Climate data (N F model)
We used 15 global bioclimatic variables from the MERRAclim (Vega et al. 2017) for the 2000s period at~4.5-km resolution in GeoTiff format (Appendix S1: Table S2). We excluded variables with artifactual values expressed as abrupt, unrealistic climatic changes between neighboring pixels (i.e., 8, 9, 18, 19, Escobar et al. 2014). Temperature layers (BIO1-BIO11) were determined in degrees Celsius (°C) and humidity layers (BIO12-BIO17) in kg of water/kg of air (Vega et al. 2017). We used principal components analysis (PCA) to generate a subset of independent variables and retained the components that explained ≥90% of the total variation in the climatic dataset (Wang et al. 2016).

Model calibration area
At the coarse scale, study sites were selected in countries that overlap the known distribution of P. lylei according to the most recent assessment for Thailand (Chaiyes et al. 2017), Cambodia, and Vietnam (Bumrungsri et al. 2008). The total study area covered 159,200 km 2 (8-16°N, 98-108°E; Fig. 2).
The BAM framework is an interplay between geographic accessibility, abiotic factors, and biotic environmental conditions (Sober on and Peterson 2005, Barve et al. 2011). Models, when calibrated, should encompass the areas where species can disperse-that is, areas that can be colonized by a species or M (Movement) in the BAM framework (Sober on and Peterson 2005). In this framework, A is the region in the geographic areas where abiotic conditions from N F are available (e.g., suitable temperature). B represents the geographic areas where biotic conditions allow the presence of the species (e.g., availability of fruits). M denotes the region that has been accessible to dispersal or colonization by the species over some relevant time interval (e.g., areas limited by the ocean or mountains). We calibrated an ecological niche model based on the species' dispersal potential, M. We determined M as the average distance among available occurrences for P. lylei and generated a buffer around occurrences using this distance (Poo-Muñoz et al. 2014, Cooper andSober on 2018). We assumed that this area formed the boundary of the areas where P. lylei have historically dispersed. We used this calibration area to limit the environmental layers and unique species' occurrences to generate an ecological niche model representing P. lylei s' N F . The N F model was developed using NicheA version 3 (see Qiao et al. 2016). NicheA assumes that species will show a Gaussian response to environmental variables so that a multivariate model will result in a minimum volume ellipsoid ). The niche centroid was estimated based on Euclidian distances in NicheA considering the shape, position, and size of the ellipsoid ).

Model evaluation
We split the occurrence coordinates in four quadrants with similar numbers of occurrences. Two off-diagonal quadrants were used for calibration and one for evaluation (Appendix S1: Fig. S7). Model evaluation was conducted shifting the evaluation and calibration quadrants for a two-way test. Evaluations included the total number of trials (total number of evaluation occurrences) and total number of successes (i.e., evaluation occurrences predicted correctly), and the proportion of area predicted by the model built with the calibration points. These values were used to develop a cumulative binomial distribution test with a significance of a = 0.05 as described elsewhere ). The final model used all the occurrences available to us.

Landscape data (N R model)
At the landscape level, our study area covered 104,100 km 2 (10-15°N, 98-103°E; Fig. 2) and included 25 provinces of Western, Central, and Eastern regions of Thailand. Landscape-level environmental variables considered as proxies for N R were assessed according to their potential influence on the abundance of P. lylei. N R variables included the normalized difference vegetation index (NDVI), annual mean temperature, percent cover of waterbodies, number of Buddhist temples, and human density to quantify their associations with P. lylei abundance. All were measured at a spatial resolution of 30 9 30 m. The NDVI reflects temporal and ❖ www.esajournals.org spatial trends and variation in vegetation productivity associated with animal resources (Pettorelli et al. 2005, Trierweiler et al. 2013. Food availability at the landscape level could affect the number of flying foxes at a roosting site (Giles et al. 2016). We used NDVI data derived from Landsat 8 satellite images at a 30-m spatial resolution, which comprised 13 tiles of the original images for the period January 2015-December 2016 (Landsat.gsfc.nasa.gov 2018; Fig. 2A).
Landscape-level temperature can affect individuals' survival and roost selection (Bronrier et al. 1999, Kerth et al. 2001). The optimal ambient temperature range, without regulatory changes in metabolic heat production of animal, lies between 24°and 35°C (Neuweiler 2000, Kingma et al. 2012). Heat stress affects flying foxes when temperatures reach ≥42°C causing mortality (NSW Government 2018). Annual mean temperature (°C) between 1981 and 2010 was obtained from meteorological stations of the Thailand Meteorological Department (TMD) and inverse distance weighting to raster format at a spatial resolution of 30 9 30 m (Fig. 2B).
Waterbodies such as rivers, streams, and lakes are crucial as a source of drinking water for bats (Rainho and Palmeirim 2011). It has been reported that P. lylei roost sites are usually found close to water (Thanapongtharm et al. 2015, Chaiyes et al. 2017. Waterbody data were obtained from the Ministry of Natural Resources and Environment (2018), Thailand, and percent of waterbodies per 1 km 2 was estimated in raster format (Fig. 2C).
Additionally, we accounted for variables representing anthropogenic disturbance (Mart ınez-Guti errez et al. 2018): number of Buddhist temples and human density. Buddhist temples are generally undisturbed animal sanctuaries according to religious practices and Thai law. Therefore, temples may provide shelter, protection, and more suitable roost trees for bats that improve the conditions for their occurrence (Boonkird et al. 2006). Locations of temples were obtained from the National Office of Buddhism (2018) (Fig. 2D). Finally, human demography data by subdistrict from the 2015 national census were obtained from the National Statistical Office (2018). A human density map was then computed by dividing the population of each subdistrict by its area. The final map represents an estimate of population density per square kilometer (Fig. 2E). We used a 23-km buffer (maximum linear distance between day roosts and foraging areas; Weber et al. 2015) around roost sites to estimate the average values of each landscape-level variable, as a proxy of the population range, which was used in the regression between abundance and genetic diversity with the landscape variables (see Abundance, genetic diversity, and environmental factors).

Species abundance
We determined the abundance of P. lylei at 30 roosting sites across Thailand based on direct counting of individuals from November 2015 to October 2016. Details regarding location of P. lylei roosts and the estimation of abundances are described in Chaiyes et al. (2017). Briefly, we geolocated each roost site and confirmed the presence of the species by morphological criteria (Francis 2008). We performed individual counts at roosting sites during daylight hours using the direct sighting technique, digital cameras, and drones, according to roost accessibility and population size. When digital cameras and drones were used, we took photographs of and around each colony and counted the number of individuals from the photographs (Fig. 2). This dataset represents one of the most comprehensive estimations of P. lylei abundances across its geographic range.

Genetic diversity
We derived genetic diversity based on haplotype diversity (Hd), average number of nucleotide differences (k), and nucleotide diversity (p) for P. lylei from Sukgosa et al. (2018) (Appendix S1: Table S1).

Abundance, genetic diversity, and environmental factors
We used a Spearman's coefficient to evaluate the correlation between P. lylei abundance and genetic diversity and the distance to the niche centroid (coarse-scale or N F model) and landscape-level variables (fine-scale or N R model). We quantified value of the relationship between each landscape-level variable and distance to the niche centroid with abundance and genetic diversity of P. lylei. Multiple regression analysis was used to quantify the linear association between abundance and genetic diversity and landscape-level variables. The stepwise model selection procedure was used to identify the best model based on the Akaike information criterion (AIC; Akaike 1973). Significance was determined at P < 0.05. All statistical analyses were performed using R version 3.2.2 (R Core Team 2016).

RESULTS
In total, we obtained 30 spatially unique records of P. lylei presence for Thailand, three for Cambodia, and three for Vietnam (Fig. 3, top). The total area considered for model calibration was approximately 4.2 9 10 5 km 2 . The first three principal components accumulated 99.6% of the original climatic information and were used in the estimation of P. lylei's N F (Fig. 3). Model evaluations indicated significant predictive power of models to anticipate independent data (P < 0.001; Appendix S1: Fig. S6). The ecological niche model revealed that areas in close to P. lylei niche centroid are mostly in the Thailand, Cambodia, and Vietnam and included area of southern of Lao.
Abundance and genetic diversity of P. lylei and distance to the centroid Correlations between abundance and genetic diversity of P. lylei and distance to the centroid were not statistically significant. Thus, we found no correlation between P. lylei abundance and the distance to the centroid in the N F model (Appendix S1: Fig. S1A; q = 0.087). Similarly, we found no correlation between any genetic diversity metrics and distance to the niche centroid in the N F model (Appendix S1: Figs. S2A, S3A, S4A; Hd q = À0.081; k q = 0.358; p q = 0.293).
Correlations between average number of nucleotide differences (k) and nucleotide diversity (p) of P. lylei and number of temples were significantly (P = 0.005), whereas the other models revealed nonsignificant results (Table 1).

DISCUSSION
We investigated spatial patterns of P. lylei demographics in relation to the center of their ecological niches. Specifically, we assessed whether the distance to the niche centroid or landscape information was correlated with geographic patterns of abundance and genetic diversity of P. lylei. Our results indicate no correlation between distance to the niche centroid from a coarse-resolution presence-only distribution model resembling N F and spatial variation in the abundance and genetic diversity of P. lylei. The landscape-level variables resembling N R (i.e., NDVI, annual mean temperature, waterbody area, density of humans, and number of temples) are thought to be important in driving the species distribution (Neuweiler 2000, Boonkird et al. 2006, Kingma et al. 2012, Trierweiler et al. 2013, Thanapongtharm et al. 2015, Chaiyes et al. 2017). However, these fine-resolution variables were not linked to P. lylei abundance and genetic diversity, with the exception of annual mean temperature and number of temples, respectively.
Despite these findings, there are obvious limitations to our methods. The climatic data used by our niche models were averages over space and time rather than the specific conditions that occur at the time and place of distribution. Despite this shortcoming, the goal of our study was to identify where general climatic (B) Ecological niche model based on P. lylei occurrences in the form of a minimum volume ellipsoid (blue ellipsoid). Environmental space (gray points) constructed using the first three principal components (X, Y, and Z). characteristics permit P. lylei occurrence, which we accomplished successfully. Additionally, logistic limitations did not allow us to include abundance and genetic diversity data for P. lylei from population in Cambodia and Vietnam. Therefore, our assessment is restricted to information from Thailand. However, the geographic coverage of our sampling across Thailand was considerable and representative of the species range (Chaiyes et al. 2017). In this study, we accounted for rivers and lakes, but streams or human-made waterbodies were neglected in our analysis, even though they are drinking water sources for bats (Rainho and Palmeirim 2011). Our simplistic model, based on ordinal-scale assumptions, attempted to capture correlative signals from the data available. However, other, more complex, regression methods may be more sensitive to capture signals of association (Mart ınez-Guti errez et al. 2018). Martin et al. (2016) reported a negative relationship between distance to the niche centroid and abundance of P. alecto and P. conspicillatus, but no relationship between distance to the niche centroid and abundance of P. poliocephalus and P. scapulatus. Other studies provide an empirical test of the niche centroid hypothesis and have found inconsistent relationships between abundance and distance to the niche centroid (Manthey et al. 2015, Dallas et al. 2017.
Our results showed annual mean temperature to be negatively correlated with P. lylei abundance. According to Dey et al. (2013), P. giganteus (Indian flying fox) abundance and ambient temperature are negatively correlated and high temperature is linked to mass die-offs. Moreover, temperature affects food sources of flying foxes, by influencing the flowering of various tropical fruits (i.e., mango and durian) influenced by the temperature (Hariyono 2017, Rajatiya et al. 2018). While we did not find an association between abundance and NDVI, others have found frugivorous bat abundance to be correlated with vegetation structure (e.g., canopy openness and tree height; Vleut et al. 2015). This suggests a plausible response of flying fox abundance to other landscape variables not considered here (Ali 2010, Dey et al. 2013 or an artifact result due to the limited period analyzed for the NDVI values. Interestingly, we found no association between P. lylei abundance and genetic diversity, suggesting that optimal conditions may not correspond to a greater population stability.
Although our analyses do not support the previous claims that species are more abundant or genetically diverse at their niche centroid (Hengeveld and Haeck 1982, Brown 1984, Holt 2009, Mart ınez-Meyer et al. 2013, Lira-Noriega and Manthey 2014, Yañez-Arenas et al. 2014, Manthey et al. 2015, Ram ırez-Gil et al. 2018, the fact that only annual mean temperature at landscape level can explain changes in the abundance of P. lylei shows the complexity and challenges behind the drivers of species abundances. For example, even when the species is highly detectable and a well-structured survey of its distribution, genetic composition, and abundance was conducted across most of the species' range (Sagarin and Gaines 2002, Tuya et al. 2008, Dallas et al. 2017, Santini et al. 2019, the environmental data did not explain the configuration of the species demography, suggesting that abundances may be unpredictable at least 0.641 Nucleotide diversity = 6.79e À01 + 2.46e À02 X 1 À3.13e À05 X 2 + 1.31e À04 X 3 À2.47e À02 X 4 + 1.36e À03 X 5 0.064 À95.713 0.870 Nucleotide diversity = 1.11e À02 À1.91e À05 X 2 0.005 À93.694 0.647 Notes: X 1 : normalized difference vegetation index (NDVI); X 2 : number of Buddhist temples; X: human density; X 4 : annual mean temperature; X 5 : percent cover of waterbodies; AIC: Akaike information criterion; R 2 : sum of squared residuals. Values in italics are significant at the P < 0.05 from the best-fit models.
under the environmental dimension we explored at coarse and fine resolutions. Instead, other biotic factors may be stronger drivers of population size of P. lylei.
We divided our analysis in coarse and fine scales by using different environmental datasets. Previous studies have mixed coarse-and finescale variables for the assessment of the niche centroid hypothesis (Santini et al. 2019). However, the classic ecological niche modeling theory cautions against mixing coarse-and fine-scale variables to facilitate model interpretation (Sob eron and Peterson 2005). We anticipate that models calibrated with climatic and landscapelevel environmental data could generate different results at the cost of challenging model interpretation (Levin 1992.
Pteropus lylei is threatened across its range due to illegal hunting and loss of roosting habitat, as existing trees are removed but not replaced (Bumrungsri et al. 2008, Chaiyes et al. 2017. In Thailand, illegal hunting areas occupied more than 50% of P. lylei foraging zones (Chaiyes et al. 2017). Pteropus lylei is able to persist in a landscape profoundly disturbed by human land use, which suggests that roost sites and food resources may not be dependent on dense forest (Weber et al. 2015). For example, P. lylei colonies commonly occur on the grounds of Buddhist temples where hunting is not allowed (Boonkird et al. 2006, Thanapongtharm et al. 2015, Chaiyes et al. 2017.
A few other studies of genetic variation in flying foxes reveal that flying foxes are essentially panmictic species. The genetic homogeneity of flying foxes over their extensive ranges concurs with evidence from direct observation of high dispersal potential of individuals (Webb and Tidemann 1996). Pteropus lylei individuals have shown a maximum distance traveled per night that ranges between 2.2 and 23.6 km in Thailand and from 6.88 to 105 km in Cambodia (Weber et al. 2015, Choden et al. 2019). However, we found significant correlations between average number of nucleotide differences (k) and nucleotide diversity (p) of P. lylei and number of temples. This finding suggests that temples play an important role in the P. lylei metapopulation dynamics and gene flow at the landscape level. Indeed, human-wildlife conflict could help explain the heterogeneous distribution of P. lylei populations and genetic structure. Flying fox's (P. giganteus) foraging activities can damage a wide variety of fruit crops, often causing considerable economic losses, which exacerbates conflicts (Aziz et al. 2017). The noise and smell from bat colonies, as well as concerns about disease transmission (e.g., Nipah virus), often result in conflicts with humans (Currey et al. 2018). These human-wildlife conflicts may shape the distribution of flying fox populations in urban and rural areas (Aziz et al. 2017, Currey et al. 2018.
In summary, this study highlights our limited ability to reconstruct the distribution of abundance and genetic diversity of P. lylei in South-East Asia. We were unable to identify relationships between patterns of abundance using coarse-scale environmental variables but were able to identify relationships at a finer scale. Further research is required to understand the effects of human presence on the distributional change of P. lylei populations. Identifying the rules explaining the abundance of species across their range is still a major challenge in modern ecology. Threats PREDICT program. We thank Paige Van de Vuurst for her comments and editions that improved the final version of this manuscript. Finally, we would like to take this opportunity to thank you for reviewers for taking the time and effort necessary to provide insightful guidance.

LITERATURE CITED
Akaike, H. 1973. Information theory and an extension of the maximum likelihood principle. Second International Symposium on Information Theory in 2nd International Symposium on Information Theory