Environmental drivers of spatiotemporal foraging intensity in fruit bats and implications for Hendra virus ecology

In the Australian subtropics, flying-foxes (family Pteropididae) play a fundamental ecological role as forest pollinators. Flying-foxes are also reservoirs of the fatal zoonosis, Hendra virus. Understanding flying fox foraging ecology, particularly in agricultural areas during winter, is critical to determine their role in transmitting Hendra virus to horses and humans. We developed a spatiotemporal model of flying-fox foraging intensity based on foraging patterns of 37 grey-headed flying-foxes (Pteropus poliocephalus) using GPS tracking devices and boosted regression trees. We validated the model with independent population counts and summarized temporal patterns in terms of spatial resource concentration. We found that spatial resource concentration was highest in late-summer and lowest in winter, with lowest values in winter 2011, the same year an unprecedented cluster of spillover events occurred in Queensland and New South Wales. Spatial resource concentration was positively correlated with El Niño Southern Oscillation at 3–8 month time lags. Based on shared foraging traits with the primary reservoir of Hendra virus (Pteropus alecto), we used our results to develop hypotheses on how regional climatic history, eucalypt phenology, and foraging behaviour may contribute to the predominance of winter spillovers, and how these phenomena connote foraging habitat conservation as a public health intervention.


Results
Data filtering and spatiotemporal aggregation of foraging areas used by our 37 GHFFs rendered 385 unique counts of foraging stops per 700 m grid cell within the study area from June 2013-Feb 2015 (Fig. 1). These count data represent counts of foraging stops per unit area and time, which we refer to as a measure of 'foraging activity' or 'foraging intensity' . Model selection criteria based on performance within the 5-fold spatial cross validation structure returned 27 models giving the mean Poisson deviance and Mean Squared Error (MSE) ( Table 1). The selected model (tree complexity (tc) = 9, learning rate (lr) = 0.001, bagging fraction (bf) = 0.6, number trees = 3400) had the lowest mean Poisson deviance (2.54) and MSE (4.11) while fitting >1000 trees. When the selected model was refitted with all data, it fitted 5100 trees. The final boosted regression tree model provided a valid fit to the data with evenly distributed residuals and an R 2 value of 0.48 (Fig. 2). Distance to roost was the most influential covariate in the model with a Relative Influence (RI) of 13.99% (Table 2). Other covariates with high RI include short-term changes in temperature and eucalypt vegetation indices such as ECARR_3mo (5.85%), AVGtmn_1mo (5.12%), ANOMtmn_1mo (4.23%), and EWDI1 (4.2%) ( Tables 3 and 2). Partial dependence plots (Fig. 3) show that the fitted function is 1-2 standard deviations higher for foraging areas within 25 km of the home roost. Predicted foraging intensity in response to all vegetation indices (EWDI1, ECBRR, and ECARR; green plots in Fig. 3) was one standard deviation below the mean fitted response when the difference in the vegetation index over 9 or 12 month time-lags was negative, but the fitted response was one standard deviation above the mean when the indices were positive, indicating an increase in photosynthetic activity or tree canopy density over these longer time-lags. Over shorter time intervals (3-month time-lags and the current month), the fitted response was up to 2 standard deviations above the mean when the difference in vegetation indices was negative, indicating that decreases in vegetation indices over the preceding 1-3 months are associated with a higher level of foraging intensity.
In terms of climatic conditions, temperature variables (red plots in Fig. 3) show that predicted foraging intensity falls dramatically when average minimum temperature (AVGtmn_1mo) rises above 20 °C, but it increases when the minimum temperature anomaly (ANOMtmn_1mo) is negative (representing an especially cool month). Precipitation variables (blue plots in Fig. 3) show that severe decreases in precipitation over 6 months (PREC_1mo_6mo) are not associated with foraging activity, and high values (>1000 mm) of cumulative rainfall over the preceding 9 months (PREC_9mo) coincide with a dramatic decrease in predicted foraging activity. Interactions between environmental variables provide further description of the fitted function. Specifically, the highest levels of foraging activity are fitted when both of the following conditions are met: 1) monthly precipitation is significantly less than the prior year and <200 mm of rainfall in the preceding 3 months, and 2) the minimum temperature anomaly of the preceding 9 months is low and coincides with a decrease in ECARR over 3 months (Fig. 4). In general, the partial dependence plots and interactions indicate that foraging activity is more intense when summer and autumn conditions are cooler than average with relatively reduced precipitation compared with previous months, and photosynthetic productivity in eucalypts diminishes. This emphasizes the importance of past conditions (3-12 months previous) that are warmer and wetter, suggesting that a moderate transition from warm/wet to cool/dry is a primary environmental driver of eucalypt flowering and flying-fox foraging intensity within our study area. When the final model is projected across the study area, foraging intensity exhibits a patchy spatial distribution that is qualitatively consistent with typical heterogeneous patterns of eucalypt flowering (Fig. 5). Similarly, when plotted as a time series and aggregated by year (Figs 6 and 7 respectively), model predictions also display interannual variability. Model predictions summarized by month for foraging areas within 30 km of the Canungra roost show that the measure of spatial resource concentration (QVC) is typically higher in mid-summer and early autumn (Dec-Mar; Fig. 7). Beyond seasonal patterns, foraging areas around Canungra exhibit some interannual variation also, with mid-summer months in 2010 and 2015 being particularly high (9). We found that independent monthly population counts of Pteropus poliocephalus observed at Canungra from 2008-2014 were significantly correlated with predicted values of monthly foraging intensity (r = 0.74, 95% CI = 0.28, 0.92, p = 0.006, Figure 1. Study area of analysis, defined by the maximum observed foraging radius (46 km for Boonah and 30 km for Canungra), is plotted with unfiltered foraging areas (red points). Location of the roosts at Boonah and Canungra are shown as a green triangle and square respectively. The four black exes indicate foraging areas that were deemed outliers based on their proximity to different roosting sites. The black line represents the maximum foraging radius for Pteropus poliocephalus and the black triangles show locations of additional roosting sites. The background shows the distribution of the Eucalypt Chlorophyll a Reflectance Ratio (ECARR) which indicates areas with dense vegetation as light yellow and areas of less-dense vegetation as light blue. This map was generated using ESRI ArcGIS Desktop: Release 10 (www.esri.com).
SCIENTIFIC REpoRts | (2018) 8:9555 | DOI:10.1038/s41598-018-27859-3 n = 12; Fig. 7), giving strong validation that the fitted model adequately captures seasonal trends in foraging intensity. Comparison of the model to population counts across years did not show a significant correlation (r = 0.31, 95% CI = −0.58, 0.86, p = 0.5, n = 5). However, from 2008-2014, the model predicts the highest foraging intensity in 2010 and 2013, the same years in which median population counts are highest. More broadly, time-lagged cross correlation reveals that model predictions are significant and positively correlated with the SOI values 3-8 months previous (n = 120; Fig. 10), as seen by slightly lagged alignment of peaks in predicted foraging intensity with peaks in SOI values (Fig. 11). To explore our model of spatial resource concentration in the context of Hendra virus ecology, we plotted our model results along with existing data of Hendra virus prevalence at the Boonah roosting site 26 and spillover events in horses 27 that allowed comparison of seasonal trends among the three phenomena (Fig. 8). The Generalized Additive Model (GAM) that we fit to QVC values (n = 120, P-value of smooth term (month) = 3.7 × 10 −9 , deviance explained = 31.4%) shows that the QVC is highest in late-summer and lowest in winter (Fig. 8a). The GAM fitted to Hendra virus prevalence data (n = 40, P-value of smooth term (month) = 1.2 × 10 −4 , deviance explained = 45.4%) showed a peak in prevalence during winter (Jun-Aug; Fig. 8b). When we plotted the total counts of Hendra virus spillover incidents, all three distance thresholds (25,50, and 100 km) showed the same seasonal pattern with the highest count occurring in June with decreasing counts in subsequent months. Within 25 km of our study roost sites there have been a total of: 3 spillover incidents which have occurred within 25 km during Jun-Jul, 7 spillover incidents within 50 km during Jun-Oct, and 20 spillover incidents within 100 km during Jun-Nov (Fig. 8c).

