Spatial-temporal dynamics of neotropical velvet ant (Hymenoptera: Mutillidae) communities along a forest-savanna gradient

Understanding how and why biological communities are organized over space and time is a major challenge and can aid biodiversity conservation in times of global changes. Herein, spatial-temporal variation in the structure of velvet ant communities was examined along a forest-savanna gradient in the Brazilian Cerrado to assess the roles of environmental filters and interspecific interactions upon community assembly. Velvet ants were sampled using 25 arrays of Y-shaped pitfall traps with drift fences for one year along an environmental gradient from cerrado sensu stricto (open canopy, warmer, drier) to cerradão (closed canopy, cooler, moister). Dataloggers installed on each trap recorded microclimate parameters throughout the study period. The effects of spatial distances, microclimate parameters and shared ancestry on species abundances and turnover were assessed with canonical correspondence analysis, generalized dissimilarity modelling and variance components analysis. Velvet ant diversity and abundance were higher in the cerrado sensu stricto and early in the wet season. There was pronounced compositional turnover along the environmental gradient, and temporal variation in richness and abundance was stronger than spatial variation. The dry season blooming of woody plant species fosters host abundance and, subsequently, velvet ant captures. Species were taxonomically clustered along the gradient with Sphaeropthalmina (especially Traumatomutilla spp.) and Pseudomethocina more associated, respectively, with cerrado sensu stricto and cerradão. This suggests a predominant role of environmental filters on community assemble, with physiological tolerances and host preferences being shared among members of the same lineages. Induced environmental changes in Cerrado can impact communities of wasps and their hosts with unpredictable consequences upon ecosystem functioning and services.


Introduction
Ecological and historical factors interact in complex ways and operate at different scales during the assembly of biological communities [1][2][3]. At wide spatial-temporal scales, the role of speciation, extinction and dispersion is more prominent, whereas at local scales community assembly is mediated by neutral processes, environmental filters and ecological interactions [4][5][6]. Environmental filters promote the co-occurrence of species that share suitable traits, either due to inheritance from a common ancestor or due to convergence, whereas competition acts in the opposite direction, limiting the coexistence of similar species [7]; because related species tend to be more similar, environmental filters can result in phylogenetic clustering, whereas competition can result in phylogenetic overdispersion [3]. Yet, contemporary coexistence theory indicates that phylogenetic clustering can also result from large fitness differences and competitive exclusion, whereas phylogenetic overdispersion may arise from environmental filters acting on species with similar traits (when traits are convergent) or lack of stabilizing niche differences between closely related species (when traits are conserved) [8,9]. The integration of phylogenies with ecological data obtained at different scales is regarded as essential in contemporary studies on community assembly [10][11][12]. However, this integration is limited by the Darwinian deficit, or the lack of knowledge on the phylogenetic relationships between species, especially in tropical regions [13]. Furthermore, the identification of relevant scales for conducting studies on community assembly is a major challenge [14,15], especially in multitaxonomic approaches, given the great variation in body size, dispersal capacity and generation time between species [16][17][18][19][20].
Since the seminal works of Whittaker [21][22][23], ecological gradient analysis is used to understand the effects of spatial and environmental factors on species distributions and community assembly, as well as community boundaries and species turnover (beta diversity) at different spatial scales [24][25][26]. In addition, the effects of spatial distances and environmental gradients upon beta diversity are relevant to understanding the sensitivity of biological communities to changing environments, such as climate change [27][28][29] or invasive species [30][31][32]. Ecological gradient analysis has become popular with the development of community ordination methods [33] and their modern alternatives [34][35][36], as well as their implementation in a modular open source platform for statistical computing and graphics [37,38].
In this work, we examined the roles of spatial distances, microclimate parameters, and taxonomic relationships upon the structure of velvet ant (Hymenoptera, Mutillidae) communities along a forest-savanna gradient in the Cerrado of central Brazil. Despite their common name, velvet ants are ectoparasitoid wasps that use larvae or pupae of other insects, mainly bees or other solitary wasps, as hosts [39,40]. Therefore, they can regulate host populations and play a fundamental role in ecosystem functioning. Velvet ants comprise about 4200 species distributed mainly in the tropical region [41,42]. They have marked sexual dimorphism, evidenced by the usually winged males as opposed to wingless females, but also reflected in body size, pattern of setae and coloration [43][44][45][46], which makes the taxonomy of the group complicated. The phylogenetic relationships among the major lineages of velvet ants are not well established [47,48] and virtually nothing is known about the relationships among Neotropical species [49,50]. Most Neotropical velvet ants are diurnal, more active early in the morning and late in the afternoon [40,51], and seem to prefer open and dry environments [52][53][54]. Temperature plays an essential role in the development and survival of ectothermic parasitoids, as well as on their hosts and the asynchronous pattern of population fluctuations of both [55]. Therefore, temporal and spatial variation in temperature should be an important environmental filter affecting the abundance and composition of velvet ant communities.
Herein, the following questions were addressed: (1) does the abundance and richness of velvet ants vary predictably throughout the year and along the environmental gradient? (2) what are the relative roles of spatial distance and microclimate parameters on the community structure? and (3) are velvet ants taxonomically structured along the environmental gradient? By answering these questions, the goal was to understand the response of velvet ant communities to spatial and environmental variation, to assess the role of environmental filters and interspecific interactions on the assembly of velvet ant communities, and to discuss the potential impacts of environmental change upon these communities.

