Density-dependent variability in an eruptive bark beetle and its value in predicting outbreaks

Several species of aggressive bark beetle in the genus Dendroctonus (Coleoptera: Curculionidae) undergo large fluctuations in population density with distinct outbreak and non-outbreak phases. We investigated attributes we hypothesized as subject to density-dependent variation (in particular, those likely to express phenotypic plasticity) in the southern pine beetle (SPB), Dendroctonus frontalis, as possible indicators of population fluctuations. These traits were morphology (body size and hindwing shape) and the sex ratio of trap-captured, dispersing SPB populations. We compared attributes of beetles from locations that ranged from having zero to high numbers (>1500) of SPB infestations (spots) at the county level for two summers. Southern pine beetle were obtained from six states in the southeastern USA and had been collected during a springtime, region-wide trapping survey used for forecasting outbreaks annually. Although we detected an expected but weak sexual dimorphism in both size and shape-related traits, no morphological differences were found between SPB collected from traps in counties with low or high densities of spots (≤10 or >10 per county, respectively). We found no relationship between numbers of SPB spots per county and trapped sex ratios in 2016, but we observed a strong trend in 2017, with about three times higher proportions of females trapped in counties with low compared with high numbers of spots. This implies that one or more known or possible factors influencing trapped sex ratios (e.g., disparities between the sexes in their responses to semiochemicals or in their propensity to disperse) can be densitydependent. Including trap-captured sex ratios in prediction models may improve current forecasting of SPB outbreaks in the southeastern USA, informing more timely and effective management of one of the most economically and ecologically important beetle species of this region.


