Do an ecosystem engineer and environmental gradient act independently or in concert to shape juvenile plant communities? Tests with the leaf-cutter ant Atta laevigata in a Neotropical savanna

Background Ecosystem engineers are species that transform habitats in ways that influence other species.While the impacts of many engineers have been well described, our understanding of how their impact varies along environmental gradients remains limited. Although disentangling the effects of gradients and engineers on biodiversity is complicated—the gradients themselves can be altered by engineers—doing so is necessary to advance conceptual and mathematical models of ecosystem engineering. We used leaf-cutter ants (Atta spp.) to investigate the relative influence of gradients and environmental engineers on the abundance and species richness of woody plants. Methods We conducted our research in South America’s Cerrado. With a survey of plant recruits along a canopy cover gradient, and data on environmental conditions that influence plant recruitment, we fit statistical models that addressed the following questions: (1) Does A. laevigata modify the gradient in canopy cover found in our Cerrado site? (2) Do environmental conditions that influence woody plant establishment in the Cerrado vary with canopy cover or proximity to A. laevigata nests? (3) Do A. laevigata and canopy cover act independently or in concert to influence recruit abundance and species richness? Results We found that environmental conditions previously shown to influence plant establishment in the Cerrado varied in concert with canopy cover, but that ants are not modifying the cover gradient or cover over nests. However, ants are modifying other local environmental conditions, and the magnitude and spatial extent of these changes are consistent across the gradient. In contrast to prior studies, we found that ant-related factors (e.g., proximity to nests, ant changes in surface conditions), rather than canopy cover, had the strongest effect on the abundance of plant recruits. However, the diversity of plants was influenced by both the engineer and the canopy cover gradient. Discussion Atta laevigata in the Cerrado modify local conditions in ways that have strong but spatially restricted consequences for plant communities. We hypothesize that ants indirectly reduce seedling establishment by clearing litter and reducing soil moisture, which leads to seed and seedling desiccation. Altering soil nutrients could also reduce juvenile growth and survivorship; if so these indirect negative effects of engineering could exacerbate their direct effects of harvesting plants. The effects of Atta appear restricted to nest mounds, but they could be long-lasting because mounds persist long after a colony has died or migrated. Our results support the hypothesis that leaf-cutter ants play a dominant role in Cerrado plant demography. We suggest the ecological and economic footprint of these engineers may increase dramatically in coming decades due to the transformation of the Cerrado by human activities.


INTRODUCTION
Species that transform habitats or create new ones are known as Ecosystem 54 Engineers (Jones et al. 1994;Jones et al. 1997), and they can have major impacts on 55 population dynamics, community composition, and ecosystem function (reviewed in 56 Kleinhesselink et al. 2014;Wright & Jones 2006). Most research on engineers to date 57 has focused on documenting the magnitude of their impacts on local biodiversity, with 58 more recent work evaluating how these impacts vary spatially (e.g., Badano et al. 2006;59 Baker et al. 2013;Dibner et al. 2015;Kleinhesselink et al. 2014;McAfee et al. 2016). An 60 emerging area of interest is identifying how the impacts of engineers vary along or even   88 populations and communities -they are major seed predators and harvest seedlings to 89 use as the substrate for their fungal gardens (e.g., Costa et al. 2017;Silva et al. 2012;90 Vasconcelos & Cherrett 1997). In addition to their direct impacts on plants, however, 91 their alteration of soil properties throughout of the landscape may also indirectly 92 influence plant growth, survivorship, or community composition (e.g., Bieber et al. 2011;93 Garrettson et al. 1998;Meyer et al. 2011;Sosa & Brazeiro 2012). To date the potential 94 for engineering by Atta to indirectly influence plant communities has primarily been 95 studied in lowland tropical forests (Farji-Brener & Illes 2000;Leal et al. 2014). However, 96 the abundance of Atta colonies can be 2-3 fold greater in the Cerrado (Costa & Vieira-97 Neto 2016), where they have the ability to completely defoliate trees (Mundim et al. 98 2012). This suggests a novel means by which this ecosystem engineer could indirectly 99 shape plant diversity -by modifying canopy cover gradients, and therefore the local 100 environmental conditions that influence seedling establishment and the subsequent 101 growth of juvenile plants (Correa et al. 2016). The magnitude of these indirect impacts 102 should vary along the gradient, however, because areas where trees are sparse will 103 already be hotter, brighter, and have limited litter on the soil surface.

