Across sites, tree cover ranged from 0–100% (mean = 27.9%, SD = 34.8%) at the smallest spatial scale we considered, a 30 m radius. Surrounding tree cover within a 1000 m radius, the largest spatial scale considered, ranged from 8.2–74.8% (mean = 33.4%, SD = 20.3%). The lowest and highest site elevations were 12 m and 1451 m above sea level, respectively (mean = 776 m, SD = 473.2 m). Mean annual temperature ranged from 18.8 to 26.4 ºC (mean = 22.6 ºC, SD = 2.3 ºC).
A total of 1,283 mosquitoes representing 48 species in 13 genera were collected from 34 sites (Fig. 2). The number of mosquitoes collected at a site ranged from 0 to 244 (mean = 35, SD = 64). Of these, 108 individuals from 12 residential and four agricultural sites were morphologically identified as Ae. albopictus, and 5 individuals from single sites within each land use category were identified as Ae. aegypti. Ae. albopictus DNA was detected in the pooled samples of molecularly identified mosquitoes from 7 sites where Ae. albopictus individuals were also morphologically identified. The five most common species were Culex quinquefasciatus, Ae. albopictus, Cx. nigripalpus, Wyeomyia adelpha/Guatemala, and Limatus durhamii (Fig. 2, Table S1).
Site-level species richness ranged from one to 19 (mean = 5.2, SD = 4.1). Overall species counts for forest, agricultural, and residential land uses were 33, 29, and 21, respectively. Ten species (21%) were observed in all three land uses. Nineteen species (40%) were shared among forest and agricultural land uses, 13 species (27%) were shared among agricultural and residential land uses, and twelve species (25%) were shared among forest and residential land uses (Fig. 2, Table S1). Eleven species (23%) were found only in forested settings, six species (13%) were found only in agricultural settings, and six species (13%) were found only in residential settings (Fig. 2, Table S1). Two species, Ae. albopictus and Culex quinquefasciatus, were common (observed at > 50% of sites) in residential settings, no species were common in agricultural settings, and three species—Culex nigripalpus, Wyeomyia complosa, and Wyeomyia adelpha/guatemala—were common in forested settings (Table S1).
At least five of the mosquito species observed are known vectors of human diseases. Three of these— the dengue and chikungunya virus vector Ae. albopictus (present at 15 sites) and the St. Louis Encephalitis virus vectors Cx. quinquefasciatus (present at 17 sites) and Cx. nigripalpus (present at 13 sites)—were the three most frequently observed species (Fig. 2, Table S1) (Reisen, 2003; Simmons et al., 2012). Rarely observed vector species included the dengue, chikungunya, yellow fever, and Zika virus vector Ae. aegypti (present at three sites spanning all three land use types) and the malaria vector Anopheles albimanus (present at one agricultural site) (Fig. 2, Table S1) (Zimmerman, 1992; Simmons et al., 2012). In contrast to Ae. aegypti, Cx. nigripalpus, and Cx. quinquefasciatus, which were observed in all land use types, Ae. albopictus was observed only in residential and agricultural settings associated with intensive human modification (Fig. 2, Table S1).
Mosquito species richness was explained by tree cover, but not by land use type. Comparisons of GLMs using tree cover calculated for radii ranging between 30 m and 1000 m surrounding each site indicated that species richness was positively correlated with tree cover at radii between 80 m and 600 m, and tree cover at a 250 m radius had the largest effect size (estimated effect = 1.29 x 10− 2, SE = 4.75 x 10− 3, z-value = 2.71, p-value = 6.65 x 10− 3) (Fig. 3a, Table S2). At the 1000 m spatial scale where both tree cover and climate data were available, the interaction between tree cover and mean annual temperature had a significant effect on species richness (estimated effect = 7.87 x 10− 3, SE = 3.65 x 10− 3, z-value = -2.16, p-value = 3.11 x 10−2 (Table S3). Specifically, at high temperatures, species richness was low even when tree cover was high (Figure S1). By contrast, Kruskal-Wallis test results indicated that species richness did not differ significantly among forested, agricultural, and residential sites (chi-squared = 2.60, df = 2, p-value = 0.28) (Fig. 3b). Notably, the highest species richness was observed at a site in the Coto Brus forest reserve (Fig. 3b).
From the Aedes disease vector survey, we present results only for Ae. albopictus because observations of Ae. aegypti were insufficient for statistical analysis. Both tree cover and land use type predicted Ae. albopictus presence. Comparisons of GLMs using tree cover calculated for radii ranging between 30 m and 1000 m surrounding each site indicated that Ae. albopictus presence was negatively correlated with tree cover at radii between 50 m and 200 m, and was best explained by tree cover at a 110 m radius (estimated effect = -4.23 x 10− 2, SE = 1.96 x 10− 2, z-value = -2.16, p-value = 3.08 x 10− 2) (Fig. 4a, Table S3). At the 1000 m spatial scale where we additionally assessed the influence of climate, Ae. albopictus presence was negatively correlated with tree cover and positively correlated with temperature (tree cover estimated effect = -7.76 x 10− 2, SE = 3.74 x 10− 2, z-value = -2.07, p-value = 3.81 x 10− 2; mean annual temperature estimated effect = 6.62 x 10− 1, SE = 3.35 x 10− 1, z-value = 1.97, p-value = 4.814 x 10− 2) (Table S5). Land use type also predicted Ae. albopictus presence (Kruskal-Wallis chi-squared = 9.58, p-value = 8.311 x 10− 3). Specifically, Ae. albopictus was significantly more likely to be observed in residential settings (present at 13/17 sites) than in forested settings (present at 0/8 sites) (Kruskal-Wallace chi-squared = 9.03, Bonferroni-adjusted p-value = 7.98 x 10− 3), and its presence in agricultural settings (present at 4/12 sites) did not differ significantly compared to either residential (Kruskal-Wallace chi-squared = 2.41, adjusted p-value = 3.6 x 10− 1) or forested (Kruskal-Wallace chi-squared = 3.04, adjusted p-value = 2.43 x 10− 1) settings. Fourteen of the 15 sites where Ae. albopictus was present were surrounded by less than 25% tree cover within a 50 m radius. The only site within a pine plantation was a clear outlier, where Ae. albopictus was present under 85% tree cover (Fig. 4B).
In contrast to species richness, community composition and dispersion were predicted by land use type. PERMANOVA results comparing all three land uses showed that land use significantly affected community composition (Sum of squares = 2.113, R2 = 0.16, F-value = 3.03, p-value = 0.001), and pairwise PERMANOVAs showed that agricultural and forest mosquito communities differed significantly from residential communities (Table 1). Wider dispersion among agricultural compared to residential mosquito communities (average distance to the median: agriculture = 0.612, forest = 0.572, residential = 0.497; Tukey test adjusted p-values: residential – agricultural = 0.0422, residential – forested = 0.303, forested – agricultural = 0.736) likely contributed to the community dissimilarity detected between these land uses (Anderson and Walsh, 2013) (Fig. 5). NMDS visualization of communities grouped by land use type suggests that more variable agricultural communities bridge relatively distinct forest and residential communities (NMDS stress = 0.12), in agreement with the statistical test results (Fig. 5, Table 1). The wider variation among agricultural sites is also evident from the species observation table: no single species was observed at more than 1/3 of all agricultural sites, whereas Ae. albopictus and Unidentified Culicidae 1 were both observed at > 70% of residential sites, and Culex nigripalpus, Wyeomyia adelpha/guatemala, and Wyeomyia complosa were each observed at > 60% of forested sites (Table S1).
Table 1
PERMANOVA results for community composition compared among land use types. Asterisks indicate p < 0.05.
Land use pair
|
Sum of squares
|
F-Value
|
R2
|
Bonferroni-adjusted p-value
|
agriculture vs. forest
|
0.591
|
1.45
|
0.0785
|
0.291 *
|
agriculture vs. residential
|
0.788
|
2.22
|
0.0846
|
0.021 *
|
forest vs. residential
|
1.56
|
4.63
|
0.181
|
0.003 *
|
Finally, generalized dissimilarity modeling (GDM) indicated that environmental gradients explained little of the species turnover among sites. Among the spatial scales for which tree cover was calculated, the model using the 50m radius explained the highest amount of species turnover among sites (Table S6). The model that included mean annual temperature, geographic distance, and tree cover at the 50m radius explained 7% of species turnover among sites. Elevation showed no relationship with species turnover. Whereas increasing tree cover was associated with a consistent increase in community turnover, increasing temperature was associated with a steep increase in community turnover up to a plateau around 22 ºC (Fig. 6).