Seasonal patterns of resource use within natural populations of burying beetles

For organisms in temperate environments, seasonal variation in resource availability and weather conditions exert fluctuating selection pressures on survival and fitness, resulting in diverse adaptive responses. By manipulating resource availability on a local spatial scale, we studied seasonal patterns of resource use within natural populations of burying beetles (Nicrophorus vespilloides) in a Norfolk woodland. Burying beetles are necrophagous insects that breed on vertebrate carcasses. They are active in Europe between April and October, after which they burrow into the soil and overwinter. Using breeding and chemical analyses, we compared the fecundity and physiological state of beetles that differed in their seasonal resource use. We found seasonal variation in carrion use by wild burying beetles and correlated differences in their reproductive success and cuticular hydrocarbon profiles. Our results provide novel insights into the seasonal correlates of behaviour, physiology, and life history in burying beetles.

Methods
PART A: EXPERIMENTAL PROTOCOLS AND METHODS

Question 1 or 3: Is there seasonal variation in the trapping frequency of burying beetles on chick and mice carrion between June and August?

Study area and trapping methods

We sampled the burying beetle population at Thetford Forest (52°20'39.5" N 0°32'14.9" W), Norfolk, UK from May to October 2017 at 10 different trap locations (See Methods table 1 below), under permit from Forestry Commission England. We used carrion-baited beetle traps (Japanese Beetle Trap Kit from Scotts Co., not treated with any pheromones), suspended in vegetation 1-2m above ground. The bottom half of the trap was filled with Miracle-Gro compost, and a small dead vertebrate was placed on the top as bait. The contents of the trap were collected at intervals and brought back to the lab for processing. The trap was then refilled and rebaited. After processing in the lab, no beetles were released back into the field.

Beetles were sampled using a paired-trap arrangement, in which we placed two beetle traps- one baited with a dead domestic chick and the other baited with a dead mouse- near each other at each trap location and recorded the number of beetles found in each trap. The traps within each experimental pair were placed 1-2 m apart. Pairs of traps were placed 200- 400 m apart from each other. With this design, beetles were given a simultaneous choice between a dead mouse and a dead chick. Each time we rebaited a trap with carrion, we rebaited it with the alternate carrion type. Therefore, if a mouse carcass had been placed in the trap previously, it was replaced by a chick carcass on the next sampling trip to ensure that the trap location itself did not bias beetle catch. The mice and chick carcasses used were matched in weight (30-40 g).

Processing field-caught beetles

At the lab, we used carbon dioxide to immobilise each beetle and brush off any mites stuck to it. We recorded the pronotum width and sex of each N. vespilloides beetle we trapped. We compared beetles collected at two different time points during the burying beetle season: the first set was collected in June 2017 after 10 days of trapping between 4 June and 14 June, and the second set was collected in August 2017 after 15 days of trapping between 4 August and 19 August. The 10 trapping locations (Methods Table 1) were the same across both sampling periods.

Question 2 or 3: Does reproductive success vary with carrion substrate and/or season?

Measuring reproductive performance

After collecting beetles from the traps, measuring and identifying them, we put each N. vespilloides individual into its own personal small plastic box (12 cm × 8 cm × 2 cm) and fed it 1 g of beef mince. The beetles were stored alone in their boxes for 7-10 days to ensure that any newly eclosed individuals had had sufficient time to become sexually mature before we measured their reproductive performance.

For breeding, we placed a pair of beetles (one male and one female) in a larger plastic breeding (17 cm x 12 cm x 6 cm) box half-filled with Miracle-Gro compost and provided with either a chick or mouse carcass. Each member of the pair had been trapped on the same type of carrion and we bred them on the same carrion they were trapped upon.  This method was used twice, once for beetles collected in June and once for those collected in August, yielding four treatments in all.

The mass of the carcass provided for reproduction was recorded and kept consistent within each treatment. We then placed the breeding box inside a cupboard so that it was shielded from light to mimic the low light conditions typically experienced by beetles as they breed below ground. Eight days after pairing the beetles (i.e., the point at which the larvae had completed development and were starting to disperse away from the remains of the carcass), we counted and weighed the surviving larvae from each pair.

We used the following measures to record reproductive success in our experiments:

Brood failure: We recorded the total number of broods that failed to produce any larvae: 0 denoted broods that failed, 1 denoted those that had at least one surviving larva at 8 days post dispersal.
Brood size:  The total number of dispersing larvae 8 days post-breeding.
Average larval mass: Total mass of the brood at dispersal (g) divided by the brood size.
Larval density: Brood size divided by the mass of the carrion used for breeding (g).
Carcass use efficiency: (total brood mass (g) divided by original carrion mass (g)) x 100%

In June, 53 pairs of beetles trapped on mice (MM) and 24 pairs of beetles trapped on chicks (CC) successfully produced broods with at least one larva. There were 4 failed broods (3 on mice carcasses and 1 on a chick carcass). In August, 16 pairs of beetles trapped on mice (MM) and 25 pairs of beetles trapped on chicks (CC) produced broods with at least one larva. There were 7 failed broods (2 on mice carcasses and 5 on chick carcasses). The failed broods were excluded from analyses of reproductive success.

Question 3 or 3: Do beetles that are attracted to different types of carrion also differ predictably and seasonally in their CHCs? 

For this experiment, we sampled a total of 63 females; 32 were trapped on chicks and 31 were trapped on mice. 40 females were collected on 23 May 2017 (“early” season- 20 on chicks and 20 on mice). 6 females were collected on 14 June 2017 (“mid” season- 3 on chicks and 3 on mice). 17 females were collected on 4 September 2017 (“late season”- 9 on chicks and 8 on mice). After removing the mites from the body of the beetles, we isolated up to two female beetles from each trap individually in a glass vial for 15-20 min before storing them in a fresh vial at -80 °C. Later, we processed the beetles for CHC extraction by allowing them to thaw at room temperature for 30 min. We then soaked them in 4 ml of solvent (99% hexane, HPLC grade) for 20 mins. We transferred the extract obtained to a clean vial and allowed it to evaporate completely in a fume hood under nitrogen gas. At this stage, the sealed vials were shipped to Prof. Patrizia d'Ettorre’s lab at Université Sorbonne Paris Nord for analysis and characterisation.

CHC analysis and characterisation

We resuspended the extract in 400 µl of pentane (HPLC grade) and added an internal standard (C18, Octadecane at 16ng/µL) to each extract. The internal standard was used to determine the absolute amount of cuticular compounds present in each sample. We then analysed 2µl of the extracts using GC-MS (Agilent Technologies 7890A gas-chromatograph coupled to a 5975C Mass Spectrometer equipped with a HP5MS GC column (30 m × 0.25 mm × 0.25 μm) and operated at 70 eV in the electron impact ionization mode). The carrier gas used was helium at 1 ml/min. The column oven was programmed as follows: an initial hold of 1 min at 70°C, then increased to 200°C at 35°C/min, to 320°C at 4°C/min (held for 20 min).

We identified cuticular hydrocarbons based on their retention times (compared to standards) and fragmentation patterns. We manually integrated the chromatograms and converted the peak areas of the total hydrocarbon fraction using the MSD ChemStation software by Agilent Technologies, Inc.

PART B: DATA VISUALISATION AND STATISTICAL ANALYSES

Field and reproductive success data

We carried out all statistical analyses to test our predictions using R (RStudio version 1.3.959) with generalised linear models (GLM) and generalised linear mixed models (GLMM) using the lme4, glmmsr, and MASS packages. Analysis-of-variance tables for model objects were calculated using the ‘car’ package. Post-hoc comparisons using Tukey’s HSD test were carried out using the package ‘lsmeans’. The asymptotic test for the equality of coefficients of variation (CV) was carried out using the ‘cvequality’ package (Feltz & Miller 1996).

Question 1 or 3: Is there seasonal variation in the trapping frequency of burying beetles on chick and mice carrion between June and August?

We calculated the average number of beetles per day by dividing the total number of N. vespilloides beetles found in a trap by the number of days the traps had been left out. We focussed on the two different time points for which we also measured reproductive outcome, namely June and August 2017, using a GLMM that included carrion type and sampling month as fixed effects, and trap ID and sampling date (to account for any differences in sampling effort) as random factors with a Poisson error structure. The total number of N. vespilloides beetles found in a trap on the sampling day was used as the response variable.

Question 2 or 3: Does reproductive success vary with carrion substrate and/or season?

We examined the effect of month-trapped, carcass type used for breeding and their interaction on the following measures of reproductive success:

brood success versus failure, using a multivariate logistic regression model with a binomial error term 
the number of dispersing larvae, using a GLM with a Poisson error term 
average larval mass using a linear model 
larval density using a linear model
carcass uses efficiency using a linear model

When arriving at a minimal model using GLMs and GLMMs to explain our results, we removed non-significant terms and interactions using stepwise elimination. When presenting the results from post-hoc analyses, we list all the terms that were tested, and their statistics at the last point when they were retained in the model.

Question 3 or 3: Do beetles that are attracted to different types of carrion also differ predictably and seasonally in their CHCs? 

To analyse the chemical profile of both sets of beetles, we selected the 17 most regularly occurring GC-MS peaks (Methods Table 2). These represented the hydrocarbons we had identified and integrated using the MSD ChemStation software.

We carried out the principal component analysis, hierarchical clustering, and visualisation of the data using ggplot2, dplyr, pvclust, FactoMineR, and factoextra packages in R (RStudio version 1.3.959).

We log-normalised the peak areas within each sample using the following formula (Aitchison 1982):

Zij = ln [Yij/g(Yj)]

where Zij is the transformed area of peak i for beetle j; Yij is the area of peak i for beetle j; and g(Yj) is the geometric mean of the areas of all peaks for beetle j.

We used the standardised area values of the 17 peaks for hierarchical cluster analysis with Ward’s classification method to classify our samples. The significance of each node in the cluster was determined by multiscale bootstrap clustering with 10,000 iterations using the ‘pvclust’ package in R (Suzuki et al. 2019). We set the confidence level for the p-value threshold to 95% and ensured that only the most significant clusters (with a p-value lower than 0.05) were highlighted.

To examine how the CHC compounds found in our samples contribute to the discrimination of these samples, we performed a principal component analysis (PCA) on the log-normalised peak areas. We implemented the PCA using singular value decomposition (Hartman et al. 2023) for better numerical accuracy. Plotting standard deviations of principal components (PCs) and proportion of variances against all 17 PCs, we used the elbow method to determine the optimal number of PCs that explained the maximum amount of variance in our data (Jolliffe 2002).  Based on the visual inspection of the elbow plots, we retained the first 10 principal components,  which capture most of the variance in the data (96.3%). We then visualised the data using a heatmap depicting the loadings of the first 10 PCs. To check the correlation between the first two Principal Components and the original variables, we calculated the squared cosine value (cos2) for each variable by squaring the cosine of the angle between the vector with the variable's coordinates and the origin of the graph.

Cluster validation of our data indicated one outlier (Sample M13E). We confirmed this visually by using a 2-dimensional scatterplot before removing the outlier. We then repeated our PCA and clustering analysis without this data point.

PART C: APPENDIX

Methods Table 1: Geographical coordinates of trapping locations at Thetford Forest

Trap location number

Latitude

Longitude

1

52.3443120’N

0.5374620’W

2

52.3428090’N

0.5382080’W

3

52.3446849’N

0.5406787’W

4

52.3435720’N

0.5398280’W

5

52.3443100’N

0.5397050’W

6

52.3451601’N

0.5406294’W

7

52.3434220’N

0.5414500’W

8

52.3449940’N

0.5430960’W

9

52.3440520’N

0.5441710’W

10

52.3471700’N

0.5435390’W

Methods Table 2: Identification of the 17 most regularly occurring peaks in the cuticular hydrocarbon profile of N. vespilloides. Diagnostic ions are provided.

Retention time (min)

Compound

Diagnostic EI ions (m/z)

Internal standard

7.47

C18 (IS) (Octadecane)

254

1

10.59

C21 (Heneicosane)

296

2

11.95

C22 (Docosane)

310

3

13.46

C23 (Tricosane)

324

4

14.57

3MeC23 (3-methyltricosane)

57, 309, 281, 323

5

16.25

C25:1 (Pentacosene)

350

6

16.34

C25:1 (Pentacosene)

350

7

16.66

C25 (Pentacosane)

352

8

17.45

5MeC25 (5-methylpentacosane)

85, 309, 281, 351

9

17.85

3MeC25 (3-methylpentacosane)

57, 337, 309, 351

10

18.42

3,9diMeC25 (3,9dimethylpentacosane)

57, 155, 252, 351, 365

11

19.57

C27:1 (Heptacosene)

378

12

19.66

C27:1 (Heptacosene)

378

13

19.92

C27 (Heptacosane)

380

14

21.13

3MeC27 (3-methylheptacosane)

57, 365, 337, 379

15

21.69

3,9 diMeC27 (3,9-dimethylheptacosane)

57, 155, 281, 379, 393

16

22.81

C29:1 (Nonacosene)

406

17

23.17

C29 (Nonacosane)

408