Alternative stable ecological states observed after a biological invasion
Study systemOur focal ecosystem is in Selvíria, state of Mato Grosso do Sul, Brazil ((hbox {20}^{circ }) (22′) (41.86”) S, (hbox {51}^{circ }) (24′) (58.90”) W), on a property owned by the São Paulo State University (UNESP). The location covers 350 ha of pasture composed of liverseed grass (Urochloa decumbens). The native vegetation was removed, pasture areas were implemented, and livestock was introduced in the 1970s, maintaining this configuration during the following 50 years. The climate of this area is categorized as equatorial savanna, with dry periods concentrated mostly during the winter, from April to August. During our sampling period (from November 23th, 1989, to November 19th, 2015), no vermifuges and insecticides that could affect negatively the community of dung beetles associated with cow pads were used1.The native dung beetle community at this site is composed of dwellers and tunnelers. Dwellers comprise the Aphodiinae subfamily, whereas all the tunnelers belong to the Scarabaeinae subfamily31. In total, there were eight species classified as dwellers (Ataenius crenulatus, A. picinus and Atanius aequalis-platensis grouped as one species, Blackburneus furcatus, Genieridium bidens, Labarrus pseudolividus, Nialaphodius nigrita and Trichillum externepunctatum) and ten native tunnelers (Ateuchus nr. puncticollis, A. vividus, Canthidium nr. pinotoides, Dichotomius bos, D. semiaeneus, D. sexdentatus, Ontherus appendiculatus, O. dentatus, O. sulcator). These species were chosen for our study because, as the invasive tunneler D. gazella (also from the Scarabaeinae subfamily), they all co-occur in pasture and exploit the same resource (cow pad)32. The initial establishment of D. gazella caused the loss of most of the native tunnelers from the community, with the invader becoming the overwhelming representative of the functional group, and an initial decrease of abundance for dwellers. Differently from native tunnelers, however, dwellers were able to recover their number a few years after invasion (Fig. 1a, Fig. S1).As reported in1, the abundance of dung beetles was significantly affected by both local minimum temperature and relative humidity. The influence of these two factors is expected, as they determine egg and larval survival and development of dung beetles. For example, because dung beetles are poikilotherms, environmental temperature is key to their development and fecundity33. One of the main dweller species, Labarrus pseudolividus, is widely found in locations with temperature averages ranging between (hbox {12},^{circ }hbox {C}) and (hbox {18},^{circ }hbox {C})34, making it tolerant to colder local temperatures. On the other hand, for D. gazella the lower developmental threshold is (hbox {15.5},^{circ }hbox {C}) (individuals cannot survive below this temperature), and the optimum temperature for population growth is (hbox {28},^{circ }hbox {C})35. For both groups, physiological growth and reproduction rates are maintained even when outside temperatures are close to the lower developmental threshold; dwellers, for example, live inside the dung pile, where temperature is higher and less variable than outside36,37. However, while tunnelers oviposit deep in the soil to protect the eggs, warmer and drier conditions reduce dweller egg viability on dung piles since they are exposed38. Low humidity conditions lead to drier dung and can cause egg and insect dessication. In addition, dwellers from our focal system have Palearctic evolutionary origins39; D. gazella’s natural distribution ranges from central to southern Africa40, presenting high physiological plasticity that allows it to tolerate high temperatures and low relative humidity better than other tunneler species41.Functional-group data collection and community structure characterizationDung beetles were collected once a week in a black-light flight intercept trap42, which guarantees the collection of coprophagic beetles. During all collection periods, climate variables were also collected from a meteorological station located within 2 km of our collecting site. See1 for the complete description of the collection process and database. For our purposes, we retained the species, number of individuals per species, and climate variables for each week sampled (Supplementary Information, SI, Figs. S1–S2).We focused first on the weekly abundance data, which we needed to process in order to avoid spurious results in our analyses stemming from the measurement protocol. Specifically, we filtered out seasonal low values associated with sampling in the coldest periods, when few beetles are captured because the reduced activity in all functional groups restricts their spatio-temporal distribution43. Including such samples would not be representative of the community and could bias the analysis since we are investigating community composition (i.e. proportions, very sensitive to low sampling). Thus, we considered only samples with a total number of beetles (that is, summing up all groups together) higher than the value of the median of all data, a conservative threshold that retains observations that allow for as much representation of the community as possible. As will become evident in the Results section and Supplementary Information, less conservative choices for the threshold did not alter our main conclusions.Following Mesquita -Filho et al.1, we categorized all sampled species into either dwellers or tunnelers. D. gazella is a tunneler and, as explained above, the native tunneler species experienced massive declines in abundance after its establishment, leaving D. gazella as almost the single representative in the tunneler functional group during the period of observation1. Thus, given the sharp contrast in community composition, we also separated the data into before and after invasion using to that end the 200th week, when D. gazella was first observed at the study site (September 11th, 1993, starting date for what we will call “after invasion”, our focal period henceforth).To describe community functional composition (i.e. system state) through time, we derived a normalized functional group ratio. First, because the abundance of each functional group spanned up to four orders of magnitude, we performed a logarithmic transformation of the number of captured insects from each group i, (log _{10}(N_{i}+K)), following Yamamura44. Here, we chose (K=1), but the value of K did not alter our results qualitatively. In addition, the original data showed random mismatches in the phenology of each group, which gave the wrong impression of extreme short-term shifts in functional group dominance within the community. To avoid such artifacts, we used nonparametric local regression (LOESS)45 to smooth the dynamics of each group46. For this smoothing, we employed the loess function in the R software 3.6.147 with a smooth parameter equal to 0.25, but other moderate values (or an optimal value calculated with Bayesian inference by the R function optimal_span) did not alter our conclusions. Finally, we extracted back from the smoothed curve the number of beetles within each functional group to calculate the fraction (f_{dwell}) that measures the relative abundance of dwellers:$$begin{aligned} f_{dwell} = frac{N_D}{N_D+N_T} end{aligned}$$
(1)
where (N_D) corresponds to the number of dwellers per week and (N_T) corresponds to the number of native tunnelers (for the period before invasion), or only the number of D. gazella observed per week (after invasion), using their corresponding smoothed curves. Including also native tunnelers after invasion did not alter our conclusions.Climate driverWe devised a single climatic driver variable that merges the weekly measurement of temperature and relative humidity over the years, abiotic factors key to the survival and reproduction of both groups (see above). We first converted minimum temperatures and relative humidity to normalized climate variables using a min-max normalization (a feature scaling that uses the total range of temperatures or relative humidity, respectively, as normalization factor):$$begin{aligned} T = frac{T_{week} – T_{min}}{T_{max}-T_{min}};;,~ ~ ~ ~ ~ ~ RH = frac{RH_{week} – RH_{min}}{RH_{max}-RH_{min}};;, end{aligned}$$
(2)
where T corresponds to the normalized temperature, (T_{week}) is the weekly temperature, and (T_{max}) and (T_{min}) are the absolute maximum and minimum temperatures observed during the whole sampling period, respectively. We used a similar notation for relative humidity, RH. Based on the information above regarding beetle response to climate, the merged climate factor c was defined as the relationship:$$begin{aligned} c = frac{T}{RH};;, end{aligned}$$
(3)
for (RHne 0). That is, higher temperatures and/or drier conditions (expected to favor D. gazella) lead to higher values for c. On the other hand, lower temperatures and/or more humid conditions (expected to favor dwellers) imply lower values for c. Intermediate values of c can represent either moderate or extreme values for both T and RH.Identifying ecological states and quantifying resilienceWith our (f_{dwell}) data as an index of community composition (i.e. system state), we calculated kernel density functions to interpolate a continuous probability distribution of the relative fraction of dwellers in the community, (p_{n}(f_{dwell})) (function density, R software 3.1.647) for a given range of climatic driver c values. We grouped the (f_{dwell}) data using ranges for c of size 0.4, to ensure a significant amount of weekly samples that allowed for the reconstruction of these probability distributions (see Table S1, first column). Note that bins with extreme values showed few data points (see first and last rows in Table S1), and thus were rejected to prevent misleading results due to reduced sampling. Also note that, for the density function, we used the default Gaussian kernel with a smoothing bandwidth adjusted to be (50%) larger than the default value (“adjust” argument set to 1.5). This conservative choice aims to reduce the effect of the different sampling across c bins and to ensure that differences among distributions across c values are not the result of spurious sampling noise.Further, we transformed the kernel density function:$$begin{aligned} V(f_{dwell}) = -ln (p_{n}(f_{dwell})) end{aligned}$$
(4)
This (V(f_{dwell})) function, called potential (e.g.48), shows by design well-defined minima for the most frequently observed values of (f_{dwell}) (i.e. configurations most frequently observed for the community, which conform the modes of the probability distribution) in a given group of data. At these points, the potential exhibits a change of trend from decreasing to increasing, and therefore its derivative shows a change of sign. Eq. (4), thus, provides a simple criterion to identify possible system states, which is a reason why potentials have been used extensively across disciplines49,50,51. Nonetheless, because the position of extrema is invariant under the transformation, using probability distributions instead would not alter our conclusions.Representing the potential obtained from all the (f_{dwell}) system states associated with a same range of climatic driver c values allowed us to identify stable community configurations associated with a specific climate. The comparison of the potentials obtained for different c ranges enabled the description of how the community changed in response to climatic variation. The location of the minima revealed which states were stable for a given value of the climatic driver; the presence of two minima, then, flagged the existence of bistability (i.e. two different community compositions possible for the same c value).These minima are materialized as wells in the potential’s landscape, which provides an easy way to understand the concept of stability: the dynamics of the system for the given value of the driver will eventually “fall” into a well (either a state dominated by dwellers or a state dominated by tunnelers), with the shape of the well (e.g. its depth) determining how difficult it is for the system to “escape” that state. Therefore, the area inside a well provides quantification of the tendency of a system to stay in that specific state, i.e. the resilience of the associated ecological state or how strong a perturbation has to be to move the system from such an ecological state to another2,3,50,51,52,53. Thus, in addition to number and location of wells, measuring their associated area allowed us to further characterize the resilience of the community. To this end, we first set a visualization window common to all potentials. Specifically, we plotted the potentials within a range for the vertical variable (the potential, V) given by ([-1.5,1.5]); the horizontal variable (fraction of dwellers, (f_{dwell})) is by definition bounded between 0 and 1. For potentials that showed one single well, the area of the well was measured as the area above the potential curve within this visualization window. For potentials that showed two wells (bistability), we measured the value of the potential at the local maximum separating the two wells, and established that value as the upper (horizontal) line closing the area of each well. To ensure all cases were comparable and eliminate any arbitrariness of the choices above, we expressed resilience as a relative area; in other words, we further normalized the well area by the total area across wells for that potential, which means that any single-well case will show a resilience (or relative area) of 1, and the resilience of the two wells when there is bistability adds up to 1.Figure 1Left: Community composition by functional group for all weeks of observation1. Green represents dwellers, blue represents tunnelers, and orange represents the invader D. gazella. Right: Sketch of responses of the community composition to the climatic driver (i.e. phase diagram) expected from the physiological and behavioral characteristics of the functional groups in the community as described in text: linear (red), or non-linear but monotonic without (blue) or with (brown) hysteresis.Full size imageIdentifying ecological transitionsMeasuring a state variable, (f_{dwell}), and a driver, c (order and control parameter, respectively, in the jargon of regime shift theory), allowed us to study how their observed behavior over time materializes in a driver-state relationship (the so-called phase diagram) defining the possible shifts in dominance (i.e. regime shifts) that the community may undergo as climate changes12. The non-monotonic temporal behavior of the components of the order parameter (i.e. dwellers and tunneler availability) and the components of the control parameter (i.e. temperature and relative humidity) makes it difficult to predict the shape of the phase diagram, and therefore whether we can expect alternative stable states in the focal example. For such cases, the dominance of the dung beetle community could (1) shift in a linear fashion toward the functional group favored by climatic conditions; (2) shift between functional groups in non-linear threshold response to climatic conditions without hysteresis; or (3) shift between functional groups in non-linear threshold response to climatic conditions with hysteresis –and thus showing bistability (see Fig. 1b, or12). Other possibilities, e.g. a non-linear shift between functional groups where one group is favored at intermediate climatic conditions12 are discarded as the invader is better suited for warmer and drier conditions. To evaluate which of these possibilities occurred, we represented (f_{dwell}) as a function of c, as well as the location of the minima shown by the potentials above. In addition to the emerging shape of this relationship, this plot can reveal the presence of alternative stable states if two or more different points occur for the same value of the control parameter, c. More