Study area
This study was conducted at Parque Municipal Mário Viana, locally known as "Parque do Bacaba" (14˚42'24" S, 52˚21'9" W), in the municipality of Nova Xavantina, Mato Grosso state, Brazil (Fig 1). With 492 ha, it is located in the western portion of the Cerrado biome, close to its contact with Amazonia [56]. The park has forest habitats, such as "cerradão" (dense woodland savanna) or gallery forest, and savanna habitats, such as "cerrado rupestre" (saxicolous vegetation) or cerrado sensu stricto [57,58]. The cerrado sensu stricto is characterized by trees 3-8 m high, abundant herb-grass stratum and a savannic appearance, whereas cerradão is characterized by trees 8-12 m high, closed canopy and a forest appearance [59][60][61]. The regional climate is markedly seasonal, with a dry season from May to September and a rainy season from October to April [62]. The annual precipitation varies between 1300 mm and 1500 mm and the mean annual temperature is 26˚C [63].

Sampling design
The data were collected in a forest-savanna gradient from cerrado sensu stricto (open canopy, warmer, drier) to cerradão (closed canopy, cooler, moister). Twenty-five arrays of pitfall traps with drift fences were installed at every 30 m along a linear transect of approximately 750 m. Each pitfall trap array consisted of four 35 l plastic buckets buried in the ground and arranged as a Y-shape, with the central bucket connected to three peripheral buckets by 6 m long and 50 cm high drift fences made of galvanized steel flat sheets. Traps, numbered 1 to 25 from cerrado sensu stricto to cerradão, were opened for one week every month from August 2015 to July 2016. When opened, traps were checked every day and velvet ants collected with forceps, stored in vials separated by traps, and preserved in 100% ethanol. Subsequently, velvet ants were mounted and identified (many to genus level or morphospecies, including several undescribed taxa). Only female velvet ants were represented in our samples. Fieldwork was conducted under permit SISBIO #335343 from Instituto Chico Mendes de Conservação da Biodiversidade (ICMBio).

Microclimate parameters
Air temperature and relative humidity were recorded at each pitfall trap array, every 10 minutes during the entire study period with dataloggers (HOBO Pro v2 Temperature/ Relative Humidity, Onset Computer Corporation, MA, EUA). Dataloggers were installed two meters away from the central bucket, 50 cm high from the ground and sheltered from rain and direct solar radiation by a PVC hood.
Temperature and relative humidity records from dataloggers were averaged by hour/day/ trap. Subsequently, monthly absolute values (minimum, maximum), means and standard deviations were estimated. This resulted in 14 microclimate parameters for each trap and month (25 traps X 12 months = 300 observations for each parameter): (1) absolute maximum relative humidity (Hmaxa), (2) absolute maximum temperature (Tmaxa), (3) absolute minimum relative humidity (Hmina), (4) absolute minimum temperature (Tmina), (5) absolute standard deviation of relative humidity (Hsda), (6) absolute standard deviation of temperature (Tsda), (7) maximum relative humidity (Hmax), (8) maximum temperature (Tmax), (9) mean relative humidity (Hmean), (10) mean temperature (Tmean), (11) minimum relative humidity (Hmin), (12) minimum temperature (Tmin), (13) standard deviation of relative humidity (Hsd), (14) standard deviation of temperature (Tsd). These parameters represent different microclimate aspects such as its central tendency and variability, as well as dominant monthly patterns and extreme values. A previous study showed that the microclimate is highly correlated with structural habitat and vegetation parameters, and is a good predictor of variation in the structure of velvet ant communities [52]. Due to technical problems with some dataloggers, microclimate parameters were not recorded in September 2015 in trap 25, in October 2015 in trap 25, in March 2016 in trap 6, and in July 2016 in trap 2 (4 in 300 or 1.33% missingness). Missing values can reduce the power of statistical analyses due to loss of information and degrees of freedom. Among the various methods available to estimate missing values [64,65], multiple imputation is considered the most robust because it requires less assumptions about patterns of missing data [66,67] and, in some cases, introduces less bias in analyses than the simple removal of missing data [68]. Missing values were estimated by averaging 100 imputations using package Amelia II [69], which employs an "expectation-maximization with bootstrapping" (EMB) algorithm suitable for time series data.

Statistical analyses
To explore the associations among microclimate parameters, and between traps and microclimate parameters, a Principal Components Analysis (PCA) of the correlation matrix [70] was implemented with package vegan [71]. Species accumulation curves based on samples (pitfall trap arrays) and individuals were calculated to evaluate sampling sufficiency [72] with package vegan. The Chao1 index [73,74] was used to estimate the asymptotic richness of the velvet ant community, while reducing the effects of subsampling and accounting for the presence of undetected species with package vegan. A Correspondence Analysis (CA) was used to explore the associations within species, and between species and traps, with package CA [75]. To reduce the influence of rare species, only those with number of captures higher than average were used [76]. Prior to analyses, species abundance data (number of captures) were transformed to log 10 . To represent the spatial-temporal variation in number of captures and richness, a surface graph was produced by fitting a local regression model using package locfit [77] with the default settings.
To identify the best environmental predictors of changes in community structure along the gradient, two approaches were used: Canonical Correspondence Analysis [78] and a combination of Generalized Dissimilarity Modelling [35] with variance partitioning [79]. The CCA was performed with package vegan and started with a null model, containing only the intersection, followed by the stepwise addition of microclimate parameters considering the Akaike Information Criterion (AIC). Model significance was determined by an analysis of variance (ANOVA) based on 1000 Monte Carlo randomizations. As in CA, only species with number of captures higher than average were used. Because GDM is not as widely used as PCA, CA or CCA in ecological studies, we provide more details about its principles and algorithm. GDM is a novel nonlinear matrix regression method that relates patterns of species turnover between pairs of sampling locations to environmental and geographic distances [80,81]. While CCA associates raw species abundances with environmental (microclimate) parameters, GDM associates dissimilarity in species composition (e.g., beta diversity) to environmental dissimilarity and spatial distance between trap arrays. It has two major advantages over classical linear methods (like matrix regression, CCA), by accommodating (i) variation in the rate of species turnover and (ii) nonlinear relationships between species turnover and environmental/spatial separation. This is achieved (i) by nonlinear transformation of each environmental predictor using maximum likelihood and monotonic I-splines and (ii) by fitting GDM models using generalized linear modelling and specifying adequate link and variance functions [27,35]. GDM was implemented with package gdm [82] using a Bray-Curtis dissimilarity matrix [83] to represent beta diversity, and the default settings.
To assess the existence of taxonomic structure in the distribution of velvet ants along the environmental gradient, a Variance Components Analysis (VCA) was used to describe the hierarchical effects of genera, subtribes, tribes and subfamilies (predictors) upon species preferences to microclimate (response). In the absence of information about phylogenetic relationships between velvet ant species, the taxonomic hierarchy was considered as its best approximation [84][85][86]. Species microclimate preferences were considered as their scores on the first axis of the CCA (above), which correspond approximately to the average of microclimate predictors weighted by species captures in each trap, i.e., their microclimate optima [78]. High components of variance indicate differences in microclimate preferences between taxonomic groups, and phylogenetic conservatism. The VCA was implemented with package VCA [87] through a linear mixed model (LME), where predictors were considered as random effects and model parameters were estimated through restricted maximum likelihood (REML) [88]. All statistical analyses were performed in R [38], using a 5% significance level.

Results
The first two PCA axes (Fig 2) jointly reduced 86.33% of the total variation in microclimatic parameters (PC1 = 60.39%, PC2 = 25.94%). This ordination described an environmental gradient of decreasing air temperature and increasing relative humidity from cerrado sensu stricto to cerradão (Fig 2, S1 Table).
The CCA model selection indicated mean relative humidity and absolute minimum temperature as the best microclimate predictors (F 2,22 = 3.62, P < 0.05), reducing 23% of the variation in species abundances between traps. The CCA ordination of traps and species resembled the CA, revealing strong community structuring along the gradient (Fig 7, S4 Table). Species of Sphaeropthalmina, mostly in the genus Traumatomutilla and especially T. geographica, T. bellifera e T. integella, were strongly associated with traps characterized by lower mean relative humidity and lower absolute minimum temperature (cerrado sensu stricto, traps 1-12), whereas species of Pseudomethocina, especially M. harpyia, H. fronticornis and P. gounellei, exhibited an inverse association (cerradão, traps [13][14][15][16][17][18][19][20][21][22][23][24][25]. The Generalized Dissimilarity Modelling (GDM) indicated that microclimate parameters and spatial distances accounted for 36% of species turnover along the gradient. Spatial distances reduced only 8% of the explained variation (Fig 8) and GDM revealed species turnover across all distances encompassed by the transect (Fig 9A-9C). High covariance was detected between beta diversity explained by microclimatic parameters and spatial distances, corresponding to 44% of the explained variation (Fig 8). Environmental conditions reduced most (48%) of the explained variation in species turnover (Fig 8), the best predictors being absolute minimum temperature, mean temperature, absolute maximum relative humidity and standard deviation of temperature (S5 Table). Species turnover was pronounced at low values of absolute minimum temperature and decreased at intermediate values (Fig 9E). In contrast, there was gradual species turnover for low and intermediate values of standard deviation of temperature, but turnover increased rapidly above this limit (Fig 9D). A continuous species turnover was observed with increasing values of absolute maximum relative humidity (Fig 9F).    The VCA indicated significant taxonomic structure in microclimate preferences above the species level, with virtually all variation concentrated among subtribes ( Table 1). Species of Sphaeropthalmina, all in the genus Traumatomutilla, showed preference for portions of the gradient with lower values of mean relative humidity and absolute minimum temperature when compared to species in other subtribes (Fig 10).

Discussion
Changes in vegetation structure along the environmental gradient from cerrado sensu stricto to cerradão gradient promoted important microclimate shifts at the ground level. Open areas are more exposed to higher solar radiation, being hotter and drier, and experiencing greater variation in temperature and humidity. Besides the spatial variation, there is a pronounced seasonality in the Cerrado, characterized by a regular and predictable variation in precipitation along the year. This spatial-temporal variation in environmental conditions is reflected in the structure of velvet ant communities, where species abundance and richness peak in the warmer and drier portions of the gradient and in the first half of the rainy season, as well as in the patterns of species turnover along the gradient. The different combinations of microclimate parameters can define favourable conditions in small spatial-temporal scales and promote the coexistence of a greater number of species when compared to more constant environments [89]. To a certain extent, this may explain the high species richness recorded in Dynamics of velvet ants in a forest-savanna gradient this study compared to other studies conducted in the Neotropical region [52,[90][91][92][93]]. Yet, these differences can be attributed to variation in sampling methods and effort. For instance, pitfall traps with drift fences seem much more efficient than Malaise or Moericke traps in capturing velvet ants [52,90,94]. Further, the proximity of the study area to the Amazonia-Cerrado transition can provide conditions and resources for species from both biomes [52,95,96]. Ecotones may exhibit greater biodiversity than adjacent areas due to additive processes and the presence of endemic species [97,98], playing an important role in lineage diversification and biodiversity conservation [99,100].
Peaks of abundance in the first half of the rainy season seem to be a pattern among velvet ants in the Neotropical region [90,94]. Because velvet ants are ectoparasitoids of larvae or pupae, mainly of other Aculeata [101][102][103][104], they should show abundance peaks coincident with the end of the provisioning period of their hosts in every generation, which results in Partitioning of total deviance in velvet ant species turnover along an environmental gradient from cerrado sensu stricto to cerradão at Parque Municipal Mário Viana, Nova Xavantina, Mato Grosso, Brazil, according to Generalized Dissimilarity Modelling (GDM) using distance matrices of microclimate parameters and spatial distance as predictors. The stacked histogram represents fractions of total deviance in species turnover explained exclusively by each predictor distance matrix, their shared contribution, as well as deviance left unexplained. https://doi.org/10.1371/journal.pone.0187142.g008 Dynamics of velvet ants in a forest-savanna gradient asynchronous abundance peaks between parasitoids and their hosts [40,105,106]. Accurate data on hosts of Cerrado velvet ants are wanting; however, most Cerrado woody plants are primarily pollinated by animals and bloom at the end of the dry season [107][108][109][110][111], providing plenty of food resources for host populations, which also peak at the same period [90,101,112,113]. Therefore, the peak abundance of velvet ants during the first half of the rainy season occurs shortly after the abundance peaks of their hosts. Yet, because velvet ants are ectotherms, the low temperatures during the dry months can promote a reduction in their activity [114,115].
The open portion of the environmental gradient (cerrado sensu stricto) concentrated greater richness and abundance of velvet ants. Like the temporal variation, this could stem  Dynamics of velvet ants in a forest-savanna gradient from species physiological tolerances to temperature and humidity close to the soil surface, as well as the distribution of host species along the gradient. The greater environmental heterogeneity and the milder microclimate of forest areas provide greater diversity of bees and wasps [116][117][118][119]; nevertheless, velvet ants seem to be more abundant and diverse overall in open areas [52][53][54]. In the Neotropics, this seems to be associated with preference for open areas of two potential host groups, Crabronidae and Sphecidae [90,102,120]. Further, parasitoid abundance and rates of parasitism are often higher in open environments or in edges resultant from habitat fragmentation [121,122]. Also, velvet ants may be better adapted to tolerate variations in temperature and relative humidity, i.e., environments with less vegetation and more exposed soil are suitable for the establishment of these organisms [53]. The strong taxonomic structuring of velvet ants along the gradient indicates a predominant role of environmental filters in community organization. Apparently, physiological tolerances to microclimatic parameters or preferences for certain host species, which determine species turnover along the gradient, are inherited and shared by members of the same lineages (subtribes). The Cerrado, a global biodiversity hotspot, largest and most diverse tropical savanna, sustains accelerated loss of natural habitats due to the expansion of agricultural activities and habitat remnants undergo periodic fires and are subject to edge effects [123][124][125]. Moreover, the biota has been affected by local climate change, characterized by reduced rainfall and increased temperature [126][127][128]. All these factors should contribute to significant changes in the structure of velvet ant communities, especially favouring those lineages more adapted to open habitats (Sphaeropthalmina), to the detriment of the lineages most associated to areas with denser vegetation cover (Ephutini, Pseudomethocina and Smicromyrmina). These changes, in turn, will certainly affect host populations with unpredictable consequences to the trophic cascade and the environmental services they provide.
Supporting information S1