104
To elucidate how gradients and ecosystem engineers interact to influence plant 105 biodiversity we used data on the distribution of over 1200 plants in a Cerrado landscape 106 dominated by the leaf-cutter ant Atta laevigata. Our study addressed the following 107 questions: (1) Does A. laevigata modify the gradient in canopy cover found in our 108 Cerrado site? (2) (Cardoso et al. 2009;Appendix A). Cerrado ralo has a dense 122 layer of grasses and herbs and sparsely distributed shrubs and trees typically <3m tall; 123 the average canopy cover in cerrado ralo is ~30%. Cerrado denso has less grass cover 124 and more abundant trees that can reach a height of ca. 6 m; average canopy cover in 125 cerrado denso is ~60%. There is large variation in the canopy cover of both vegetation 126 types, however, so there can be strong gradients in canopy cover in landscapes where 127 they abut. At Panga Station, for instance, the canopy cover gradient in the Cerrado ralo 128 /Cerrado denso mosaic ranges from 0-95% (Mean = 52% ± 33.1 SD; Figure 1A).

129
Our focal ecosystem engineer is Atta laevigata, whose large nest mounds are 130 formed by workers depositing excavated soil around the main entrance to the nest. Atta 131 laevigata is the most common Atta species in both cerrado ralo and cerrado denso 132 (Costa & Vieira-Neto 2016); A. sexdens is also found at Panga Station, but primarily in 133 closed-canopy forest. In 2010 we haphazardly selected 10 active A. laevigata nests in 134 each vegetation type (surface area of the N=20 nests: 7-37 m 2 , mean = 16.7 m 2 ± 6.7 135 SD). We then established three 1x2 m plots around each nest -one on the center of the 136 nest mound, one immediately adjacent to the mound, and one 10 m from the mound 137 edge (Appendix A). In each plot, we recorded the identity and height of all woody plants 138 <1.2 m tall as well as data on environmental conditions (described below). Although 139 there are some abandoned nests in our site, we restrict our analyses to active colonies 140 because the effects of time-since-abandonment on environmental variables is unknown.

