Shifts in mosquito diversity and abundance along a gradient from oil palm plantations to conterminous forests in Borneo

Deforestation precipitates spillover of enzootic, vector-borne viruses into humans, but specific mechanisms for this effect have rarely been investigated. Expansion of oil palm cultivation is a major driver of deforestation. Here, we demonstrate that mosquito abundance decreased over ten stepwise distances from interior forest into conterminous palm plantations in Borneo. Diversity in interior plantation narrowed to one species, Aedes albopictus, a potential bridge vector for spillover of multiple viruses. A. albopictus was equally abundant across all distances in forests, forest-plantation edge, and plantations, while A. niveus, a known vector of sylvatic dengue virus, was found only in forests. A. albopictus collections were significantly female-biased in plantation but not in edge or forest. Our data reveal that the likelihood of encountering any mosquito is greater in interior forest and edge than plantation, while the likelihood of encountering A. albopictus is equivalent across the gradient sampled from interior plantation to interior forest.


Introduction
Novel pathogens, particularly enzootic, vector-borne viruses, have emerged at an increasing rate in recent decades (Taylor et al. 2001, Huang et al. 2019. Land cover and land use change (LCLUC) are frequently implicated as one of the drivers of this acceleration (Foley et al. 2007, Patz et al. 2008, Lambin et al. 2010, Guo et al. 2019. LCLUC may effect changes in host and vector abundance and diversity, formation of novel vector-host interactions, and shifts in microclimate (Foley et al. 2007, Patz et al. 2008, Lambin et al. 2010, Rogalski et al. 2017, Faust et al. 2018) that engender spillover, the transmission of a pathogen from a reservoir to a novel host (Althouse et al. 2018). To anticipate, and ideally forestall, future emergence events, it is critical to determine the specific mechanisms by which LCLUC may precipitate pathogen spillover, particularly in biodiversity hotspots where the impact of LCLUC is likely to be acute (Han et al. 2016, Marklewitz andJunglen 2019).
Sarawak, Malaysian Borneo (Fig. 1A) encompasses enormous biodiversity and is experiencing unprecedented LCLUC, primarily for oil palm production (Koh et al. 2011, Gaveau et al. 2014. Malaysia is a primary producer of oil palm with over one million hectares of oil palm planted in Sarawak alone (Gaveau et al. 2016, Wong 2018. Insertion of oil palm plantation into forests consistently reduces species richness and is particularly impactful in Borneo, where forests house hundreds of mammal and bird species (The Borneo Project n.d.). Hypotheses for the impact of LCLUC and associated depletion of biodiversity on spillover of arthropod-borne pathogens, reviewed in (Marklewitz and Junglen 2019), range from the dilution effect, wherein spillover potential is lowest at the lowest levels of LCLUC due to high host and vector diversity (Daszak 2000, Keesing et al. 2006, to the intermediate disturbance hypothesis, wherein host and vector biodiversity, interspecies contacts, and spillover potential peak at intermediate levels of LCLUC (Loaiza et al. 2017, Faust et al. 2018, to the amplification effect, wherein spillover potential is highest at high levels of LCLUC due to the increased relative abundance and diversity of susceptible hosts or competent vectors (Keesing et al. 2006, Randolph and Dobson 2012, Faust et al. 2017, Luis et al. 2018. Microclimate also varies substantially between forests and plantations, as temperature is generally higher and more variable (Luskin and Potts 2011, Hardwick et al. 2015, Manoli et al. 2018) and humidity is lower in plantations than forests (Hardwick et al. 2015, Manoli et al. 2018, and plantations vary more in ground cover (Luskin et al. 2017), evapotranspiration, and water runoff than forests (Manoli et al. 2018). All of these differences could impact the life history, vector competence and mobility of arthropod vectors (Reisen 2010, Kramer 2016. Furthermore, millions of people are employed in oil palm agriculture in Southeast Asia, including hundreds of thousands of workers in Sarawak (SALCRA Sarawak Land Consolidation And Rehabilitation Authority n.d., Sarawak Oil Palms Berhad Group of Companies 2018), and local people use the plantations as routes to access the forest (K.I. Young, personal observation).
As pathogen spillover per se is rare, stochastic, and difficult to capture via surveillance, the current study instead investigated the impact of LCLUC on the abundance and diversity of mosquitoes, one of the major vectors of viruses. Specifically, we addressed the hypothesis that oil palm insertion into forests and resulting edges between forests and plantation alters the abundance, diversity, and distribution of mosquitoes in ways that enhance the potential for spillover of sylvatic arthropod-borne viruses (arboviruses), particularly dengue virus (DENV). DENV emerged into the human population on at least four separate occasions from sylvatic cycles maintained in transmission among non-human primates by arboreal, primatophilic Aedes species in tropical forests in Southeast Asia Vasilakis 2009, Vasilakis et al. 2011). Spillover of sylvatic DENV is perilous because each spillover event has the potential to launch new virus strains into human circulation. Unfortunately, several cases of sylvatic DENV spillover into humans have occurred in Borneo in the past decade, all in individuals who had contact with the forest (Cardosa et al. 2009, Teoh et al. 2010, Liu et al. 2016, Pyke et al. 2016).
The impact of oil palm on vector communities has been investigated only to a very limited extent in Borneo (Chang et al. 1997, Brown et al. 2018. Two mosquito species are particularly likely to act as bridge vectors of sylvatic DENV spillover from non-human primates to humans in Borneo: A. albopictus, which is known to feed promiscuously and to vector human-endemic DENV (Hawley 1988, Gratz 2004, and A. niveus s.l., a poorlyknown species complex of primatophilic, canopy-dwelling mosquitoes that is the only known vector for sylvatic DENV maintenance in Asia (Rudnick 1986).
To quantify the impact of oil palm insertion on mosquito communities, we sampled mosquitoes in fourteen oil palm plantations and conterminous forests at different distances from the edge between the two land covers in Sarawak, Malaysian Borneo ( Fig. 1). We used these data to test four specific predictions about mosquito diversity and abundance in oil palm plantations relative to adjoining forests: (1) mosquito diversity and abundance would decline in oil palm plantations, similar to most other taxonomic groups; (2) mosquito diversity would be greatest at the forest edge, consistent with the intermediate disturbance hypothesis (Marklewitz and Junglen 2019); (3) A. albopictus would occur in equal abundance in plantations and forests due to its opportunistic behavior in both bloodfeeding and oviposition choice (Appendix S1: Fig. S1); and (4) A. niveus would be most abundant in forests as it is a canopy-dwelling mosquito (Rudnick 1986; Appendix S1: Fig. S1).

Study area and study sites
The study was undertaken in plantations of oil palms (Elaeis guineensis jacq.) and adjacent forests (Holdridge life-zone classification: Subtropical wet forest; Leemans 1990) within the Bau district and Siburan and Padawan sub-districts of the Kuching administrative division in Sarawak, Malaysian Borneo (central coordinates 1.5533° N, 110.3592° E; Fig. 1A). Insertion of an oil palm plantation into forested land proceeds through three stages, clearance and planting, growth of young tree stands and development into mature tree stands Potts 2011, Dislich et al. 2017; Appendix S1: Fig. S2); a single plantation may contain a mosaic of all three classes. This study was conducted in young stands, albeit old enough to produce fruit (n = 23 sites) and mature stands (n = 113 sites) of 14 different plantations inserted into 10 separate forests. The forests themselves were a mix of secondary (n = 189 sites) and primary (n = 15 sites) growth; none were protected areas. Adult mosquitoes were collected using a backpack aspirator from a total of 340 randomly selected sites at 10 distances from the interface of oil palm and forest, including the edge between forest and plantation (0 m), four distances extending from the edge into interior oil palm (−10, −20, −50, −100 m), and five distances extending into interior forests (+10, +20, +50, +100, +500 m; Fig. 1B).

Land cover classification, sampling design and site selection
We digitized the boundary between oil palm and forest at a scale of 1:3,000 using GoogleEarth satellite imagery from July 25, 2015 and October 2, 2015 in ArcGIS version 10.3.1 (ESRI, Redlands, California). Permission to access plantations was granted from the Sarawak Oil Palm Berhad Group board of directors and forest access was granted from Sarawak Forestry Corporation (permit no.: NCCD.e07.4.4(J LD. 1 3)-299). Oil palm plantations for which permission to access was not granted, or that were <15 hectares, or that lacked adjacent forest were excluded from further analysis, leaving 14 oil palm plantations and 10 adjacent forests.
Isodistances were drawn from the boundary of forest and oil palm (0 m) at five distances into interior forest (+10, +20, +50, +100, and +500 m) and four distances into interior oil palm (−10, −20, −50, −100 m; Fig 1B). Because oil palm is a monoculture, the sampling distance from edge was limited to 100 m; this included 85% of the total sampling area in all oil palm plantations. In contrast, sampling locations were distributed up to 500 m into forests, the maximum distance that was logistically feasible given the rugged terrain, dense vegetation, and the fact that sites were not located along trails. All sampling sites were randomly distributed rather than using transects as mosquitoes may track humans moving along a linear transect. We randomly selected 34 sampling locations at each of the 10 isodistances across all oil palm plantations and forests (N = 340 total sampling sites; Fig.  1B).

Mosquito collections
We chose aspiration as our sampling method, even though an aspirator captures fewer mosquitoes than many other methods, because it is the least biased mosquito collection method available to date (Silver andService 2008, Vazquez-Prokopec et al. 2009). Collections were conducted consistently from June 2016 through January 2017, with an average of seven forest and five plantation sites sampled per working week. Sarawak experiences an equatorial climate where all twelve months receive rainfall (Appendix S1: Fig. S3). Mosquitoes were collected at each site on a single visit beginning no later than 08:00 and ending no later than 15:00. While it was not logistically possible to randomize site visits across all sites, for a given sampling day a random site would be chosen and then a subset of nearby oil palm and forest sites would be sampled the same day to avoid biasing sampling in a given land cover by season. A backpack aspirator (Bioquip, Rancho Dominguez, California, USA) was used to collect mosquitoes from all vegetation up to head height as well as bare ground and leaf litter within a 5 m radius from a sampling point. Sampling was considered complete after the entire area was swept with the backpack aspirator and was not limited to a specific time. Mosquitoes were identified morphologically using taxonomic keys (Walter Reed Biosystematics Unit n.d., Reuben et al. 1994, Huang and Rueda 1998, Rueda 2004, Rattanarithikul et al. 2005, 2006, 2010 to genus and, when possible, species.

Diversity estimation
Alpha diversity measures, Shannon-Wiener index (H) and evenness (H/ln(richness)) (both excluding sites at which no mosquitoes were collected), and genera richness per site (including sites with 0 mosquitoes), were calculated and compared using Kruskal-Wallis tests and pairwise Wilcoxon rank sums tests were used to identify significant pairwise differences. Community dissimilarity was analyzed between distances by calculating the Jaccard index for beta-diversity and the significance of differences among distances were determined using a permutational multivariate homogeneity of group dispersions analysis followed by Tukey's HSD test. All analyses were performed in R version 3.4.4 using the vegan (Oksanen et al. 2018) and multcompView packages using individual sites as the unit of measurement (Graves et al. 2019).

Mosquito abundance and sex ratio analyses
A Pearson chi-square goodness of fit test indicated mosquito abundance data best fit a negative binomial distribution. A generalized linear model (GLM) with negative binomial distribution was then used to assess variation in total mosquito, all A. albopictus, and female A. albopictus abundance at the distances sampled, followed by an adjusted Tukey-Kramer post hoc analysis to identify specific differences among distances. All negative binomial GLM analyses were performed in SAS 9.4 using proc glimmix. Following these analyses, global and local spatial autocorrelation analyses were used to assess spatial patterns in mosquito abundance and occurrence (presence/absence). Moran's I (Moran 1950) was used for the global analysis and Local Indicators of Spatial Association (LISA; Anselin 1995) with false discovery rate (FDR) control was used for the local analysis. Both analyses were run in ArcGIS version 10.7.1 using a 95% confidence level and 999 permutations.
The sex ratio of all mosquitoes combined as well as A. albopictus was analyzed using a contingency table analysis of total counts of female and male across the distances sampled. A Fisher's exact test was used to test whether the observed proportion of all mosquitoes and of female A. albopictus sampled in each land cover was different then an expected ratio of 50:50.

Relationship of climate and vegetation variables to total mosquito and A. albopictus presence
We collected data about multiple environmental variables. Total precipitation for each site sampled was measured remotely via the TRMM 3B43 satellite product. Composite threehour precipitation measurements were gathered for each site at 0.25° spatial resolution for each day 10 days prior to sampling to two days post sampling. Mean temperature and mean percent humidity was recorded every 30 minutes using ibutton Hygrochron data loggers (Maxim Integrated Products, San Jose, California, USA) for three days at each site beginning at sampling and ending 2 days after sampling. The elevation of each site was recorded via GPS unit. The slope of each site was derived from 30-meter spatial resolution SRTM (Shuttle Radar Topography Mission) DEM (Digital Elevation Model) data using ArcGIS 10.6 (ESRI, Redlands, California, USA). Quadrat sampling at each site was used to estimate ground vegetation cover and leaf litter depth. Briefly, five 1 m × 1 m quadrats were set at the center and the approximate four corners of each 10 m diameter site. Vegetation cover was estimated in even categories of 0-25%, 25-50%, 50-75%, and 75-100% total cover. The mode categorization of the five quadrats was used to determine the total ground vegetation cover of the site. At each of these five quadrats, the leaf litter depth was measured using a ruler and the average was calculated. The percent canopy cover was measured from each site by converting photographs of the canopy directly above the center of the site to binary in ImageJ version 1.51 (NIH, Bethesda, Maryland, USA). The total number of open and closed pixels was summed, and the percent of closed pixels was calculated from each picture to represent canopy cover. Soil moisture was estimated from a single sample at each site using the feel and appearance method (U.S. Department of Argiculture 1998). The modified Xu Normalized Difference Water Index (NDWI; Xu 2006), which characterizes water content, was derived from a radiometrically corrected 30-meter spatial resolution Landsat 8 Operational Land Imager (OLI) satellite imagery acquired over the study area on 10 August 2016.
Environmental variables were standardized by subtracting the mean and then dividing by the standard deviation. Correlations between continuous environmental variables were performed using Spearman's rank-order correlation test. Relationships between categorical and continuous environmental variables were evaluated using a non-parametric Kruskal-Wallis test. The pairwise Wilcoxon rank sum was used when significance differences were detected. These tests were then used to assess variation in continuous environmental variables by distance, month, and land cover.
The effect of a subset of environmental variables on total mosquito and A. albopictus presence or absence was tested with a logistic regression with log-link function. A Wald's test was used to assess goodness of fit of the model and the AUC-ROC curve for the model was determined. All analyses were performed in R version 3.4.4 using the MASS (Venables and Ripley 2002), lsmeans (Lenth 2016), Hmisc (Harrell and Others 2019), aod (Lesnoff and Lancelot 2012), and pROC (Robin et al. 2011) packages.
Moran's I (Moran 1950) and LISA (Anselin 1995) with FDR were used as described above to assess global and local spatial autocorrelation in elevation, average leaf litter depth, slope, and NDWI.

Environmental gradients from interior forest to interior oil palm
A complete suite of 11 factors known to influence mosquito occurrence and life history, including temperature, humidity, daily rainfall, elevation, leaf litter depth, canopy cover, slope, elevation, NDWI, soil moisture and ground vegetation cover, were measured at 307 of the 340 collection sites; ibutton data loggers failed at 33 sites. As expected, many of these variables were significantly associated with each other (Appendix S1: Fig. S4 and Table S1) and based on these associations we focused subsequent analyses on temperature, humidity, slope, leaf litter depth and rainfall. Consistent with previous studies of oil palm insertion (Hardwick et al. 2015), mean temperature increased significantly and mean humidity decreased significantly (Fig. 2, Appendix S1: Table S2) and the range in temperature and humidity per site increased (DF = 9, P < 0.0001 for both comparisons) moving from interior forests to interior oil palm. Both temperature and humidity differed significantly by month as well (Appendix S1: Table S3). Leaf litter depth decreased significantly with distance into oil palm (Appendix S1: Table S2) and differed significantly among months (Appendix S1: Table S3). Considering land cover classes, mean daily temperature and leaf litter depth were significantly lower in forest, intermediate at edge and highest in plantation, and mean daily humidity was significantly lower in forest and edge than plantation (Appendix S1: Table S4).

Declines in diversity and shifts in community composition of mosquitoes with increasing distance into oil palm
A total of 883 mosquitoes were collected from 173 sites (51% of 340 total sites: 40 oil palm sites, 133 forest and edge sites) from which at least one mosquito was detected, of which the proportion female was 57%. This sample comprised 19 genera (Fig. 3A) distinguished by morphological identification and one group of eight individuals (<1% of mosquitoes collected) designated as unknown genera; excluding this category did not change the significance of any subsequent analyses. The most frequently collected genus was Aedes (52% of all mosquitoes collected), followed by Culex (17%) and Armigeres (16%). The remaining genera each occurred at <7% of the collection. Seven genera, Aedes, Anopheles, Armigeres, Culex, Topomyia, Uranotaenia, and Verrallina, were collected in both oil palm plantations and forests. The remaining genera were collected only in forests; no genus was collected only in oil palm plantations. Aedes was the only genus collected at all distances sampled.
The number of genera sampled decreased between forest and plantation (Fig. 3A). Importantly, the only genus collected at the most interior oil palm sites, −100 m, was Aedes and 100% of these 9 specimens were A. albopictus (Fig. 3B). Rarefaction curves indicated that the total richness of genera expected was sufficiently characterized only at −20 m into oil palm and +10 m into forests (Appendix S1: Fig. S5), signifying that more genera remain to be detected in the other distance categories. Generally, the diversity of mosquitoes sampled at individual sites (α diversity) was greater in forests than oil palm plantations. Mosquito genera richness was significantly greater at sites sampled within forests compared to oil palm (DF = 9, P < 0.0001; Fig. 4A). Additionally, forests had significantly greater Shannon-Wiener diversity estimates (H) compared to plantations (DF = 9, P < 0.0001; Fig.   4B). Evenness of the mosquito genera sampled also decreased significantly moving from interior forest to plantation, and evenness estimates were significantly greater in forests than plantations (DF = 9, P < 0.0001; Fig. 4C). Diversity was highest at +100 m into forests with a total of 14 genera detected. In all cases, post-hoc pairwise comparisons indicated that diversity at distances sampled within a given land cover type (forest or plantation) were more similar to each other than to the other land cover type, and that edge was not different from either land cover type (Fig. 4).
Mosquito community composition between the distances sampled (β diversity) was compared using Jaccard distances. Communities from plantations and forests differed from each other, but neither plantation nor forest communities differed from edge communities (Table 1; DF = 9, Permutations = 999, P = 0.001).

Mosquito abundance decreased with increasing distance into oil palm
The percent of sites in which at least one mosquito was collected increased at each distance from 15% at −100 m in plantation to 71% at +10 m into forest, after which it levelled off to between 65% and 71% (DF = 9, P < 0.0001). Mosquito presence was significantly spatially clustered at the global level (Moran's index = 0.14, P = < 0.0001). Significant spatial autocorrelation of mosquito presence was detected locally at 9 sites: 5 forest sites with significant Low-Low clustering and 3 forest and 1 oil palm site with significant Low-High clustering (Appendix S1: Fig. S6).
Mosquito abundance, measured as mean number of mosquitoes sampled per site, declined significantly moving from forest into plantation (F(9,330) = 6.86, P < 0.0001). In particular abundance declined sharply and significantly at the deepest distance into oil palm plantation (Fig. 5). Mosquito abundance was significantly higher in June, August, and November than December (F(7,332) = 6.33, P <0.0001). Mosquito abundance was significantly spatially clustered at the global level (Moran's index = 0.04, P = <0.0001). At the local level, spatial autocorrelation of mosquito abundance was detected at 35 sites: 32 with significant low-low clustering (19 in forests and 13 in oil palm), 1 (forest) with significant low-high clustering, 1 (forest) with significant high-low clustering, and 1 (oil palm) with significant high-high clustering (Appendix S1: Fig. S7). There was no significant difference in the proportion of sites with spatial autocorrelation among land cover types (X 2 (2, 340) = 0.16, P = 0.92).
Excluding sites with spatial autocorrelation in the GLM analysis did not change the significance differences among distances in mosquito abundance.
Within plantations, 25 mosquitoes were collected from 35% of the 23 sites sampled in young stands; 102 mosquitoes were collected from 28% of the 113 sites sampled in mature stands. All mosquitoes collected in young stands were Aedes species, while species of Aedes, Anopheles, Armigeres, Culex, Topomyia, Uranotaenia, and Verrallina, were collected in mature stands. Given the small sample size from young stands we did not attempt statistical analysis of these data, and the two age groups are lumped as oil palm in all other analyses.
A logistic regression to test the impact of distance category, month, mean temperature, rainfall, and slope per site on the presence of mosquitoes per site revealed a significant effect only of distance and month on mosquito presence (AUC = 0.793; Appendix S1: Table S5): all distances had a significant positive effect on mosquito presence when compared to −100 m into plantations as a reference, and the proportion of sites in which mosquitoes were captured increased significantly in November compared to June as a reference (Appendix S1: Table S5). We also repeated this analysis replacing mean temperature and rainfall with the average monthly temperature and precipitation for each site from CHELSA (Climatologies at High Resolution of the Earth's Land Surface Areas) satellite data (Karger et al. 2017). Inclusion of these variables in the model recapitulated the impact of distance on mosquito presence, while no other variable had a significant effect independent of distance. Three environmental variables, elevation, slope, and NDWI, clustered significantly at the global level (Moran's index = 0.10, P = 0.0009, Moran's index = 0.06, P = 0.04, Moran's index = 0.09, P = 0.005 respectively), but no significant spatial autocorrelation detected for any of the environmental variables at the local level. Young et al. Page 8

Shifts in mosquito sex ratio between interior forest and interior oil palm
From the perspective of disease transmission, female mosquitoes are of greater importance than males because only females take a bloodmeal. Mosquito sex ratio differed across the distances sampled (DF = 9, P < 0.0001; Appendix S1: Fig. S8A). Female mosquitoes dominated the sample at −100 m in plantations, then the ratio of females declined steadily moving to the forest-plantation edge and rebounded in the forest to ≥50% of the sample. This result was accentuated when analyzing by land cover where the mosquito sex ratio was female biased in oil palm and forest, 64% and 60% of the collection respectively, and male biased at edge (35%; DF = 9, P < 0.0001). Moreover, these ratios differed significantly from a null expectation of 50% female in all three land covers (P = 0.03 in plantation and edge and P = 0.0005 in forest). Sex ratio differed across distances for A. albopictus as well (DF = 9, P < 0.0001; Appendix S1: Fig. S8B), with a pattern generally similar to all mosquitoes combined. The sex ratio of A. albopictus was significantly different from 50% female in plantations (P = 0.01) but not in edge (P = 0.06) or forest (P = 0.38).

Land cover generalism and specialization by A. albopictus and A. niveus, respectively
A total of 339 A. albopictus were collected from 89 sites at all distances sampled, with no statistically significant difference between the mean number collected among distances (F(9,330) = 1.71, P = 0.09; Appendix S1: Fig. S9). Repeating this analysis using the 164 A. albopictus females collected did not change this result (F(9,330) = 1.20, P = 0.30). However, when analyzing these data by land cover type, there was a statistically significant difference in the mean number of A. albopictus collected (F(9,330) = 3.96, P = 0.02), with abundance highest at the edge, lowest in oil palm, and forest having statistically similar abundance to both. When this analysis was limited to just female A. albopictus, this difference was no longer significant (F(9,330) = 0.81, P = 0.45).
A logistic regression was used to test in impact of the environmental variables mean temperature, rainfall, slope per site as well as distance and month on the presence of A. albopictus at individual sites. Distance, rainfall, and slope had significant effects on A. albopictus presence. All distances had a positive effect on A. albopictus presence when compared to −100 m into plantations, with the differences at −50 m and −10 m in plantations, edge (0 m), and 20 m into forests being significant. Only rainfall and slope had significant effects on the presence of A. albopictus. Rainfall was positively associated with A. albopictus presence, while this species was more scarce in more steeply sloped sites (AUC = 0.691; Appendix S1: Table S6).
Only 14 A. niveus were collected, precluding statistical analyses; however, all A. niveus were collected 10 m or more into forests, including two old growth forests and five secondary forests. A. albopictus and A. niveus were collected at the same site only once throughout the study, +500 m into an old growth forest. We also collected 9 other species of Aedes, whose capacity to vector arboviruses is unknown (Fig. 3B).

Discussion
Here we investigated how conversion of forest to oil palm plantation affects the abundance and diversity of mosquitoes in order to gain insight into how this form of LCLUC may impact spillover of mosquito-borne pathogens. We sampled mosquitoes at ten different distances moving from interior plantation to interior forest. These distance categories captured significant variation in microclimate and vegetation: temperature increased and humidity decreased moving from forest to oil palm plantation, as shown in previous studies (Hardwick et al. 2015, Jucker et al. 2018) and leaf litter depth increased moving from plantations into forests. Overall, we found that the shift from native forests to oil palm agriculture drastically reshaped mosquito communities, but we did not find that forestplantation edge supported disproportionately high mosquito abundance or diversity. Importantly, A. albopictus abundance was even across the entire gradient from interior oil palm to interior forest, and A. albopictus was the only mosquito species captured at the deepest distance into oil palm plantations.

Mosquito diversity and abundance in forests and oil palm plantations
We structured this study around four predictions. Our first prediction, that diversity and abundance of mosquitoes would decrease from interior forest to interior oil palm, was borne out. Mosquito diversity dropped dramatically moving from the forest into plantation. Nineteen genera were detected in the forest, while only a subset of these (7 genera) were detected in plantations. Despite the high diversity detected, rarefaction analyses indicated that our sampling did not capture the complete diversity for most distances sampled, similar to other studies of forest and rural habitats in highly biodiverse areas (Thongsripong et al. 2013, Young et al. 2017. Clearly, additional sampling is needed to describe the diversity of mosquitoes across these land cover types. However, as rarefaction curves tended to asymptote or flatten in oil palm and close to the palm-forest edge, we believe that further sampling would amplify the patterns of diversity that we detected in the current study. A caveat to our findings is that we collected only during the day and only at ground level; as mosquito diversity is stratified both vertically and diurnally (Galindo et al. 1950, Brant et al. 2016), a study expanded in time and space would doubtless capture more diversity. However, our study did reflect both the stratum (ground level) and time period (daytime) where humans are most likely to interact with mosquitoes.
Mosquito abundance was even across the five distances sampled within the forest and the forest-plantation edge but declined sharply just −10 m within plantations. This result mirrors other reports of reduction of mosquito abundance in oil palm plantations compared to adjacent forests (Chang et al. 1997, Brant 2011, Zahouli et al. 2017. Two substantial clusters of low mosquito abundance were detected via analysis of spatial autocorrelation. One of these was comprised of mature oil palm and lowland inundated secondary forests sites and the other was comprised of mature oil palm and secondary forests sites at the western end of a highland region. These regions were sampled across different months, and likely represent an underlying, albeit unidentified, environmental variable that disfavors mosquitoes in this region; none of the environmental variables that we measured showed similar patterns of autocorrelation. Spatial autocorrelation decreased significantly when analyzing presence, but the lowland cluster described above of low mosquito abundance persisted. Spatial clustering patterns of vector species is frequently used to assess environmental drivers of vector presence and pathogen transmission risk (Zhou et al. 2007, Diallo et al. 2012, Richman et al. 2018, so future studies of the factors that drives these lacunae in mosquito abundance would be valuable.

Mosquito sex ratio between interior forest and interior oil palm
In addition to the abundance of a given species, its sex ratio, biting rate, host range and vector competence also shape its role in virus transmission. In this study we found that the mosquito sex ratio varied across distances with a female bias in plantations and forests but a male bias at the edge between the two. Chang et al. (1997) reported that human bites by Anopheles species declined as intact forests in Sarawak were replaced with oil palm monoculture, but the biting rate of A. albopictus and C. tritaeniorhynchus remained constant during this shift. Brown et al. (2018) amplified 38 bloodmeals from mosquitoes collected in Borneo from domestic, peridomestic, and forest land covers and identified two hosts, chicken (Gallus gallus) and Homo sapiens; the majority of bloodmeals were from Culex mosquitoes in domestic sites, with only four samples derived from palm and rubber tree plantations and none from forests (Brown et al. 2018). In a separate study, we identified the hosts of five blood-engorged A. albopictus collected within oil palm plantations from the current study and three blood-engorged A. albopictus collected in forests from this study as well as in nearby forests in 2017: all were human (Young et al. 2020). Vector competence of A. albopictus for various arboviruses, and the impact of temperature on competence, has been tested in multiple studies (Alto and Bettinardi 2013, Vega-Rúa et al. 2015, Murdock et al. 2017, Liu et al. 2017, Alto et al. 2018, Reinhold et al. 2018. Generally, competence is thought to have a bimodal relationship with temperature (Mordecai et al. 2017). The effect of temperature on competence has been shown experimentally to depend on the temperature of the larval habitat, diurnal temperature ranges experienced by the adult, thresholds for viability, and interactions among viral genotypes, mosquito genotypes, and environmental conditions (Alto and Bettinardi 2013, Zouache et al. 2014, Vega-Rúa et al. 2015, Murdock et al. 2017, Alto et al. 2018). In the current study, interior oil palm was, on average, 1°C hotter and 4% less humid than interior forest, and the diurnal range in temperature and humidity was significantly greater in plantations than forests. Thus, vector competence of A. albopictus could differ substantially between these land covers. Experimental tests of the competence of other sylvatic mosquito species for arboviruses under a range of microclimate conditions is lacking and would enhance assessment of spillover risk and routes, but laboratory colonization of most of these species has not been achieved.

Mosquito diversity at the edge of forest and oil palm plantations and environmental drivers of mosquito presence
Our second prediction, that mosquito diversity would be greatest at the edge between forest and plantation, was not supported. Diversity at edge was intermediate to forest and plantation, and alpha diversity increased from interior oil palm to interior forests. Moreover, edge communities were not distinct from either oil palm or forest. Similarly, Loaiza et al. (2017) found that mosquito diversity was highest in old growth forests in Panama and declined with increasing forest disturbance. We found no clear environmental driver of mosquito presence across the oil palm and forest interface, suggesting further investigation is required. However, our study was not designed to address potential shifts in mosquito diversity and abundance through time because we did not resample at our study locations. Considering mosquito abundance and presence was significantly different across the months of sampling, conducting a study in which sites were sampled repeatedly across seasons would be informative.

Abundance and distribution of key vector species of sylvatic dengue virus in forests and oil palm plantations
While a large body of empirical , Ostfeld 2017, Marklewitz and Junglen 2019 and theoretical (Luis et al. 2018) studies have investigated the impact of host diversity per se on the prevalence and spillover of vector-borne pathogens, there has been less attention paid to potential effects of vector diversity (but see Marklewitz and Junglen 2019 and studies reviewed therein). Instead, research has generally focused on the abundance and distribution of particular species known to transmit key pathogens. In this vein, we also predicted that two key vectors of sylvatic DENV, A. albopictus and A. niveus, would be habitat generalists and specialists across the distances studied, respectively (Chang et al. 1997, Brant 2011. While the sample size of A. niveus precluded statistical analysis, this species was only captured in forests and at the ground level, despite being described as a canopy-dwelling primatophilic mosquito (Rudnick 1986). This is the second study in Borneo to demonstrate that A. niveus will descend from the forest canopy when enticed by a potential primate bloodmeal at ground level (Rudnick 1986, Young et al. 2017. While this finding suggests that A. niveus could act as a bridge vector in forests, further studies of the vertical stratification of A. niveus are needed. A. albopictus was collected at each of the distances sampled at statistically indistinguishable abundances; however, at the land cover scale we found an increase in abundance at the edge of plantations and forests. Two previous studies have also found a broad distribution of A. albopictus across land covers in Borneo (Young et al. 2017, Brown et al. 2018; in these two studies, the abundance of A. albopictus was slightly but significantly higher in forest edges and agricultural land (not oil palm), respectively, than in other land covers. Importantly, as collections moved deeper into oil palm, the percentage of the collection composed of A. albopictus increased until, by −100 m into plantation, 100% of mosquitoes collected were A. albopictus. Further, the A. albopictus sample from plantation was significantly female-biased. Thus A. albopictus have high potential to act as a bridge vector in oil palm plantations, where people are relatively abundant.

The distribution of non-human primate hosts in forest and agricultural landscapes
The distribution and infection status of hosts is as important as that of vectors for arbovirus spillover. Non-human primates are the hosts of sylvatic DENV in Southeast Asia and Sub-Saharan Africa (Hanley et al. 2013). Recent studies have revealed that non-human primate abundance and diversity is diminished by anthropogenic land cover change, particularly agricultural disturbance (Almeida-Rocha et al. 2017). Oil palm insertion into native forests has been associated with the reduction of non-human primate abundance and diversity in Northern Borneo, Indonesia, the Brazilian Amazon, and Africa (Danielsen and Heegaard 1995, Ancrenaz et al. 2015, Bernard et al. 2016, Mendes-Oliveira et al. 2017, Strona et al. 2018, Seaman et al. 2019. Moreover, the majority of non-human primate detections in oil palm plantations are near forest edges (Ancrenaz et al. 2015, Bernard et al. 2016, Seaman et al. 2019. While non-human primates will sometimes enter interior plantations (Marchal andHill 2009, Ancrenaz et al. 2015), in Borneo only two of the documented 13 Bornean nonhuman primate species (World Wildlife Federation n.d.), orangutans and pig-tailed macaques, have been detected in oil palm plantations (Bernard et al. 2016). Thus, from the perspective of hosts, oil palm plantation would seem to pose a low risk for spillover of sylvatic DENV.

Potential impacts of mosquito community structure on sylvatic DENV spillover in forests and oil palm plantations in Borneo
Based on our findings for mosquito community structure in conjunction with the literature on primate distributions, spillover of sylvatic DENV would be more likely in forests than plantations. A. albopictus occurred at equal abundance in plantation and forest and A. niveus occurred only in this habitat, consistent with previous studies (Rudnick 1986, Young et al. 2017. Most primate species are highly arboreal in Borneo (Phillipps andPhillipps 2016, Bernard et al. 2016); although we did not sample the forest canopy, A. niveus is known to prefer this region (Wells et al. 2004) and A. albopictus has been detected 9 m up in the canopy in Borneo (Brown et al. 2018). Further studies of the vertical stratification of these and other potential vectors are urgently needed to understand changes in spillover risk during tree felling.

Supplementary Material
Refer to Web version on PubMed Central for supplementary material.  Mean and variation in temperature decline while mean humidity increases and variation in humidity declines moving from interior plantation to interior forest: (A) hourly temperatures (°C) and (B) hourly % humidity recorded per site within the distances sampled. See text and Appendix S1: Table S3 for statistical analyses of these data. The relative abundance of mosquitoes collected at 10 distances in oil palm (−100, −50, −20, −10 m), forest edge (0 m), and adjacent forest (+10, +20, +50, +100, +500 m). A total of 340 sites were sampled (n = 34 per distance). (A) The proportion of mosquito genera collected. Samples that could not be identified but were clearly different from identified genera are grouped into a single unknown category. Interior forest had the most unique genera collected, and mosquito diversity tended to decrease with distance into forest edge and oil palm, with concomitant enrichment of Aedes species. (B) The proportion of Aedes species collected. Samples that could not be identified were grouped into a single A. unknown category. Aedes species diversity tends to increase with distance into interior forests and A. albopictus was the dominant species collected at all distances except +50 m into forests. The number of genera (g) and Aedes species (s) detected at each distance are listed above each bar.  The mean number of mosquitoes collected per site per distance differs significantly between oil palm and adjacent forests (F (9,330) = 6.86, P < 0.0001). Statistical analysis of these data is described in the text; significant differences derived from a Tukey-HSD post hoc test are indicated by different letters above bars.