Seed dormancy explains plant response to mass mortality events

Mass mortality events (MMEs) are large-scale, rapid die-offs resulting in extreme inputs of carrion biomass. Recent work demonstrates the effects of increasing carrion biomass on plant communities modulated by vertebrate scavengers and herbivores. However, the mechanisms underlying plant response to MMEs remain unclear. We hypothesized that carrion decomposition would interact with vertebrate herbivory and scavenging to generate distinct ecological filters on plants grouped by three seed dormancy classes (no dormancy, physiological dormancy, and physical dormancy). We designed a replicated field experiment crossing two levels of carrion biomass (~30 kg and ~360 kg) with three levels of vertebrate exclusion (open/no exclusion, scavenger exclusion, and herbivore exclusion) to quantify plant extirpation and colonization, plant performance, and the response of simulated seed bank and seed rain. We measured carrion decomposition rate, plant tissue nutrients (N, P, and K), seed survival, plant height, flower production, and plant community changes over 3 years. Carrion biomass levels associated with mass mortality increased plant tissue nutrients and plant extirpation and decreased seed bank survival, likely promoting plant colonization from seed rain. Vertebrate exclusion determined colonization probability in different ways, depending on seed dormancy classes. Vertebrate scavenger exclusion delayed decomposition, and the resulting environment favored the colonization of plants with impermeable seed coats (physical dormancy). In contrast, vertebrate herbivore exclusion promoted rapid germinators (no dormancy) that arrived as seed rain and quickly capitalized on the nutrient-rich exposed soil. With both functional roles intact (open/no exclusion), vertebrate scavenging ameliorated the negative effects of decomposition on seeds and plants, and selective herbivory on nutrient-rich plant tissue reduced overall plant height and flower production. The activity of both functional roles favored plants of the dominant physiological seed dormancy class and those with physical dormancy. Our results suggest that seed response to carrion decomposition and vertebrate functional activity can inform plant community response to MMEs.

Methods

Materials and Methods

Replication statement

We conducted simulated mass mortality events of herbivores and scavengers by crossing exclusion and carrion biomass treatments. Several response variables related to vertebrate and vegetation response enabled measurement of these treatments. To facilitate the interpretation of replication and scale, we present these response variables in the context of scale of inference, scale of treatment, and number of appropriate replicates.

Response variable

Scale of inference
Scale at which factor of interest is applied
Number of replicates at the appropriate scale

