Resource allocation to growth, reproduction, and body maintenance varies within species along latitudinal gradients. Two hypotheses explaining this variation are local adaptation and counter-gradient variation. The local adaptation hypothesis proposes that populations are adapted to local environmental conditions and are therefore less adapted to environmental conditions at other locations. The counter-gradient variation hypothesis proposes that one population out performs others across an environmental gradient because its source location has greater selective pressure than other locations. Our study had two goals. First, we tested the local adaptation and counter-gradient variation hypotheses by measuring effects of environmental temperature on phenotypic expression of reproductive traits in the burying beetle, Nicrophorus orbicollis Say, from three populations along a latitudinal gradient in a common garden experimental design. Second, we compared patterns of variation to evaluate whether traits co-vary or whether local adaptation of traits preclude adaptive responses by others. Across a latitudinal range, N. orbicollis exhibits variation in initiating reproduction and brood sizes. Consistent with local adaptation, (1) beetles were less likely to initiate breeding at extreme temperatures, especially when that temperature represents their source range; (2) once beetles initiate reproduction, source populations produce relatively larger broods at temperatures consistent with their local environment. Consistent with counter-gradient variation, lower latitude populations were more successful at producing offspring at lower temperatures. We found no evidence for adaptive variation in other adult or offspring performance traits. This suite of traits does not appear to coevolve along the latitudinal gradient. Rather, response to selection to breed within a narrow temperature range may preclude selection on other traits. Our study highlights that N. orbicollis uses temperature as an environmental cue to determine whether to initiate reproduction, providing insight into how behavior is modified to avoid costly reproductive attempts. Furthermore, our results suggest a temperature constraint that shapes reproductive behavior.

Nicrophorus orbicollis populations

We derived the laboratory beetle populations used for this experiment from wild-caught beetles captured with baited pitfall traps near Big Falls, Wisconsin (high latitude; HL; 44.6165°N, -89.0161°W), Waveland, Indiana (mid latitude; ML; 39.9417°N, -87.0917°W), and Spavinaw, Oklahoma (low latitude; LL; 35.3704°N, -95.0486°W) (Fig. 2a) in May-June of 2014 and 2015. We housed all N. orbicollis in individually marked plastic containers (15 x 11 x 7 cm) in an environmental chamber at 21^oC with a 14:10 h light:dark (L:D) cycle and fed chicken liver ad libitum. These conditions simulated the natural light/dark pattern and temperature consistent with the beetles’ summer breeding season in their natural environment (Cook et al. 2019). We used these wild-caught beetles to establish the first generation (F1) laboratory populations used for experiments. We bred wild-caught beetles by placing a male and female with a fresh mouse carcass in plastic containers (18 x 15 x 10 cm) filled two-thirds full with top soil. We removed the wild-caught males when larvae first appeared on the carcass, and we removed wild-caught females when larvae dispersed from the carcass. We left F1 larvae undisturbed until eclosion (approximately 28-30 days), and then we maintained them in individual plastic containers as described above until used for the experiments at sexual maturity (21-28 days old post eclosion).

Common garden experiment

To assess genetically based latitudinal patterns, we tested representatives from each of the three source locations at each of five constant air temperatures (12-13°C, 15°C, 20°C, 25°C, or 27-28°C) in a common garden experimental design. We began each trial by randomly selecting a laboratory-reared, sexually mature (21-28 days old) male and female from the same population, and from different parental lines. At the beginning of each reproductive bout, we weighed females and males and measured their pronotum width. We placed each pair of beetles in a plastic container (18 x 15 x 10 cm) filled two-thirds full with commercially purchased topsoil and given a freshly thawed 30 g (± 1 g) mouse carcass (Fig. 3). Because either parent can successfully raise a brood if their partner is handicapped or removed (Smiseth, et al. 2005, Creighton et al. 2015), we removed males after 48 hours to allow sufficient time for mating to occur, but to minimize male impact on female life history characteristics thereafter (Rauter and Moore 2004, Smiseth et al. 2005, Creighton et al. 2015, Smith et al. 2015). We randomly assigned pairs from each of the three populations to a treatment of one of five temperatures: 12-13 ºC (T_min), 15°C (average three month daily low temperature for Wisconsin), 20 ºC (average three month daily temperature for Indiana), 25ºC (average three month daily high temperature for Oklahoma), or 27-28 ºC (T_max); (L:D;14:10; Fig. 3). Temperatures represent the average daily high and low temperatures experienced within the range of the collection sites for the known breeding months of May-August. We determined experimental breeding temperatures by calculating a three-month daily mean air temperature for each location using ten-year temperature data sets obtained from the National Climatic Data Center (NCDC) Annual Climatological Summaries of U.S. station data (ncdc.noaa.gov accessed on 9/10/2013). We checked, photographed, and monitored broods daily to measure response variables. We monitored all broods until beetles completed the reproductive cycle or until we determined brood failure (i.e., carcass preparation stopped). Each temperature treatment by location combination had at least 10 replicates, but more broods were required at the intermediate temperatures to allow for analysis of offspring traits. We summarized sample sizes for each temperature treatment by location combination in Table 1.

Parental reproductive performance

We measured both parental and offspring response variables from each brood to characterize variation in reproductive performance across the temperature gradient. To measure parental performance, we evaluated the following variables: reproductive success, hatching asynchrony, offspring number, and developmental timelines. Each of these response variables yielded one measure per brood, so sample sizes are equivalent to those in Table 1 for the three intermediate temperatures.

Carcass preparation and reproductive success