Discussion
We developed a spatiotemporal model of flying-fox foraging intensity that is validated by independent data of seasonal bat population counts and provides a quantitative link between bat foraging dynamics and cycles of the ENSO. Specifically, increased rainfall and minimum temperatures (typical of La Niña events and favorable for the growth of eucalypts) in preceding seasons drive flowering and spatially concentrated foraging intensity in the following autumn. We also found that in the same region as our study area, both observed prevalence of Hendra virus in bat populations, and occurrence of spillover to horses, align with seasonal and annual periods when our model indicates a low concentration of nectar-based resources (Fig. 8). Nectar scarcity may drive bats to feed in proximity to horses or may modulate bat immune status leading to virus shedding 13,14 . Further, future climate changes will likely alter climatic conditions during important growth phases of eucalypts 28 , affecting the spatiotemporal distribution of nectar-based foraging resources for flying-foxes, which could have implications for both Hendra virus ecological dynamics and strategies for habitat conservation.
Previous research into the movement ecology of bats has shown important aspects of mobility and habitat preference. Telemetry and GPS tracking studies have looked at the scale of movement [29][30][31][32][33] , habitat preference 34,35 , and seasonal differences in foraging behavior between wet and dry seasons 36 in several species of flying-fox, which emphasize that mobility and generalist foraging behavior allow population redistribution in response to food resources. Previous studies have also modeled the spatial extent of bat habitat using roost locations [37][38][39] . However, models based on roost locations are limited when describing foraging habitat. Giles et al. 40 modeled population occupancy and abundance at day roosts with proxies of eucalypt flowering. They addressed the discrepancy between roosting and foraging habitat by taking a spatially weighted sample of model predictors within