Response of vertebrate herbivores (# of observations)
Plot (70 m2)
Plot (70 m2)
4, 4 blocks with a 3x2 factorial design (n = 6) + 1 reference plot (6 treatment plots, 1 reference plot); block is the level of replication

Response of vertebrate scavengers (# of observations)

Carrion decomposition rate

Plant colonization (binary 0/1)

Plant extirpation (binary 0/1)

Plant tissue nutrients (effect size)

Seed survival (binary 0/1, effect size calculated by comparing to reference plot)
Microsite
Plot (70 m2)
4, 4 blocks with a 2x2 factorial design (n =4) nested within 3x2 factorial design (n = 24) + 2 packets per reference plot (6 treatment plots, 1 reference plot); block is the level of replication

Plant reproduction (# of flowers)
Plant
Plot
13, 26 total plants (opportunistic design does not account for location of plants in different clusters/plots)

Overview and Design

Location

We conducted this experiment on Oswalt Ranch, a managed rangeland (i.e., a mix of mostly native grasses, shrubs, and trees along riparian corridors) owned by Noble Research Institute (NRI) that is underlain by poorly drained and shallow soils (34°11'8.75" N, 97°5'12.42" W). Oswalt Ranch is in a humid subtropical region of southcentral Oklahoma with average temperatures of ~17 C° and average rainfall of ~870 mm (1981–2010; Oklahoma Climatological Survey, 2022). Plants on site with no seed dormancy are often annual species; plants with physical dormancy include legumes (Fabaceae), bindweeds (Convolvulaceae), and animal-dispersed seeds that require gut scarification; and most of the remaining plant species, including native and improved pasture grasses (Poaceae), are often physiologically dormant (Supplemental materials 1). White-tailed deer (Odocoileus virginianus) are common herbivores. Black vultures (Coragyps atratus), turkey vultures (Cathartes aura), and coyotes (Canis latrans) are dominant scavengers. Non-native and domesticated vertebrates include wild pigs (Sus scrofa) and cattle (Bos taurus). Noble Research Institute staff culled wild pigs to protect native ecosystems and property (Gaskamp et al., 2021). After a scheduled culling event, they deployed the resultant carcasses that we monitored for our experiment. Given that this culling event occurred distinct from our experiment, we did not require IACUC approval. In addition, we required no permit to conduct research on this private land.

Experimental design

NRI staff constructed four blocks (0.5 ha) spaced at least 750 m apart. Each block consisted of six plots randomly assigned a combined treatment based on two factors: biomass (single or mass mortality) and exclusion (open, scavenger, and herbivore; Figure 1A). Outside of each block, we also constructed a fenced reference plot (i.e., no carrion and no herbivory) for baseline plant tissue nutrient and seed survival measurements. We did not use netting to restrict avian access to reference plots. MME simulation plots contained ~12 times the biomass (~360 kg) as plots simulating background mortality (~30 kg). This biomass density does not approach that of some natural MMEs and experiments (e.g., 36 kg per m2; Steyaert et al., 2018; Baruzzi et al., 2022a) and thus likely represents a conservative biomass density. However, we had almost 1,200 kg of biomass in each block with this density of biomass (~7 kg per m2, ~50 carcasses per block), which is difficult to accomplish with simulated MMEs (Lashley et al., 2018; Fey et al., 2019).

We simulated MMEs of vertebrate herbivores and scavengers using carcasses with exclusion fencing and netting deployed in two phases, where each phase aligned with the period when the functional role exerts influence on the ecosystem. Beginning with carrion deployment (2019-04-15), we installed cages topped with netting around scavenger exclusion plots, leaving all other plots open. After active decomposition (2019-05-05), we removed the cages from the scavenger exclusion plots and established fencing in the herbivore exclusion plots to exclude vertebrate herbivores from the plots during plant recovery. The cages remained on herbivore exclusion plots for more than three years after carrion deployment until the conclusion of the sampling period (2021-07-18). Because this site actively grazes cattle in pastures, NRI employees also installed a two-strand polywire electric fence to exclude cattle from the entire block while still allowing vertebrate wildlife to enter the block (Figure 1A).

Software and analysis

All analyses conducted in R studio version 4.2 (R Core Team, 2023). We cleaned and visualized using the tidyverse package (Wickham et al., 2019). Linear model structures are guided by experimental design and hypotheses rather than model selection. In some cases, complexity was reduced for model convergence. We visualized the distributions of data prior to selecting a distribution and assessed model assumptions using statistical tests and visualizations of normality, outliers, independence, and homoscedasticity using the DHARMa package (Hartig, 2021). Unless otherwise noted, we present estimated marginal means and effect sizes from emmeans, which corrects for multiple comparisons using the Tukey method (Lenth, 2022). For models with interaction terms, we conducted type 3 analyses using the Anova function from the car package (Fox & Weisberg, 2019)

Measuring Ecosystem Response to Carrion 

1.1 Measuring response of vertebrate functional roles
We measured vertebrate activity in the plots to (1) determine the effect of carrion biomass on herbivore and scavenger activity and (2) ensure the efficacy of the exclusion efforts on functional role suppression. At each plot, we mounted a camera (Bushnell Trophy Cam HD; Bushnell Outdoor Products, Overland Park, KS) ~1 m from the plot edge facing the same cardinal direction (Figure 1B) with the following settings: high sensitivity, 5-minute delay, 15-second video. We deployed cameras 3 days prior to carrion deployment and retrieved them after the first year of the experiment. Our measure of activity results from counting all individual vertebrates in each video (i.e., the same individual animal in consecutive videos would count as two observations.

1.2 Response of vertebrate functional roles analysis

After processing the camera trap data, we grouped vertebrate observations (counts) into scavengers (obligate and facultative), herbivores, and other animals. Then, we analyzed total scavenger and herbivore observations during active decomposition (2019-04-15–2019-05-05) and plant recovery (2019-05-05–2019-10-06). We used the MASS package to conduct a generalized linear model with a negative binomial distribution for total scavenger observations during the active decomposition period (Venables & Ripley, 2002) and base R to conduct Poisson regression for total herbivore observations during the plant recovery period (R Core Team, 2023). These models would not converge with any degree of complexity. Both models thus include only treatment (carrion biomass and vertebrate exclusion combined) and site as predictor variables.

2.1 Measuring carrion decomposition rate

We photographed carcasses each day in the week following deployment, for another three times in 3-day intervals, and once more 7 weeks later (n = 11). In MME treatments, we randomly selected three carcasses during each sampling event. We analyzed pictures using the Total Body Score method (sensu Keough et al., 2017).

2.2 Carrion decomposition rate analysis

We used the glmer function from the lme4 package to conduct two general mixed effects models evaluating the number of days elapsed before a carcass entered advanced and skeletal decomposition stage (Bates et al., 2015). We used biomass and functional impairment as fixed effects and block as a random effect (i.e., intercept).

3.1 Nested seed survival experiment

We crossed timing of seed arrival (seed bank or seed rain) with location (carrion contact or 0.5 m from carrion) by placing a single sealed pouch in each position (10 x 10 cm, 1 mm mesh screening). We placed pouches simulating seed rain on top of carrion or on top of the soil and simulated seed bank pouches beneath the soil (Figure 1C). In both cases, we affixed a metal wire to the pouch and a pole at the center of the plot to prevent removal and enable rediscovery. Having no carrion, the reference plot did not need the distance comparison used in single carrion and MME plots and received only two pouches (seed rain and seed bank). Each block received a unique set of random species belonging to one of three functional roles (Table 1), and the number of seeds inside the pouch varied by species due to seed size and availability (Figure 1D). 

Table 1. Simplified classes reflect the commonly encountered types of seed dormancy.

Physical dormancy

Physiological dormancy
No dormancy

Malva sylvestris 
Triticum aestivum
Sinapsis alba

Rhus glabra

Syringa spp.
Beta vulgaris

Bixa orella
Hordeum vulgare
Pisum sativum

Glycine max
Capsicum annum
Pinus sylvestris

 We deployed 2,258 seeds from 12 distinct species 1 day after carrion deployment. Then, we collected the pouches 21 days after deployment and recorded the number of seeds recovered, perished, and survived (including living germinated seeds). We recovered 84% of the seeds. After this initial assessment, we attempted to germinate any remaining uncategorized seeds (70% of recovered seeds) to evaluate viability (e.g., perished or surviving). Seeds with no dormancy germinated with moisture exposure; seeds with physiological dormancy required stratification before moisture exposure; and seeds with physical dormancy required scarification and stratification before moisture exposure (Baskin & Baskin, 2004). We repeated this process three times on the remaining potentially dormant seeds. A small subset of seeds remained undetermined (~1%).

3.1 Nested Seed Survival Experiment analysis

We constructed a binomial generalized linear model with a binomial distribution using the number of surviving or perished seeds out of total seeds in each packet as a response variable (i.e., the number of successes and failures out of total trials) with the glmer function from the lme4 package (Bates et al., 2015). The predictor variables for the model included carrion biomass, exclusion treatment, timing (seed bank or seed rain), location (carrion contact or 1 m from carrion), and seed dormancy class, with block as a random effect. Our observations (n = 276) did not permit a full interactive model, so we included the interactions that best represented our hypotheses regarding the differential responses of seed dormancy classes to MMEs and the response of seed banks and seed rain to carrion decomposition.

Seed survival ~ Carrion biomass * Exclusion treatment * Dormancy class + Dormancy class * Timing * Carrion biomass + Location + (1| Block)

We did not include reference plots in this model but separately estimated the MLE of survival probability as a baseline to visualize the results. Then, we calculated the difference in survival compared to the baseline at each block and summarized the mean odds ratio for each treatment. Then, we further summarized these means across groupings of treatments (i.e., mean for all seed bank treatments):

4.1 Measuring extirpation and colonization of seed dormancy classes

We sampled vegetation before carrion deployment (2019-03-17), each month during the initial growing season (2019-04-21, 2019-05-05, 2019-06-02, 2019-07-21), and once in each of the following 2 years (2020-07-20, 2021-07-18) using the point-intercept method at 1 m intervals for 4 m in each cardinal direction. At each point, we marked the presence of species (i.e., multiple species may occur at each point but not multiple individuals belonging to a single species). We also recorded bare ground, litter, dead plants, and carrion. For each plot and sampling date, we calculated the cover value in the resulting community matrix for each species by dividing the number of observed points by the total possible points (16 total points per plot). Then, we classified the surveyed plants according to seed dormancy using Baskin & Baskin’s (1998) seed dormancy database and published literature (Supplemental materials 1).

4.1 Extirpation and colonization of seed dormancy classes analysis

We converted the plant community time series into two datasets describing the probability that individual plants belonging to seed dormancy classes (1) vanished from transect points after the initial surveys or (2) colonized transect points after the experiment. Colonization or extirpation thus refers to an individual point, and we make no effort to confirm absence or presence across the transect or plot. This approach also makes a key assumption that whether we detected a species indicates the actual presence or absence of a species at the transect point. For the extirpation model, we subsetted the data into transect points with species belonging to physiological dormancy in the survey before carrion decomposition. Then, we determined whether absences of physiological dormancy occurred at that point to represent extirpation. We only used physiological dormancy for this model because nearly all individuals present during the initial pre-treatment survey (2019-03-17) belonged to this dormancy class. For the colonization model, we tracked each transect point and recorded any time a dormancy class absent during the initial pre-treatment survey colonized by the final sampling period (2021-07-18), again resulting in a binary response variable for colonization. We used the glm function to conduct binomial generalized linear models analyzing extirpation with treatment as a predictor and colonization with dormancy class and treatment as interactive predictors (R Core Team, 2023). We did not include further complexity, including block as a predictor, due to convergence issues.

5.1 Plant tissue sampling

We randomly collected aboveground tissues from similar plant morphotypes in each block to meet biomass requirements. The Mississippi State University Soil Laboratory (Starkville, MS, USA) conducted plant tissue analysis. We focused on three essential plant macronutrients that (A) leach from carcasses and (B) may affect plant selection in white-tailed deer: nitrogen, phosphorous, and potassium (Barton et al., 2013a; Dykes et al., 2018).

5.2 Plant tissue analysis

We analyzed plant tissue nutrients using the adonis function from the vegan package to conduct a permutation-based multivariate ANOVA (Oksanen et al., 2022). The model included treatment and date as interactive predictor variables, block as a random effect, and Hellinger-transformed (e.g., standardized) macronutrient values for the response variable. To visualize the effects of biomass and exclusion on nutrients, we calculated the log ratio of treatments compared to reference plots as an effect size:

6.1 Opportunistic Plant Fitness Sampling

We developed an opportunistic sampling regime (stratified random selection) upon visually observing dramatic increases in plant fitness correlates in simulated herbivore MMEs (Younginger et al., 2017). During the final year of the experiment (2021-07-18), we measured the height and inflorescences of a dominant annual or biennial plant (Lactuca serriola) on individuals evenly distributed within and outside herbivore MME plots (n = 26). To analyze the opportunistic measurement of plant performance metrics, we conducted separate Welch’s t-tests using base R to determine the effects of treatment on height and inflorescences (R Core Team, 2023).