We measured reproductive success by two response variables – degree of carcass preparation and probability of producing offspring that survived to adulthood. We scored carcasses for degree of preparation and assigned a score of zero to four to represent the stage of carcass preparation (Table 2). We evaluated the final stage of carcass preparation for each brood as a measure of reproductive success at different temperatures. To characterize the probability of producing offspring at a given temperature, we scored each brood as a success (score = 1) if any adult offspring were produced, or a failure (score = 0) if no adult offspring were produced.

Hatching asynchrony

We calculated brood hatching asynchrony as the spread of hatching dates within a brood (i.e., the number of days between the first and last arriving first instar larvae on the prepared carcass; Aparicio 1999). While there may be an adaptive value to increased asynchrony, for this analysis we assumed that less asynchrony in larval arrival to a carcass is advantageous.

Parental developmental timelines and reproductive output

To determine parental developmental timelines we used two variables: time in days, 1) to fully prepare the carcass, and 2) for larvae to fully develop and consume the carcass (Scott 1998). To characterize reproductive output, or number of offspring produced, we determined the number of final (3^rd) instar larvae that dispersed from the carcass.

Offspring performance

To measure offspring performance, we evaluated the following variables: offspring growth rate, adult offspring body size, developmental stability (as measured by degree of fluctuating asymmetry of newly eclosed adults), and percentage offspring body fat.

Offspring growth rate and brood size

To evaluate offspring growth rate, we subtracted the average mass of an individual larvae from the first day they were present on the carcass from the average mass of an individual larvae on the final day they were present on the carcass. We then divided this value by the total number of days that larvae were on the carcass to account for asynchronous arrival and dispersal to the carcass by larvae. The response variable was mean offspring growth rate per brood. To measure offspring body size, we measured the adult offspring body mass from all offspring in the brood at the time of eclosion. We used mean offspring body size as the response variable for analysis.

Offspring developmental stability

To measure developmental stability and percentage body fat, we randomly selected one adult male and one adult female offspring from each brood. This selection resulted in two replicates of developmental stability and percent body fat per brood for analysis. These individuals were pinned for the purpose of taking photographs which were then used to measure fluctuating asymmetry. We used three variables: 1) an anterior to posterior transect through the beetles upper orange elytra spot, 2) a basal to distal transect from the elytra edge to the lower orange elytra spot, and 3) a basal to distal transect from the “pronotum cleft” to the edge of the pronotum (Fig. 2b). We used the difference between the left and right side for each of these variables as a measurement of developmental stability (Van Valen 1962). Additionally, we determined, percentage body fat on these same individuals following the lipid extraction techniques used by Marden (1989).

Data analysis

No beetles produced offspring at either the highest or lowest of the five temperatures in the experimental design. For this reason, only degree of carcass preparation was analyzed at all five temperatures. Because the response variable was discrete (i.e., scored as an integer from 0-4) we used a generalized linear model with a log link function and a Poisson distribution (Neter et al. 2004). The model had two fixed factors: temperature treatment (5 levels) and location (3 levels). We included the interaction between temperature treatment and location, and female and male body size as covariates. We used Proc GENMOD in SAS for the analysis (SAS 9.3 SAS Institute, Cary, North Carolina, USA).

We modeled the probability of producing offspring (i.e., reproductive success) as a binomial response (0 or 1) with a probit link function in a generalized linear model framework with a Poisson distribution (Neter et al. 2004). Because no beetles produced offspring at the highest and lowest temperatures in the experimental design, we used only the three intermediate temperatures for this analysis. The model had two fixed factors: temperature treatment (3 levels) and location (3 levels). We included the interaction between temperature treatment and location, and female and male body size as covariates. We used Proc GENMOD in SAS for the analysis (SAS 9.3 SAS Institute, Cary, North Carolina, USA).

Number of offspring produced (i.e., final brood size), hatching asynchrony, and developmental timing (2 response variables) were represented as count data, so we used a generalized linear model with a log link function and a Poisson distribution (Neter et al. 2004). The model had two fixed factors: temperature treatment (3 levels) and location (3 levels), and we included the interaction between temperature treatment and location. For number of offspring produced, time preparing the carcass, and time until dispersal of offspring, we included female and male body size as covariates. For the model of hatching asynchrony, we included female mass and brood size as covariates. We used Proc GENMOD in SAS for each of these analyses (SAS 9.3 SAS Institute, Cary, North Carolina, USA).

Growth rate of offspring and size of adult offspring (measured as mass in g) were continuous response variables so we used a general linear model and a Poisson distribution (Neter et al. 2004) on untransformed data for analysis. Data consisted of means calculated for each brood so we have one replicate per brood. Raw data exhibited normally distributed residuals and fit the assumptions of the model. The model had two fixed factors: temperature treatment (3 levels) and location (3 levels), and we included the interaction between temperature treatment and location. For the growth rate model we included female combined parental body size and brood size as covariates, and for the offspring size model we included female body size, male body size, and brood size as covariates. We used Proc MIXED in SAS for the analysis (SAS 9.3 SAS Institute, Cary, North Carolina, USA).

To evaluate percent body fat and developmental stability (three measures of fluctuating asymmetry; S1) of adult offspring we used a general linear model and a Poisson distribution (Neter et al. 2004). The models had three fixed factors: temperature treatment (3 levels), location (3 levels), and sex (2 levels), and we included all interactions among temperature treatment, location, and sex. For the model evaluating percent body fat, we included two covariates – offspring size (pronotum width) and brood size. For the models evaluating fluctuating asymmetry, we included brood size as a covariate. Raw data exhibited normally distributed residuals and fit the assumptions of the model. We included the ID number of the brood as a random effect in these models because we measured two individuals from each brood. We used Proc MIXED in SAS for the analysis (SAS 9.3 SAS Institute, Cary, North Carolina, USA).