ECARR_3mo
Change in ECARR over 3 months 5.85

AVGtmn_1mo
Average minimum temperature of preceding month 5.12

ANOMtmn_1mo
Cumulative minimum temperature anomaly over preceding month 4.23

ECBRR_6mo
Change in ECBRR over 6 months 3.70

PREC_1mo_6mo
Change in PREC_1mo over 6 months 3.29

EWDI1_3mo
Change in EWDI1 over 3 months 3.08

EWDI1_12mo
Change in EWDI1 over 12 months 3.08

PREC_9mo
Cumulative  areas known to contain diet species of bats. However, robust models of foraging habitat require explicit data of occupied foraging areas. The high-resolution tracking data collected in Westcott et al. 4 allowed us to identify these foraging areas with high confidence. When combined with spatial proxy variables of eucalypt phenology from Giles et al. 40 we built a model of foraging intensity at the regional scale, which serves as a starting point for developing landscape scale models. Environmental variables that predict bat foraging intensity are consistent with known climatic triggers of eucalypt flowering. We found that, in our study area, high predicted values of foraging intensity were associated with long-term increases (9-12 months) but short-term decreases (1-3 months) in eucalypt vegetation indices (green plots in Fig. 3). This suggests that there is a sequence of key seasonal time periods prior to foraging when phenological changes in eucalypts occur (e.g. growth phase, bud production, anthesis), leading to flowering and nectar availability in the following months. We propose then, that a slowing of vegetative growth in the eucalypt canopy is an important precursor to foraging activity, which is consistent with the reproductive phenology of eucalypts. Generally, growth and flowering are separate processes, where vegetative growth is suppressed in favor of bud development and floral initiation 41,42 . In many species, these changes from one phenological stage to the ANOMtmn_1mo Cumulative minimum temperature anomaly over preceding month      43,44 ].
In terms of climate, we found that foraging intensity was higher during periods that are cooler and drier than average, especially if they are preceded by conditions that are warmer and wetter than average. Interaction plots (Fig. 4) show that flowering and foraging is favored when precipitation slows over the preceding year, indicating that previous rainfall followed by a transition from wet to dry is important. Similarly, Law et al. 8 observed that, for 20 species of eucalypts in northern New South Wales, the greatest flowering in a 10 year period occurred 9 months after an above average rainfall event. Further, it is interesting that partial dependence plots (blue plots in Fig. 3) indicate that in the region of our study area, changes in precipitation must be relatively mild, and that extreme transitions (e.g. an exceedingly wet period followed by an average one, or an average period followed by a drought) may not be amenable to flowering. We speculate that during wet periods, eucalypt forests may forgo  When considering temperature, we found that flowering and foraging is associated with low minimum temperatures (below 19 °C) in the preceding month (red plots in Fig. 3). Likewise, Law et al. 8 showed that prior cool temperatures were the most consistent climatic trigger for flowering in 9 of 20 eucalypt species in northern New South Wales. In more southern species, Specht and Brouwer 44 note that in E. regnans the growth phase is halted when mean temperatures drop below 16-18 °C, followed by flowering [typically in autumn for this species 43  Beyond seasonal climatic drivers, we observed a positive correlation between the spatial concentration of nectar foraging and values of the SOI at a temporal lag of 3-8 months (Fig. 10), indicating that previous climatic conditions associated with La Niña have a cumulative effect on flowering phenology and resulting foraging by flying-foxes. It has been known for some time that the state of the ENSO can predict spring/summer rainfall for eastern Australia 20 , so the observed correlation is not entirely surprising. SOI values above zero are associated with La Niña cycles, where increased ocean surface temperatures increase cloud cover, humidity, rainfall, and minimum temperatures for southeastern Queensland 46 . Some species of eucalypt may invest more energy into the growth phase during the amenable conditions of La Niña, until drier and cooler conditions of a neutral phase trigger flowering. However, it should be noted that each species of eucalypt has a unique seasonal schedule for growth and flowering phenology that is influenced by both photoperiod and time-lagged climatic effects 47 . Therefore, physiological responses to the broader climatic conditions brought by La Niña may influence eucalypts within a region by means of the Moran effect 48 , causing periodic deviations from strict seasonal dynamics. We hypothesize, accordingly, that the timing of ENSO cycle changes (i.e. La Niña to neutral phase) is a general driver of the notorious inter-annual periodicity of eucalypt growth and flowering phenology, and the erratic patterns of redistribution characteristic of flying-fox populations.
Seasonal and annual trends in model predictions indicate that nectar-based foraging for flying-foxes has important implications for Hendra virus disease ecology. As shown in Fig. 8a,b, seasonal trends in the model predict peaks in spatial resource concentration in late-summer, followed by lower values in autumn/winter. Coincidentally, Field et al. 26 and Paez et al. 49 show that peak Hendra virus prevalence in bat populations also occurs in the winter, which may indicate epidemiological drivers of Hendra virus prevalence in bat roosts in response to resource scarcity. It should be noted that our model of resource distribution is based on the foraging patterns of GHFF whereas BFF is thought to be the primary driver of Hendra virus transmission in bat populations and excretion of the virus into the environment 15 . However, Paez et al. 49 found that measures of climate in preceding seasons (similar to those used in our model) in combination with the number of BFF could predict peaks in Hendra virus prevalence. Specifically, they found that when hot and dry conditions in   the previous spring/summer coincide with larger numbers of BFF, peaks in Hendra virus are higher and occur earlier in the winter season. Based on our model and its correlation with the SOI, warmer and drier conditions in the spring/summer (typical of the El Niño phase of ENSO), appear connected with both reduced flowering and foraging in the following autumn/winter, and larger peaks in Hendra virus prevalence in the winter season. Following Plowright et al. 14 , a period of resource scarcity could impact viral prevalence in bat populations through a decreased allocation of energy to immune defenses. Resource scarcity could also increase contact rates and connectivity between local roosts due to more dispersed foraging behavior in the region, thus affecting viral dynamics.
Spillover events from flying-foxes to horses in the region follow a similar pattern with winter occurrence (Fig. 8c) 13 . Plowright et al. 13 also note that there was an unprecedented cluster of spillover events in the winter of 2011 in southern Queensland and northern New South Wales, a year in which our model predicts the lowest concentration of nectar-based resources around Canungra (Figs 7 and 9). We, therefore, hypothesize that periods of a low spatial concentration of nectar resources result in dispersed foraging behavior among resident GHFF and BFF individuals, which increases the number of bats foraging in horse paddocks. In this manner, the type of population distribution across the landscape (e.g. aggregated during resource plenty or dispersed during resource scarcity) determines where the virus can be excreted into the environment, providing the first enabling factor for spillover to horses 50 .
The role of foraging behavior in bat population distribution and Hendra virus ecology emphasizes the importance of coherent management actions which preserve or enhance flying-fox foraging habitat. This is especially timely because a 2013 amendment to the Queensland Vegetation Management Act 1999 loosened regulations on remnant and old growth forest clearing. A subsequent resurgence of vegetation clearing [395,000 hectares per year in 2015-16 51 ], put Queensland on a trajectory to return to the clearing rates of the 1990s when it was labeled a world deforestation 'hotspot' 52,53 . In response, the Queensland government reinstituted legislation this year to  protect remnant vegetation 54 , which is a promising step to reduce future clearing rates. However, the impacts of clearing on bat populations and risk of spillover may not be observed until subsequent foraging bottlenecks occur due to annual or inter-annual cycles in climate. For instance, cleared foraging areas for bats may not cause significant changes in nectar resource availability and foraging behavior until the next El Niño phase when warmer and dryer spring/summer conditions cause poor flowering in the following autumn. The time-lagged effects of deforestation and habitat fragmentation, therefore, present significant challenges when studying the ecological impacts of landscape change. Hence, habitat conservation strategies should employ analyses on the dynamic changes in flowering and foraging over time to identify areas that: 1) offer consistent winter flowering species, and 2) contain species that produce large amounts of nectar on an inter-annual basis, and seek to preserve these areas, or even establish new areas to plant important diet species. Such a conservation strategy targets both temporal and spatial gaps in nectar-based resources 24 .
Our model is subject to three important limitations inherent in the data that must be considered. First, the data that we used to define foraging areas of bats came from 37 individuals roosting at Boonah and Canungra between June 2013-February 2015. Although we have confidence in our definition of foraging areas, we must assume that the areas used by our 37 individuals are representative of a larger population. Second, there is likely spatial autocorrelation in the foraging area data (even after we clustered foraging stops into areas). We, therefore, designed a modeling approach that reduces this influence by aggregating foraging areas to more coarse spatial and temporal scales, utilizing the 'bagging' technique within the BRT algorithm, including distance from the home roost as a predictor, and by determining model meta-parameters using a 5-fold spatial cross validation structure. Even with appropriate precautions, it is difficult to completely remove the confounding effects of autocorrelation in the data. Third, there is heterogeneous species richness in eucalypts within the study area, which may lead to diverse responses to climate. There is likely variability in foraging dynamics that is not included in our model because the model predictions are summarized over coarse spatial and temporal scales in order to draw broader inferences. This is a trade-off that must be made in order to describe the general environmental drivers of flowering and bat foraging within the region of our two study roosts over time. Therefore, our conclusions are relative to observed foraging activity within our study area, and differences in eucalypt species richness and regional climate regimes may lead to different results in other regions.
So then, are we any further along than Ratcliffe in implementing 'better-informed conservation'? Recent dispersal actions 55 indicate that, in terms of implementation, the answer is no. However, we are making progress in understanding the dynamics of bat population distribution, which is difficult given the spatial scale at which they operate. Here we showed that the climatic responses of eucalypt phenology drive bat foraging, which has important implications for bat population distribution, Hendra virus ecology, and conservation strategy. This provides additional ecological context to previous work, and connects ecosystem phenology to broader climatic drivers such as the ENSO, allowing speculation about the future impacts of climate changes. For example, under future climate predictions, the frequency of El Niño cycles are expected to increase, with more irregular timing 56 . Based on our data and model, this would result in less reliable temperatures and rainfall in spring/summer seasons, contributing to a decrease in the spatial concentration of flowering and foraging in our study area. Therefore, future research that takes the concepts presented here, and applies them to the landscape scale could produce a useful spatiotemporal map defining priority conservation areas that provide sufficient foraging habitat over time. Additionally, mechanistic models of bat distribution that include this information would help to understand the spatiotemporal distribution of bat populations in a human-dominated landscape, and the potential for food resource phenology to drive Hendra virus ecological dynamics.