142 Environmental data
Litter biomass in each plot was measured by collecting all litter from a randomly 144 selected half of each plot once during the 2010-2011 rainy season, drying it at 50° C for 145 72 h, and weighing it with a microbalance. Similarly, we dried and weighed all grasses 146 from a randomly selected half of each plot to estimate above-ground grass biomass. 147 Canopy cover above each plot was estimated using photos analyzed with Adobe 148 Photoshop (Adobe Systems Inc., San Jose, California, USA). In our analyses we used 149 the average canopy cover in two photos taken during the same rainy season. Photos 150 were taken with a Nikon Coolpix 950 from a height of 50 cm (modified from Engelbrecht 151 & Herz (2001) in either the early morning (6h) or early evening (18h).

152
At the end of the 2011 dry season we estimated surface soil moisture content in 153 plots by collecting a sample of the top 20 cm of soil from two points separated by 100 154 cm. These samples were bulked, weighed, dried at 50° C for 96h, then weighed again 155 to estimate percent moisture content. As a proxy for soil compaction we used soil 156 penetrability: we dropped a 1 m long x 5 mm diameter iron rod vertically from a height of 157 50 cm at three haphazardly selected points in each plot, then measured the depth to 158 which the rod penetrated the soil at each point. We used the average of these values in 159 our analyses; these data were recorded at the end of the 2010-2011 rainy season.  During the 2011 rainy season we selected N = 5 nests in Cerrado ralo and N=5 165 nests in Cerrado denso nests for analyses of soil chemistry. For each nest we collected 166 soils in the plots on nest mounds and the plots 10m from nests. We collected 5 soil 167 samples of ~100 g samples each from each plot: one from the plot center and one from 168 each corner. The 5 samples from each plot were bulked into a single sample and taken 169 to the Soils Analysis Lab of the Federal University of Uberlândia (UFU), where pH, P, 170 K + , Ca 2+ , Mg 2+ , Al 3+ , and total organic matter were measured using their standard 171 protocols (EMBRAPA 1997).  We used Principal Components Analyses (PCA) to summarize environmental 186 conditions in each plot because many of the biophysical variables we measured were 187 highly correlated (Appendix B). The complete suite of environmental data was only 188 collect in a subset of N = 10 nests, so we conducted two separate PCAs. The first was Statistical analyses: Do A. laevigata and canopy cover act independently or in 214 concert to influence plant abundance and species richness? 215 We used two sets of Generalized Linear Mixed Models with Poisson error 216 distributions to determine if plant abundance and species richness in plots were best 217 explained by ant-related factors (proximity to leaf-cutter ant nests, environmental 218 conditions in plots), canopy cover, or a combination of both. The two analyses that were 219 identical except for the PCA scores used to summarize local environmental conditions: 220 the first group of models used the axis scores from 'PCA-1' (i.e., all nests and plots but 221 no data on soils), while the second used the axis scores from 'PCA-2' (all environmental 222 variables but fewer nests and plot locations). Juvenile plant abundance or richness were 223 the dependent variables. Once again preliminary analyses indicated including colony 224 area did not improve the fit of models. Nest identity was again included as a random 225 effect; because of significant overdispersion the models for plant abundance also 226 included a random per-observation term.

232
We surveyed N=1257 individual plants from N=66 genera. We were able to 233 identify 89% of these stems, with the remainder assigned to one of N=35 234 morphospecies (Appendix C). Most individuals were trees (62.6%) or shubs (26.6%), 235 with only 10.8% of stems being woody vines. The most common species recorded were 236 Miconia albicans (Melastomataceae, N=239), Tapirira guianensis (Anacardiaceae, 237 N=98), Matayba guianensis (Sapindaceae, N=66), Cordiera myrciifolia (Rubiaceae, 238 N=65) and Serjania erecta (Sapindaceae, N=57); these five species represent 42% of 239 the stems sampled (Appendix C). Average plant height was 16.6 cm ± 20.9 SD; 83% of 240 the stems were <30 cm tall and 56% were <10 cm tall, suggesting most were seedlings 241 or relatively recent recruits. The model that best fit the data on the amount of canopy cover over a plot is the 245 one including only the random effect of nest identity (Table 1). This indicates A. 246 laevigata colonies alter canopy cover around their nests, but not in a systematic way 247 related to nest size, and that there is no predictable change in canopy cover as a 248 function of proximity to ant nests ( Figure 1B). Plots on nest edges and those far from nests overlapped in ordination space, 253 indicating they had very similar environmental conditions ( Fig. 2A). However, there was 254 almost no overlap in ordination space between either of these locations and the plots 255 located in the middle of A. laevigata nest mounds ( Fig. 2A), even when the number of 256 nests was reduced to include soil data (Fig. 2B). In 'PCA-1' the first axis explained 257 45.6% of the variance and was positively correlated with litter biomass and soil moisture 258 content. The second axis explained an additional 29.6% of the variance; it was 259 negatively correlated with grass biomass and soil penetrability (Table 2). In 'PCA-2' the 260 first axis explained 42.9% of the variance and was positively correlated with litter and 261 grass biomass, soil moisture content, and soil P, Al 3+ , and organic material (Table 3). 262 The second axis explained 21.4% of the variance and was positively correlated with all 263 other environmental variables measured. In light of these results we used the scores 264 from the first principal components as the dependent variable in subsequent analyses.

265
When using the results of 'PCA-1', canopy cover over a plot was positively 266 correlated with a plot's PCA1 score (ρ = 0.44, t = 3.77, df = 58, p < 0.001), suggesting 267 an association between canopy cover and local environmental conditions. However, the 268 best-fit model included both canopy cover and plot location. This indicates leaf-cutter 269 ants also influenced environmental conditions, but that the magnitude of the effect 270 varied with plot proximity to nests (Table 4, Fig. 3A). When data on soils were included 271 in the PCA, however, there was no longer a correlation between canopy cover over a 272 plot and that plot's score on the 1 st PCA axis (ρ = 0.30, t = 1.35, df = 18, p = 0.20).
273 Furthermore, the model that best fit the data on environmental conditions in a plot only 274 included the proximity of a plots to ant nests and the random effect of nest identity 275 (Table 5, Fig. 3B). In other words, the impact of ants on environmental conditions 276 influencing establishment far outweighs that of canopy cover, but this is only revealed 277 once data on soil properties are included in the analyses. We found 20.95 ± 18.14 SD recruits (range=0-85) in each 2 m 2 study plot. 282 However, the mean number of recruits plot -1 decreased as one moved closer to the 283 center of nests: plots far from nests had on average 29.55 ± 19.39 SD recruits in them 284 vs. 24.9 ± 15.62 SD recruits plot -1 on nest margins and 8.4 ± 11.93 SD recruits plot -1 in 285 the center of nest mounds. The mean number of species per plot was also lowest in 286 plots on the center of nests (3.2 ± 2.9 SD) with almost four-fold higher species richness 287 in plots on nest margins (10 ± 4.4 SD) and 10 m from nests (11.8 ± 4.6 SD).

288
Plant abundance in plots appears to be primarily influenced by ant-related factors 289 ( Figure 4A,B). The best fit models include the proximity of plots to ant nests and 290 environmental conditions in plots (Tables 6,7), though note dAIC values for the model 291 including canopy cover was <2. Species richness in plots, however, appears to be 292 shaped by both canopy cover and Atta-related effects ( Figure 4C,D). These results 293 were consistent whether the PCA used to summarize environmental conditions in plots 294 included data on soils or not (Tables 6,7). The significant effect of nest identity also 295 indicates that some nests exert larger or smaller effects on local the abundance and 296 diversity of recruits than others of similar size. Both ecosystem engineers and environmental gradients are known to exert 300 strong effects on biodiversity, but it is unknown if in general they act independently or in 301 concert. This is because empirical studies simultaneously evaluating the relative 302 influence of engineers and gradients remain rare (e.g., Badano & Marquet 2009;303 Kleinhesselink et al. 2014). We quantified the abundance and diversity of woody plant 305 gradient, as well as data on environmental conditions influencing plant establishment, 306 growth, and survivorship (Fig. 2). Our results suggest that plant diversity in plots is 307 shaped by both leaf-cutter ants and canopy cover. However, seedling abundance in 308 plots is primarily driven by the ecosystem engineer, which both harvests plants and 309 alters demographically relevant environmental conditions.

337
It is notable that the impacts of A. laevigata on the abundance and diversity of 338 plant recruits appears restricted primarily to the nest mound itself, which may limit the 339 spatial extent of an individual colony's impact. However, a salient feature of many 340 engineers is that their localized impacts can often persist long-term (Hastings et al.     Generalized Linear Mixed Model selection for the effect of plot proximity to leaf-cutter ant (Atta laevigata) nests on canopy cover in plots The significance of plot proximity was assessed by comparing the model including only the random effect of nest identity (model 1) with models including this random effect, plot proximity to ant nests, and nest mound area as a covariate (model 2: no plot location x covariate interaction; model 3: main effects of plot location, the covariate, and a plot location x covariate interaction). All models used a Gaussian distribution with an identity function; nest mound area was not included as a covariate because preliminary analyses indicated it did not improve the fit of models. Considering the location of plots or nest mound area does not improve the fit to the data, indicating canopy cover is independent of proximity to ant nests and nest mound size. The best model is noted in bold.

Model Factors
Resid. Df Resid.  Factor loadings for the four principal components axes summarizing environmental variables measured in study plots located in Brazilian Cerrado.
The cumulative proportion of the variance explained by these axes = 100%. The variables included in this PCA (referred to as PCA-1 in the text) were litter biomass, soil penetrability, grass biomass, and soil moisture content. Data for PCA-1 were collected in plots on the center of, adjacent to, and 10m from the edge of N = 20 all Atta laevigata nest mounds.  Factor loadings for the first four principal component axes summarizing environmental variables measured in study plots in Brazilian Cerrado.
The summed proportion of the variance explained by these axes is 84.9%. The variables included in this PCA (referred to as PCA-2 in the text) were litter biomass, soil penetrability, grass biomass, soil moisture content, soil pH, several soil macronutrients, and soil organic material, and the data were collected in plots in the center of and 10m from N = 10 Atta laevigata nest mounds.     Model selection for effects of canopy cover vs. nest mound area, plot location, and local environmental conditions (i.e., PCA-2 scores) on seedling abundance and species richness.
The significance of these factors was assessed by comparing the models including only the random effect of nest identity and per-observation random effects (model 1) with models including these random effects and canopy cover (model 2), random effects and those related to ants (model 3), or random effects and canopy-cover and ant-related variables, and local environmental conditions (axis 1 scores from PCA-2), which analyses indicated were influenced by both canopy cover and proximity to ant nests (model 4  . Canopy cover over plots on Atta laevigata nests (blue circles), adjacent to nests (blue triangles), and far from nests (gray squares). Canopy cover is independent of plot proximity to the N=20 nests (Table 1) Figure 1 (B). Canopy cover over plots on Atta laevigata nests (blue circles), adjacent to nests (blue triangles), and far from nests (gray squares). Canopy cover is independent of plot proximity to the N=20 nests (Table 1), indicating ants are not responsible for or modifying the canopy cover gradient in our study site. The variables included in this PCA (referred to as PCA-2 in the text) were litter biomass, soil penetrability, grass biomass, soil moisture content, soil pH, several soil macronutrients, and soil organic material. Symbol size indicates the percent canopy cover over that plot.