INTRODUCTION
Population densities can alter availability of resources such as food and breeding substrate and influence a range of other pressures, and adaptive density-dependent phenotypic plasticity of relevant traits can offer a selective advantage under either low or high densities (Simpson et al. 2005). Phenotypic plasticity can increase an organism's fitness through adaptive modifications of morphological, behavioral, and other characteristics to address the shifting demands of variable environments (Stearns 1989, Leclaire andBrandl 1994). Phenotypic plasticity is more likely to evolve in species with fluctuating populations (Svanbäck et al. 2009), and consequently, density-dependent phenotypic plasticity occurs in eruptive pest insect species (Anstey et al. 2009, Whitman andAgrawal 2009). Non-adaptive plasticity in traits may also occur in response to density (Stearns 1989).
The dispersal capacity of individual insects is substantially determined by their morphology and behavioral traits, whereas changes in population density can create pressures for dispersal (Travis et al. 1999). Plasticity in dispersal-associated traits can be advantageous, as investment of energy and resources in dispersal can cause a reciprocal loss in reproductive output (Simmons and Thomas 2004, Homburg et al. 2013, Fraser et al. 2014. Morphological traits associated with flight capacity, and therefore dispersal ability, include body size and mass, wing loading (the ratio between body mass or size and wing area), and wing shape (Angelo and Slansky 1984, Arribas et al. 2012, Hassall 2015, Kalberer and Kölliker 2017. Insect species with density-dependent polyphenism in traits associated with flight capacity have been identified in several orders of insects, including Diptera (BitnerMathé and Klaczko 1999), Hemiptera (Fujisaki 1986, Mori andNakasuji 1991), Hymenoptera (Kölliker-Ott et al. 2003), Lepidoptera (Angelo and Slansky 1984), and Orthoptera (Nishide and Tanaka 2013). Within Coleoptera, species that have been shown to exhibit differences in wing morphology associated with population phases include Callosobruchus maculatus Fabricius (Utida 1972), Hylobius abietis (L.) (Tan et al. 2010), Oreina cacaliae (Schrank) (Kalberer and Kölliker 2017), and Diabrotica virgifera virgifera LeConte (Mikac et al. 2019). Evidence for density-dependent plasticity has been reported in bark beetles (Coleoptera: Curculionidae: Scolytinae) (Bentz et al. 2011, Klutsch et al. 2020, and differing behavioral responses to tree defenses are known to occur in bark beetles from either outbreak or endemic populations (Lorio 1986, Wallin and Raffa 2004, Boone et al. 2011. Bark beetle behavioral responses to semiochemicals may exhibit plasticity depending on local population densities (Sullivan et al. 2011, Klutsch et al. 2020, and in the case of mountain pine beetle (Dendroctonus ponderosae Hopkins), individuals are less deterred by host pine defenses when at high densities (Wallin and Raffa 2004). However, the density dependence of flight-associated traits has received limited investigation in bark beetles (Sallé et al. 2005, Jones et al. 2019.
Certain species of bark beetles produce population eruptions that result in devastating, widespread mortality of their host trees (Turchin et al. 1991, Birt 2011, Jarvis and Kulakowski 2015. Colonization of healthy, high-quality hosts by aggressive bark beetles occurs predominantly in association with these eruptions, when beetles are present in sufficient numbers to mass attack such hosts and overcome their defenses Berryman 1983, Martinson et al. 2013). Rates of reproduction are high on these high-quality hosts, which can contribute to rapid population growth (Martinson et al. 2013). However, at low population densities, bark beetles are restricted to hosts with compromised defenses that are typically widely dispersed and low-quality (Clarke 2012, Raffa et al. 2016. Hence, alternative phenotypes (e.g., plasticity in wing morphology, body size, and responses to semiochemicals) may be better adapted for survival and reproduction under the specific conditions that result from either high or low population densities. For example, when beetle populations are low, individuals may have to fly longer distances to locate rarer and widely distributed hosts rendered susceptible by microclimate, disease, or localized disturbances (Raffa and Berryman 1980). During such periods, beetles possessing a wing size and wing loading that confers greater flight capacity should be at an advantage (Kausrud et al. 2012, Brown et al. 2017, Jones et al. 2019. Likewise, responses to semiochemicals may differ during outbreak and non-outbreak periods (Miller et al. 2005). Differing concentrations and combinations of semiochemicals may signal a host's susceptibility to colonization (Erbilgin et al. 2003, Klutsch et al. 2020, which differs between periods of high or low density of aggressive bark beetle species.
Southern pine beetle (SPB), Dendroctonus frontalis Zimmermann, is an aggressive bark beetle (Coleoptera: Curculionidae: Scolytinae) that occurs in the southeastern, northeastern, and southwestern USA, Central America, and Mexico, and colonizes and kills mature Pinus spp. for feeding and reproduction (Johnson et al. 1982, Hain et al. 2011. During outbreaks, mortality typically occurs in spatially discrete, expanding clusters of contiguous, infested trees called infestations or spots (Thatcher et al. 1980, Price et al. 1992. This native pest has caused significant destruction of pine trees in the southeastern USA, particularly on loblolly (Pinus taeda (L.)) and shortleaf pine (P. echinata Mill.), causing over a billion dollars in economic losses over the past three decades (Clarke andNowak 2009, Pye et al. 2011). In recent years, SPB has expanded its range northward, taking advantage of warmer winter temperatures and attacking novel Pinus spp. hosts in the northeastern USA (Lesk et al. 2017, Dodds et al. 2018. Regional SPB outbreaks in the southern USA occur periodically but are largely unpredictable (Pye et al. 2011, Weed et al. 2017, although factors including weather, predation, and host availability have been extensively investigated as possible drivers of population cycles (Beal 1933, Turchin et al. 1991, Friedenberg et al. 2008, Clarke and Nowak 2009, Reeve 2018. To address this unpredictability, a model for same-year forecasting of SPB outbreaks was developed that utilizes data from a network of pheromone-baited traps deployed each spring throughout the southeastern USA (Billings 1988(Billings , 2011. More recent, ongoing investigations have resulted in the development of a new, updated model for predicting SPB outbreak risk based on current year SPB trap catches, spot abundances the year before last, and the previous year's trap catches of the primary predator of SPB, Thanasimus dubius (Fabricius) (Coleoptera: Cleridae) (Aoki 2017). However, predictions for the data in this study (2016 and 2017) were based on the methodology outlined below. Data on previous year spot abundances and trap catch data of both SPB and its predator T. dubius were entered into the model to generate a rough prediction of spots per county for the coming summer and predict if levels of SPB are increasing, static, or declining (Billings and Upton 2010). Advance warning of an outbreak provides forest managers additional time for preparations such as scheduling of spot detection flights and acquisition of resources necessary for anticipated management efforts. Between 1999 and 2005 in 12 surveyed states, the mean prediction accuracy for outbreak occurrence in a particular county was 82%, and for change in infestation, abundance was 74% (Billings and Upton 2010). However, forecast reliability was low when SPB populations were increasing or high, which is also when accurate foreknowledge is most critical to resource protection. This insufficiency has spurred recent methodology changes aimed at improving efforts to monitor SPB populations, and these include adoption of a revised lure in the monitoring traps (Sullivan et al. 2011) and a longer trapping interval. We posit that monitoring of density-correlated traits of beetles in the trap catches may offer a novel basis for detecting and forecasting SPB population eruptions, potentially augmenting or superseding the data currently used for predicting outbreaks.
Bark beetle sex ratio is a trait subject to variability that is potentially influenced by or correlated with population densities (Lobinger 1996, James et al. 2016. Although the sex ratio of emergent brood SPB has been reported to be approximately 1:1 (Stephen 2011), sex ratios of SPB caught in pheromone-baited traps are typically male-biased and vary considerably (Sullivan 2005, Staeben et al. 2015, and may deviate substantially from the sex ratio of adults in the local population (Hughes 1976, Moser and Browne 1978, Cronin et al. 2000. This disparity is apparently due at least in part to behavioral differences between the sexes in response to the aggregation attractant (Hughes 1976, Moser and Browne 1978, McCarty et al. 1980. One of the authors (BTS) has observed that sex ratios of SPB in pheromone-baited traps vary between years and over relatively long distances, indicating trap-captured sex ratio may be a population-associated trait with interannual fluctuations possibly linked to density. Further, sex ratios of trapped SPB are more male-biased inside SPB infestations than in the surrounding forest (Sullivan et al. 2011), implying that a correlation exists between trapped sex ratio and local abundance of colonized trees.
The objective of our study was to identify traits of SPB trap catches-beyond the number of trapped beetles-that are correlated to population v www.esajournals.org densities, as these might have value for predicting outbreaks. We investigated the quantitative associations between the summertime number of spots per county and two qualities of catches in the spring SPB trapping survey: (1) beetle morphology (body size and hindwing size and shape), and (2) sex ratios. These qualities were chosen because they are easily quantifiable and are likely to be plastic and influenced by population density.

Data
Trapping survey.-For this study, we examined SPB specimens collected during the annual springtime SPB trapping survey coordinated by the USDA Forest Service Forest Health Protection (FHP). The survey uses 12-unit Lindgren multiple-funnel traps (Lindgren 1983) baited with devices releasing the SPB pheromone component frontalin and the pine odors alpha-and beta-pinene Upton 2010, Billings 2011). Starting in 2017, the survey lure also included the pheromone component endo-brevicomin, which increases attractiveness of the lure by approximately an order of magnitude (Sullivan and Mori 2009, Sullivan et al. 2011. Traps were deployed in hardwood inclusions within mature pine stands and at least~10 m from any live pines to avoid initiating an SPB attack (Billings and Upton 2010). Traps used for insect collection in our study were dispersed within boundaries of National Forests, and individual traps (3-6 per National Forest district with minimum 1.6-km spacing) were at approximately the same locations during both years of our study. Trapping was initiated at each location at the start of blooming for eastern redbud (Cercis canadensis L.), and catch collections were made on a weekly basis for 4-6 weeks. A total of 79,982 SPB specimens were collected, identified, sexed (according to Wood 1982), and counted.
Selection of SPB adults.-We used specimens collected in Alabama, Georgia, Mississippi, North Carolina, South Carolina, and Tennessee during 2016 and 2017 to assess morphological variation in SPB (Table 1). We sampled five beetles, at a minimum, from each trap within each year; if fewer than five SPB were collected in a trap, we excluded the trap from the study due to insufficient sample size. We randomly sampled beetles until two of each sex were obtained for each trap/year; the fifth beetle was randomly selected regardless of sex. In traps with fewer than two female SPB, we sampled additional males to obtain five total specimens. For assessment of morphological variation in SPB populations throughout the southeastern USA, we sampled 148 male and 86 female SPB from the 2016 trap catches, and 135 male and 102 female SPB from the 2017 catches (Table 1).
Number of SPB spots per county.-We obtained data on the SPB spots that occurred during the summer through fall of 2016 and 2017 from FHP, which documents the total number and locations of multiple-tree infestations (≥10 trees) at the county level at the end of each calendar year. Southern pine beetle spots are initially detected through repeated aerial surveys that begin in early summer, are identified aerially as clusters of fading or dead foliage in the forest canopy, and are recorded as infestations only after confirmation through ground-checks. We used the number of spots per county as the dependent variable in our analyses to be consistent with SPB outbreak risk prediction and preventative management efforts which are carried out at county level. The SPB specimens we examined were collected from 31 counties in 2016 and 34 counties in 2017. For the counties where we collected SPB samples for our study, the number of SPB spots per county ranged from 0 (numerous counties) to 277 (Franklin County, Mississippi, USA) spots in 2016, and from 0 (numerous Modeling number of SPB spots per county.-For both sampling years, the residuals of spots/ county based on the predictor variables were non-normally distributed, making it impossible to develop prediction models with a linear regression. The high occurrence of zeros in the data (i.e., counties with zero spots), as well as some counties with a high number of spots, are a typical example of a zero-inflated response variable that requires more complex error assumptions. To address this problem, we applied a hurdle model (Cameron and Trivedi 2005) to simultaneously model the likelihood of outbreak occurrence and the magnitude of such an occurrence with separate processes. This modeling framework is tailored to both provide more accurate estimations of error for use in risk assessment and improve fitting properties without violating regression assumptions. The hurdle model is a two-part model that includes a truncated count component, typically Poisson or negative binomial, for modeling positive counts and a hurdle component, either a binomial model or a censored count distribution, for modeling zero counts (Mullahy 1986, Zeileis et al. 2008). The hurdle model features a count data model f count (y|x, β) left-truncated at y = 1 and a zero hurdle model f zero (y|z, γ) right-censored at y = 1 and takes the form (Cameron andTrivedi 2005, Zeileis et al. 2008): where the parameters β, γ, and any additional dispersion parameters θ (for negative binomial distributions) are estimated using a maximum likelihood procedure. The hurdle modeling framework is more flexible than zero-inflated models, as the count component is only used where the hurdle for modeling the occurrence of zeros is exceeded (i.e., only one source of zeros) (Zeileis et al. 2008, Chen et al. 2016. Due to overdispersion of the residuals for spots/county for both years, a truncated negative binomial model provided the best fit for modeling the positive counts (i.e., counties with ≥1 SPB spot), and a binomial model with a logit link for the absence or presence of spots was used for the hurdle (i.e., counties with zero SPB spots; Mullahy 1986, Chen et al. 2016). The negative binomial distribution for y i |x i can be parameterized using the following probability density function (Cameron andTrivedi 2005, Zeileis et al. 2008): where μ is the mean, θ is the shape (dispersion) parameter, Γ(•) is the gamma function, and maximum likelihood is used to estimate all parameters.

SPB sex ratio
SPB spots per county.-Due to the non-normal distribution of the residuals of spots/county, we used hurdle models (Eqs. 1, 2) to assess whether trapped SPB sex ratio was a predictor of the number of SPB spots per county during 2016 or 2017.
SPB caught per trap per day. -We used linear regressions to determine whether SPB sex ratio was a predictor of the mean number of SPB caught per trap per day (i.e., trap catch rates). Linear models take the following form (Seal 1967): where y i is the response variable (SPB/trap/day), β 0 is the intercept, βs are the slope coefficients, x values are the predictor variable(s) (sex ratio), and ϵ i is the error term. The residuals of SPB/ trap/day were non-normally distributed for both years; therefore, we applied a log base 10 transformation to improve the fit of these models (Barrera-Gómez and Basagaña 2015). We ran analyses in R 1.0.136 using the pscl, readxl, and stat packages (Zeileis et al. 2008, Wickham and Bryan 2017, R Development Core Team 2019).

SPB morphological trait variation
We used measurements of three different morphological characteristics (length of right elytron,  (Knapp and Knappová 2013). We measured these features with an ocular micrometer to the nearest 0.1 μm and converted to mm for analyses. We averaged SPB size measurements at the trap-level prior to analyses (n = 48 traps in 2016; n = 47 traps in 2017). We analyzed the two years of data separately due to the addition of the pheromone component endo-brevicomin to the trapping survey lure in 2017, as this has been shown to increase lure attractiveness Mori 2009, Sullivan 2016) and could potentially influence which SPB phenotypes are attracted. We used Pearson's correlation coefficients to evaluate allometry of trap-level elytra length, hind femur length, and hindwing length for both years of samples (Pearson 1895). We ran hurdle models (Eq. 1) to assess whether any of the average trap-level SPB size attributes were predictors of spots/ county. To examine the relationships between body size traits and catch rates, we used linear models (Eq. 3) to regress log-transformed SPB/ trap/day against the mean trap-level size measurements.
To account for any sexual size dimorphism that could bias our assessment of differences in morphological traits with sexes combined, we also averaged individual body size measurements by sex for each trap. We modeled relationships between male and female size measurements and spots/county using hurdle models (Eq. 1), and between male and female size measurements and SPB/trap/day using linear models (Eq. 3; 2016 n = 48 for males, n = 42 for females; 2017 n = 47 for males, n = 46 for females; numbers of replicates varied between sexes because some traps had no female SPB).
To test for sexual dimorphism, we regressed size measurements of individual beetles (n = 234 in 2016; n = 237 in 2017) against sex; all residuals were normally distributed, and thus, linear regression (Eq. 3) was appropriate to model the body size variables. We used dummy coding to include SPB sex as a predictor in these models (i.e., 1 and 0 were assigned for male and female, respectively; Cohen et al. 2003). Individual SPB were the unit of replication for these analyses, as trap-level means of size variables could not be used due to missing data (i.e., several traps had no females). We ran all analyses in R 1.0.136, and used the ggpubr, readxl, pscl, and stat packages (Zeileis

Hindwing landmark and semilandmark acquisition
We utilized landmark-based geometric morphometrics (GM) to quantify variation in hindwing shape among our specimens (Bookstein 1986, Webster and Sheets 2010, Zelditch et al. 2012. The GM technique summarizes shape as a configuration of landmarks (LMs; discrete, homologous points that can be identified in every specimen) that are captured in an image and plotted on Cartesian axes. Geometric morphometrics is a commonly used tool for detecting and assessing subtle morphological differences among individuals and populations (Bookstein 1986, Webster andSheets 2010). The membranous wings in insects (the hindwings in Coleopterans) are particularly well suited for GM; their two-dimensional shape reduces possible distortions due to angle and depth-of-field that are possible with 3-dimensional subjects when photographed (Webster and Sheets 2010), and they possess numerous suitable LMs at the terminations and intersections of veins. Membranous wings have been shown to be a particularly useful feature for identifying morphological variation among insect populations (Oguz et al. 2017, Wilk-da-Silva et al. 2018. We made slide mounts with the right hindwings by compressing them under a coverslip following immersion in 95% ethanol. We took micrographs with a Leica S6D digital camera, and the total hindwing lengths (to the nearest 0.001 mm) were measured using the Leica software's measurement tool (Leica Microsystems, Wetzlar, Germany; Fig. 1b). In studies on hindwing morphology in Coleoptera, LMs have typically been chosen along the costal, radial, medial, and cubitus veins (Bai et al. 2011, Benítez et al. 2014, Su et al. 2015. We selected ten type I LMs (i.e., occurring at the discrete juxtaposition of features; in this study, points of termination or intersection along the medial, cubitus, and radial veins; LM1-9; Fig. 2a), and one type II LM (i.e., geometric constructs including minima/maxima of curves; LM10 was the most distal point of the wing margin; Fig. 2a; Zelditch et al. 2012).
We also assigned semilandmarks (SLMs) in order to provide a shape analysis of the curvature along the outline of the wing (Zelditch et al. 2012, Armendáriz-Toledano et al. 2014b. SLMs are points placed at evenly distributed increments along a curve or outline of a given feature and can be incorporated into a landmark-based analysis to provide a more robust characterization of shape (Webster andSheets 2010, Zelditch et al. 2012). We used MakeFan8 (Sheets 2014) to superimpose over each micrograph a grid of 10 evenly spaced, parallel lines (a comb) oriented perpendicular to a line between LM1 and LM10. We assigned 16 SLMs where the gridlines intersected the anterior and posterior margins of the wings, beginning near the base of the costa and ending near the terminus of the cubitus 1 vein ( Fig. 2b; Hopkins 1909, Sheets 2014. We randomly ordered the hindwing images to reduce bias and then manually assigned the locations of both LMs and SLMs on each digital image to obtain their Cartesian coordinates using tpsDig2 version 2.30 software (Rohlf 2017a).
Analysis of hindwing shape.-To remove nonshape variation (i.e., size, position, and rotation of the hindwings in the micrographs) associated with the raw LM data, we ran a full Procrustes superimposition on these coordinates (Webster and Sheets 2010, Klingenberg 2011, Zelditch et al. 2012, Mikac et al. 2013). We generated a covariance matrix from this data set (Klingenberg 2011, Klingenberg et al. 2012) and derived a correlation matrix for multivariate analyses to account for heterogeneity in the variances.
We ran a principal component analysis (PCA) on the Procrustes coordinates of 27 LMs and SLMs to examine overall patterns of wing shape variation in the sampled SPB populations for 2016 (n = 234) and 2017 (n = 237; Klingenberg et al. 2012). Additionally, we ran a discriminant function analysis (DFA) using the hindwing LMs and SLMs to determine whether these data could partition sampled SPB into low spot density (≤10 spots/county) and high spot density (>10 spots/county) groups based on hindwing shape (Konigsberg and Frankenberg 2018). The DFA used Fisher's cutoff point (i.e., a score of zero) to separate discriminant scores into two groups corresponding to spot density, and Hotelling's T 2 statistic was used to test for significant differences in discriminant scores (Lachenbruch 1967, Klingenberg 2011, Mikac et al. 2013). We performed a permutation test with 10,000 iterations to test the null hypothesis of equal group means, and the reliability of low spot density-high spot density discrimination was determined using leave-one-out cross-validation (Lachenbruch 1967). We applied an identical DFA and error rate estimation to test for the presence of sexual dimorphism in wing shape.
To test for allometry between hindwing shape and length, we ran a multivariate regression of shape against the centroid size (the square root of the sum of squared distances of landmarks from the object centroid) of each hindwing (Klingenberg et al. 2012, Mikac et al. 2013. Scores associated with changes in hindwing shape (i.e., direction of the shape change vectors) were the dependent variable, and centroid size was the independent variable (Drake andKlingenberg 2008, Klingenberg et al. 2012).
We used individual SPB hindwing centroid size as a metric of overall hindwing size to assess whether this attribute was a predictor of spots/ county using hurdle models (Eq. 1). We used logtransformed centroid size to address test assumptions. To determine whether hindwing centroid size differed between male and female SPB, we ran linear regressions (Eq. 3) using sex as the predictor. Hindwing photos (with LMs and SLMs assigned) were consolidated into datasets using tpsUtil32 version 1.74 software (Rohlf 2017b), geometric morphometric and multivariate analyses (including Procrustes superimposition and DFA) were performed in MorphoJ version 1.06d (Klingenberg 2014), and all models and PCAs were executed in R 1.0.136 using the pscl, readxl, and stat packages (Zeileis et al. 2008, Wickham and Bryan 2017, R Development Core Team 2019).

Sex ratio in trap catches
The trap-level sex ratios (female counts/male counts) ranged from 0 to 0.64 in 2016, and from 0.012 to 0.75 in 2017. In 2016, we found no relationship between the captured sex ratio and spots/county (Table 2). However, in 2017, there was a strong negative relationship between trapped sex ratio and SPB spots/county for both components of the hurdle model; traps that had a higher proportion of female SPB were associated with lower numbers of spots (P < 0.001, Table 2; Fig. 3). Similarly, trap-captured sex ratio was not a good predictor of SPB/trap/day in 2016 (Table 3). However, sex ratio increased as the mean number of SPB/trap/day decreased in 2017 (linear model, F 1,45 = 7.14, P = 0.01, Table 3; Fig. 4).

Size variation
In both 2016 and 2017, we found strong positive Pearson correlations among hindwing length, elytra length, and femur length (for all three pairwise correlations, P < 0.001; Table 4). Mean trap-level elytra length, hindwing length, v www.esajournals.org and hind femur length were not significant predictors of spots/county in 2016 or 2017 either with data for sexes combined or analyzed separately (Tables 2, 5). Similarly, for both years, linear models did not show a relationship between the number of SPB/trap/day and the mean traplevel elytra length, hindwing length, or hind femur length with data for the sexes combined or analyzed separately (Table 3).

Hindwing shape variation
In the 2016 PCA on LMs and SLMs, the first five principal components explained 68.7% of the overall variance in SPB hindwing shape. Plots of the first two components displayed nearly complete overlap of specimens from low spot density (≤10 spots/county) and high spot density (>10 spots/county) counties, indicating the absence of a strong association between hindwing shape and infestation abundances (Appendix S1: Fig. S1). The DFA revealed that there were differences in average shape between SPB captured in low and high spot density areas in 2016 (T 2 = 122.0, P = 0.001; Fig. 5a), with hindwings of SPB trapped in low spot density counties being slightly wider and longer. However, a relatively small Procrustes distance (0.00633) between the average hindwing shapes for the two groups indicated minimal overall shape variation corresponding to infestation density, and leave-one-out cross-validation of the DFA resulted in the misclassification of 36% (56 of 155) of low spot density SPB as high, and 44% (35 of 79) of high spot density SPB as low (Fig. 5b).
The DFA analyzing hindwing shape by sex in 2016 showed there were differences in average shape between female and male SPB (T 2 = 138.4, P < 0.001; Fig. 5c), with female hindwings being slightly longer and wider. However, the small Procrustes distance (0.00822) between average hindwing shapes for female and male SPB suggests minimal sexual dimorphism in this feature. Leave-one-out cross-validation misclassified 39.5% (34 of 86) of female SPB as males, and 33% (49 of 148) of male SPB as females (Fig. 5d).
In the 2017 PCA, the first five principal components explained 69% of the overall variation in hindwing shape. As in 2016, graphical assessment of the principal components indicated extensive overlap in the hindwing shapes of SPB collected from counties with high and low spot densities (Appendix S1: Fig. S2). The DFA revealed differences in average hindwing shape between SPB captured in low and high spot density counties in 2017 (T 2 = 171, P < 0.001; Fig. 6a), with SPB trapped in low spot density counties having slightly wider and longer hindwings; however, the Procrustes distance (0.00323) between the average shapes was small. Cross-validation resulted in the misclassification of 29.5% (38 of 129) of low spot density SPB as high, and of 33% (36 of 108) of high spot density SPB as low (Fig. 6b).
For 2016, multivariate regression analysis assessing allometry within the hindwing data revealed that centroid size explained 25.2% of the total variation in hindwing shape (R 2 = 0.25; Fig. 3. Relationship between southern pine beetle spots per county (all data) and trap-level sex ratio (female counts/male counts) in 2017. Plot includes a smoothed line of fit (loess method) and 95% confidence interval. There was no relationship in 2016 (P > 0.05). Notes: SPB, southern pine beetle. Significant results are shown in bold (α ≤ 0.05).
† Data log-transformed for analyses. ‡ Morphological measurements averaged at the trap level; sex ratio has trap as unit of replication.
v www.esajournals.org P < 0.001). The linear regression of centroid size on SPB sex did not show a difference in hindwing centroid sizes between males and females in 2016 (P = 0.32). Using the hurdle model with positive spot counts (i.e., for y > 0; Eq. 1), we detected a marginally significant relationship between spots/county and hindwing centroid size, with SPB having slightly larger centroid sizes in counties with lower numbers of spots (P = 0.05; Table 2), but this relationship was not significant (P = 0.49) for counties with zero spots (i.e., for y = 0; Eq. 1).
For 2017, hindwing centroid size explained 35.2% of the total variation in hindwing shape (R 2 = 0.35; P < 0.001). Contrary to 2016, we found differences in centroid size between the sexes in 2017 (F 1,235 = 5.44, P = 0.02); hindwing centroid size for males was 98% as large as female centroid size (male: 1802.6 AE 10.1, female: 1840.1 AE 12.8). The hurdle model showed no relationship between SPB spots/county and hindwing centroid size in 2017 (Table 2).

Relationship between SPB trait variation and the number of SPB spots per county
To our knowledge, there are no studies comparing multiple phenotypic traits of SPB adults collected from areas with differing levels of SPB damage, although variation has been observed in attributes of SPB and other Dendroctonus species that vary with population density (e.g., female SPB body size, oviposition rate, gallery length and density, and development time; Bentz et al. 2011, Stephen 2011, Foelker and Hofstetter 2014. However, we did not find compelling evidence for variation in trapped SPB body size (elytra length, hind femur length, or hindwing length) that was associated with the number of spots in surrounding forests. The size of our data set may have limited our ability to detect density-dependent variation. We had only two years of SPB trapping survey data available, and because the trap lure components differed between 2016 and 2017, we were unable to integrate morphological trait data and infestation numbers from these two years. However, the very small differences in size parameters associated with areas with high and low spot densities do not suggest that these differences are biologically relevant or adaptive.

Relationship between SPB trait variation and SPB caught per trap per day
The number of SPB captured in baited funnel traps has been shown to be significantly but inconsistently correlated with other indicators of local population density such as proximity of attacked trees (Sullivan et al. 2011; also recently demonstrated in the congeneric D. ponderosae, Klutsch et al. 2020). Raw SPB catches are broadly correlated with spot abundances and used currently as one indicator for forecasting SPB damage (Billings and Upton 2010). Intuitively, Table 4. Correlation coefficients (r), t statistics, and P values from Pearson's correlations between measured body size traits averaged at the trap level.  springtime trap catches should be associated with SPB damage potential, since a threshold density of dispersing, pheromone-responding SPB is necessary for beetles to organize mass attacks to overcome host defenses and initiate spots (Stephen 2011, Schowalter 2012). However, we did not find a relationship between SPB body size attributes and trap catch numbers, providing corroborating evidence that beetle size is not a useful indicator of beetle population density and damage potential.

Relationship between hindwing shape and the number of SPB spots per county
Although we found in 2016 that SPB hindwing centroid size was correlated with spot numbers for counties that had at least one spot (Table 2), this relationship was only marginally significant and did not hold when counties with zero spots were included in the analysis. Combined with the absence of any trend in 2017, these data imply that wing centroid size would not be useful in predicting outbreaks and is similar to our 15 † Trap-level data for both sexes. ‡ None of the models were statistically significant (α ≤ 0.05). § The square root of the sum of squared distances of hindwing landmarks from the wing centroid; it is a measure of overall wing size. Data were log-transformed.
other size measurements in this respect. There were subtle average hindwing shape differences between SPB trapped in low vs. high spot density counties in 2016 and 2017, with hindwings of SPB from low spot density counties being slightly wider and more elongated. However, cross-validation indicated that hindwing shape scores were not a reliable metric for distinguishing SPB captured from low vs. high spot density counties in either year (63% and 70.5% accuracy in classifying low spot density SPB, and 56% and 67% accuracy in classifying high spot density SPB in 2016 and2017, respectively;Figs. 5b, 6b). It is important to note that to assess differences in hindwing shape using DFA, the trapped SPB samples were partitioned into two, somewhat arbitrary groups (counties with >10 or ≤10 spots). However, the lack of clustering of beetles into distinct groups in the PCA indicates the absence of polymorphism (Appendix S1: Figs. S1, S2), and the strong overlap of clusters for our two spot density classes suggests that alternative partitioning of spot densities would likely not have permitted stronger group separation. Additionally, the Fig. 5. Score frequencies from discriminant functions calculated to separate beetles trapped in 2016 as being from either (a, b) high spot density (scores < 0) or low spot density counties (scores > 0) or (c, d) female (scores < 0) or male (scores > 0) based on hindwing shape (27 landmarks and semilandmarks). High and low spot density counties had >10 or ≤10 spots, respectively. In subfigures a and c, scores are calculated with a single function derived from the entire data set, whereas subfigures b and d display scores calculated from a leave-oneout cross-validation.
v www.esajournals.org small difference in average shape and significant overlap in hindwing shape scores between SPB captured in lower vs. higher spot density locations (Figs. 5a, 6a) suggests that this phenotypic trait would have minimal if any utility in detecting shifts in population density.

Sexual dimorphism
Our results provide support for the existence of sexual size dimorphism in SPB, with female beetles being slightly larger than males (longer elytra in both years; longer hindwings and larger hindwing centroid sizes in 2017). Results from the DFA assessing hindwing shape variation also showed subtle differences between male and female SPB, with female hindwings being slightly more elongated and wider in both years as compared to males. However, the cross-validation on the hindwing shape discrimination indicated that shape was not a consistent indicator of sex (60.5% and 65% accuracy in classifying female SPB, and 67% and 63% accuracy in classifying male SPB in 2016 and 2017, respectively; Figs. 5d, 6d). Previous studies have shown that Fig. 6. Score frequencies from discriminant functions calculated to separate beetles trapped in 2017 as being from either (a, b) high spot density (scores < 0) or low spot density counties (scores > 0) or (c, d) female (scores < 0) or male (scores > 0) based on hindwing shape (27 landmarks and semilandmarks). High and low spot density counties had >10 or ≤10 spots, respectively. In subfigures a and c, scores are calculated with a single function derived from the entire data set, whereas subfigures b and d display scores calculated from a leave-oneout cross-validation.
v www.esajournals.org female beetles, on average, are larger than males for several species in Dendroctonus (Bentz et al. 2011, Graf et al. 2012, Foelker and Hofstetter 2014, Lachowsky and Reid 2014, Liu et al. 2017. For example, there may be greater selective pressures on females than males for larger size, as the females initiate attacks on a new host and need adequate resources to locate the best environment for brood feeding and development (Elkin andReid 2005, Bentz et al. 2011). The successful detection of expected sexual dimorphism with our data indicated the suitability of our morphometric procedures for exploring density-dependent variation in morphology. Furthermore, the very small differences between the sexes in our examined morphological traits suggest that it was unlikely that sex was a confounding factor in detecting morphological correlates of spot densities.

Captured SPB sex ratio
For one year of trapping data (2017), female SPB were trapped in substantially higher proportions in counties with lower infestation densities; the average trap-captured sex ratio was 0.294 in low spot density counties (≤10 spots/county) vs. 0.097 in high spot density counties (>10 spots/ county). This trend was absent in 2016. The major procedural difference between years was the addition of endo-brevicomin to the trapping survey lure in 2017. Although endo-brevicomin has not influenced captured sex ratios of SPB in previous tests (i.e., it enhances attraction of both sexes similarly; Mori 2009, Sullivan 2016), it is possible that an undiscovered, density-dependent sexual dimorphism in effects of endo-brevicomin exists. Additionally, the average number of SPB spots per county was substantially higher in 2017 (~211 spots/county) than 2016 (~45 spots/county), and it is possible that a trend in trap-captured sex ratios was not detectable in 2016 due to an insufficient quantity of data from counties with high spot densities. Other studies have indicated an association between sex ratios in traps and local population densities in Dendroctonus. Catches in traps located inside of active SPB infestations (where beetle densities were ostensibly high) were substantially more male-biased than catches in traps deployed 100 and 200 m away among uninfested trees (Sullivan et al. 2011). In one of two years and with specific lure combinations, higher population densities of D. ponderosae were associated with a higher proportion of males in traps (Klutsch et al. 2020), as was observed with SPB in this study.
In 2017, we found that the proportion of males in traps increased with increasing trap catch rates (Fig. 4). This relationship is consistent with the correlation between spot abundance and trapped sex ratios observed in 2017 (Table 2, Fig. 3), since trap catches of SPB and spot abundances are known to be correlated, as discussed previously. Density dependence in sex ratios of trap-captured SPB could result from changes in the sex ratio within the dispersing population or in the sex-specific propensities of beetles to be caught in traps; however, these possibilities cannot be distinguished with our data. Sex ratios of emerging brood adult and reemerging parent adult SPB have not been found to differ significantly from 1:1 (Cooper and Stephen 1978, Coulson et al. 1979, Cook and Hain 1985. This suggests that skewed sex ratios in brood production or sexual differences in mortality rates in the host could not be factors contributing to strongly skewed sex ratios in the dispersing population or in trap catches. During outbreaks, traps are competing with a greater number of sources of natural attractant (i.e., mass-attacked trees), and traps and attacked trees likely differ in their relative attractiveness to the sexes. Male and female SPB differ in their responses to some semiochemicals. Two male-produced compounds in SPB, verbenone and acetophenone, reduce attraction of males more than females (Salom et al. 1992, Sullivan 2005, and likely other such semiochemicals exist. The odor blends of trap lures and trees are substantially different (Pureswaran and Sullivan 2012). Hence, if females have a stronger attraction than males to mass-attacked trees relative to traps, the higher numbers of these trees associated with outbreaks could produce male-skewed trap catches. Additionally, a semiochemical background in the environment can alter SPB responses to baited traps (B. T. Sullivan, unpublished data), and increased background of SPB semiochemicals within outbreak areas could affect male and female responses to trap lures differently.

Applications of findings for management and monitoring of SPB
The trend we documented in 2017, with lower proportions of female SPB in traps occurring in v www.esajournals.org counties that experienced higher numbers of beetle infestations, has important implications for managing SPB and predicting outbreaks. Inclusion of trap-captured sex ratios in the total data set currently used for outbreak forecasting (i.e., SPB trap catches, T. dubius catches, and spot densities of previous years) could contribute to the accuracy of these forecasts at local and regional levels. This potential will be assessed with trial collections of sex ratio data during the annual survey for the next several years. If implemented as part of the survey, collection of sex ratio data would likely require a more efficient means of sexing large numbers of insects than inspection of individuals under a dissecting microscope. However, relatively affordable, high-throughput molecular techniques such as qPCR could be used to determine sex ratios of beetles within homogenates of catch samples by quantifying proportions of sex-linked genes (Belousova et al. 2019). Hence, it may potentially be feasible to incorporating sex ratio data in SPB outbreak forecasting within the ongoing trap catch analyses, and these data may potentially contribute to longer-term sustainability and conservation of resources in the southeastern USA.

ACKNOWLEDGMENTS
We acknowledge the financial support of the D.B. Warnell School of Forestry and Natural Resources, the USDA Forest Service Forest Health Protection, and the USDA Forest Service Southern Research Station. We thank Ronald Billings, Anthony Elledge, Wood Johnson, James Meeker, and Richard Spriggs (USDA Forest Service Forest Health Protection) for providing the samples, spatial data with locations of SPB infestations, and answering many inquiries related to trapping protocol and sampling methodology. JoAnne Barrett (USDA Forest Service Southern Research Station) assisted with the morphometrics work, and Brittany Barnes (University of Georgia) helped with the Leica microscope and software. BHM designed the experiment, collected data, performed analyses, and was the primary author of the manuscript text. BTS assisted with conceptual background, experimental design and procedure, geometric morphometric analyses, and refining scope and structure of the manuscript. HLM and CRM contributed to statistical analyses and data interpretation. JTN provided background on SPB trapping survey and assistance in developing SPB sampling scheme. CV helped to refine manuscript text. KJKG assisted with the development of conceptual background, experimental design, and structuring and editing the manuscript.