Methods
Environmental data. We developed a suite of environmental data layers for the rectangular extent encompassing the maximum foraging radii of two study roost sites (Boonah and Canungra; Fig. 1) from 2006-2015 at 8-day intervals (spatial resolution of 700 m in Albers Equal Area projection). These include Eucalypt-based vegetation indices, average monthly precipitation, cumulative precipitation, average monthly temperature extrema, and temperature anomaly. The Eucalypt-based vegetation indices are spatial adaptations of laboratory measurements from Datt 57,58 , which represent chlorophyll and water content of Eucalypt leaves. For a more extensive discussion of the indices and climate variables, see Giles et al. 40 . Given the significant influence of prior environmental conditions on the timing and intensity of Eucalypt flowering 8,28 , we also derived time-lagged differences of these variables at 1, 3, 6, 9, and 12 month time lags. Also, we added distance from the home roost as a predictor variable to help model spatial autocorrelation of the foraging data not addressed by data aggregation and modeling techniques such as bagging and spatial cross validation (discussed below). For a complete list of variables explored in the model, see Table 3. We then aggregated the raw data on two levels: first, by clustering movement data into foraging stops, and second, by counting the number of points aggregated to the spatial and temporal resolution of the environmental data. Fine-scale movement data was filtered to identify times and locations of significant foraging activity. We considered bats to be in a 'stopped' state when the speed of the device was recorded as <2 ms −1 . We defined a successful foraging event (hereafter referred as a foraging 'stop') as an area where a group of stop points was within 10 meters of one another. Accordingly, we produced a dataset of foraging 'stops' by clustering points within 10 m and taking the centroid, rendering 2997 points with a spatial distribution shown in Fig. 1. Although the data have an even distribution across gender (female = 1440, male = 1557), season (summer = 1318, winter = 1679), and the two home roosts in our study area (Boonah = 1590, Canungra = 1407), movement data was not retrieved from any bats in Mar/Apr and Sep/ Oct, and the bats exhibit a preference for foraging near Canungra in the summer and Boonah in the winter (Fig. 12). For more detailed specifications of the GPS technology developed for the tracking collars, see Sommer et al. 59 , and for additional description of their implementation to collect movement data, see Westcott et al. 4 . To remove the potential confounding effect of generalist foraging behavior compared with foraging on nectar-based resources in Eucalypt forests, we removed points within 500 m of urban and peri-urban land use types as this was expected to contain high proportions of non-native fruiting trees, leaving 2416 foraging stops that are likely driven by Eucalypt phenology (Fig. 13). To attain counts of foraging stops per unit area and time, we aggregated the data to the temporal and spatial resolution of the environmental data. We aggregated foraging stops temporally by assigning them to the nearest 8-day interval of the environmental data, and then spatially by counting the number of foraging stops within each grid cell, resulting in aggregated counts for 385 unique locations and times.
To account for sampling bias and provide information on the environments of foraging sites that are available, but not used, by our 37 individual bats, we implemented a use-availability sampling design 60,61 . Specifically, we compiled a data set of 'available' foraging sites, for which the response variable (number of foraging stops) was set to zero, by drawing a large number of randomly stratified points within the maximum observed foraging distance of each roost (46 km for Boonah and 30 km for Canungra; see Fig. 1), for all dates with observed foraging activity. We adjusted the total number of points so that an 'available' foraging site was sampled once per 5 km 2 within each of the 25 unique dates (approximately 19000 points total). This sampling design provides a thorough representation of the environmental space at available foraging sites within the foraging radius of the home roosts at times when individuals were roosting there 60 ; however, it results in imbalanced data classes that can lead to underestimation of fitted values and hinder model assessment 62 . We assumed that both 'used' and 'available' foraging areas are equally accessible based on observed foraging radii. So, before modeling, we applied a correction that balances the weight of the two data classes. Following King and Zeng 62 , we weighted 'used' samples as τ = w y / 1 , and 'available' samples as w y Where τ is set to 0.5, and y is the frequency of 'used' foraging areas in the data.

Modeling.
To estimate spatiotemporal foraging intensity, we explored models using boosted regression trees (BRT) with a Poisson error distribution and a log link function. This gives a functional approximation for the expected density of foraging stops within a grid cell at a given time relative to available foraging areas within the foraging radii around our study roosts. We used BRT because regression trees provide an additive model structure that can fit a flexible non-monotonic response that ignores weak covariates 63,64 . The BRT algorithm also reduces the impact of autocorrelation that may result from spatial clustering in the data via 'bagging' , where models at each iteration of the algorithm are built with a randomly selected proportional subset of the data 65 . BRT is especially useful when modeling environmental patterns associated with flying-fox foraging areas because they forage on a wide range of Eucalypt species and forest types, each with a potentially unique pattern of climatically driven flowering. The flexible response fitted by BRT can accommodate the complex relationship between Eucalypt phenology and bat foraging behavior within a framework similar to classic Poisson regression. Moreover, recent literature has elucidated the relative proportionality among Poisson regression, presence-background regression, and inhomogenous Poisson point processes when modeling the intensity of a point process [66][67][68] .
To identify an optimal model, we followed the recommendations of Elith et al. 64 , where we explored the performance of models built with various combinations of meta-parameters such as tree complexity (tc = 5, 7, 9), learning rate (lr = 0.005, 0.001, 0.0005), and bagging fraction (bf = 0.5, 0.6, 0.7). Out of the 27 combinations, the optimal meta-parameter set was determined based on performance on independent test data within a 5-fold spatial cross validation structure. Cross validation folds were assigned using a 5 km spatial grid placed over the study area (defined by maximum foraging radii). Points were assigned to one of the 5 folds based on random assignment of the grid cells into one of the 5 folds (similar to Renner and Warton 67 ; Fig. 14). To ensure reliable metrics in cross validation tests, we selected a random seed that ensured even proportion of used and available foraging areas across the 5 folds. We tested model generality by calculating the mean Poisson deviance and mean squared error (MSE) among the 5 folds. The final selected model was refitted with all data using the meta-parameter combination that returned the lowest mean deviance and MSE, while fitting >1000 trees (as suggested by 64 ). Then we re-fitted the final model with all the data to determine the optimal number of trees on the full data set. To describe the fitted function, we selected 12 influential variables (shown in bold in Table 3) and constructed partial dependence plots by predicting the model over the range of the variable in question while holding all other variables in the model constant. We also explored the influence of interactions among environmental variables and plotted the two most influential interactions. To assess the spatiotemporal pattern of GHFF foraging intensity across the study area, we projected the final model onto the spatial data at 8-day intervals from 2006-2015. We then aggregated predictions to monthly summaries, by taking the mean at each grid cell within each month. We summarized temporal trends in two ways. First, we plotted a time series of violin plots that show the distribution of all cell values for the monthly spatial predictions. Second, we plotted a univariate time series that gives the quartile variation coefficient (QVC; a measure of relative dispersion suited to non-normal count data). Here, high values of the QVC correspond to spatial predictions that produce high values of foraging intensity surrounded by comparatively low values (indicative of high-quality nectar production events), and low values correspond to those with a more homogeneous distribution of low or mid-range values (indicative of poor foraging conditions and dispersed foraging behavior). In this manner, the QVC can be interpreted as a measure of spatial resource concentration.
To assess temporal trends in spatial resource concentration, we plotted model values by year and month and validated them with independent data of population counts. Specifically, we calculated time series of QVC values for spatial predictions within the maximum foraging radius of Canungra (30 km) and compared them to the median value of monthly population counts observed at Canungra from 2008-2014 using a Pearson's product moment test. We focus on the Canungra roost here because over this time period, consistent monthly population counts of Pteropus poliocephalus were performed by the Queensland Department of Environment and Heritage Protection as part of the Southeast Queensland Flying-fox Monitoring Program. Population counts were also performed at the Boonah roost from 2011-2014, however, this roost was subjected to repeated dispersal attempts during this time, leading up to its removal in 2014. We also explored the connection between QVC across the whole study area from 2006-2015, and a broad climate indicator, the Southern Oscillation Index (SOI 69 . In this case, we calculated the cross-correlation function between the temporal summary of the model and the 3-month moving average of monthly SOI values at time-lags of up to 18 months. To explore hypotheses on how eucalypt phenology and observed foraging activity are linked to Hendra virus dynamics, we gathered data on Hendra virus prevalence and spillover occurrence within our study area and plotted seasonal trends with our model of nectar-based foraging resources. We acquired data on Hendra virus prevalence in fruit bats from the field sampling study conducted by Field et al. 26 . As described by Field et al. 26 prevalence values are calculated from pooled urine samples collected under the roost. Although both GHFF and BFF were present at this roost at the time of sampling, Hendra virus excretion is primarily attributed to BFF 15 . Within our study area, Hendra virus prevalence data was only available for the Boonah roosting site. The times and locations of Hendra virus spillover events to horses came from the Queensland Government's public database of documented Hendra virus infection in horses 27 . To summarize seasonal patterns in the data, we fit Generalized Additive Models (GAMs) with cubic cyclic spline smoothing terms to monthly values of spatial resource concentration (approximated with the Quartile Variation Coefficient (QVC)) and Hendra virus prevalence using Gaussian and Quasi-Binomial link functions respectively. The spline terms used integer values corresponding to each month as the covariate. We calculated the total number of monthly spillover events which occurred within three distance thresholds (25,50, and 100 km) of the two study roost sites and plotted all three data sources together. All data processing and modeling was performed in the R statistical language 70 using the mgcv, gbm, dismo, raster, and spatstat packages [71][72][73][74][75] . Data and R code can be found in the online supplementary information.