When gene flow accompanies speciation, recombination can decouple divergently selected loci and loci conferring reproductive isolation. This barrier to sympatric divergence disappears when assortative mating and disruptive selection involve the same “magic” trait. Although magic traits could be widespread, the relative importance of different types of magic traits to speciation remains unclear. Because body size frequently contributes to host adaptation and assortative mating in plant-feeding insects, we evaluated several magic trait predictions for this trait in a pair of sympatric Neodiprion sawfly species adapted to different pine hosts. A large morphological dataset revealed that sawfly adults from populations and species that use thicker-needled pines are consistently larger than those that use thinner-needled pines. Fitness data from recombinant backcross females revealed that egg size is under divergent selection between the preferred pines. Lastly, mating assays revealed strong size-assortative mating within and between species in three different crosses, with the strongest prezygotic isolation between populations that have the greatest interspecific size differences. Together, our data support body size as a magic trait in pine sawflies and possibly many other plant-feeding insects. Our work also demonstrates how intraspecific variation in morphology and ecology can cause geographic variation in the strength of prezygotic isolation.

Geographic variation in host needle width and adult body size: Because Neodiprion eggs must be completely embedded in pockets that adult females carve within pine needles, egg size and female ovipositor size are likely constrained by needle width (Bendall et al. 2017). Given these constraints, we predicted that geographic variation in adult body size would correlate with geographic variation in host needle width. Although the ideal way to test this prediction would have been to collect host needles and sawfly adults from the same sites, we only had access to pine and sawfly samples that were collected at different times and without this hypothesis in mind. We therefore analyzed geographic variation in host needle width and adult body size separately, with the following expectations: N. lecontei adults and host needles will be bigger than those of N. pinetum and geographic trends in body size would mirror those observed for needle width (i.e., have slopes with the same sign and similar magnitude).

To characterize geographic variation in pine needle width, we sampled needles from multiple trees and sites for each of 10 Pinus species used by N. lecontei (P. taeda, P. palustris, P. echinata, P. elliottii, P. clausa, P. glabra, P. virginiana, P. rigida, P. resinosa, and P. banksiana) and the single host species used by N. pinetum (P. strobus). For each pine species, we collected clippings from 60-100 trees sampled from 6-10 sites (9-10 trees per site). Using digital calipers (Mitutoyo CD-6’PMX), we measured the width of 3 randomly sampled needles per tree and averaged these measurements to obtain a single needle width value per tree. In total, we sampled 878 individual pine trees from 88 sites (Table S1). To determine how the needle width varies as a function of host species and latitude, we used linear regression to model needle width as a function of pine species, latitude, and their interaction. We used a Type III Analysis of Variance (ANOVA), implemented in the car v3.1-0 package (Fox and Weisberg 2019), to evaluate the significance of model terms and the emmeans v1.8.0 package (Lenth 2021) for post-hoc tests with false discovery rate (FDR) correction for multiple testing. These and all other statistical analyses were performed in R version 4.1.0 (R Core Team, 2021). Because we were specifically interested in comparing body size clines between sawflies with needle width clines from their pine hosts, we also used linear regression to estimate the relationship between (1) needle width and latitude for all N. lecontei hosts (regardless of host species) and (2) needle width and latitude for P. strobus, the only N. pinetum host pine.

To characterize geographic variation in body size for adult females and males of both species, we collected mid- to late-instar N. lecontei and N. pinetum larvae across the eastern United States between 2015 and 2021 (Table S2 and Fig. S1) and reared the immature stages to adults in the lab using host plant clippings and standard lab protocols (Harper et al. 2016; Bendall et al. 2017). Upon emergence, live adults (which are non-feeding) were either preserved immediately or stored at 4°C to prolong life until needed for propagating lab lines or use in experimental assays. All adults were ultimately placed in 100% ethanol and stored at -20°C. We used a Neiko Tools Digital Caliper (model 01407A) to measure body length (i.e., from the tip of the head to the tip of the abdomen) to the nearest hundredth of a millimeter (mm) for 1,080 preserved adults (N. lecontei females: N = 328; N. lecontei males: N = 243; N. pinetum females: N = 279; N. pinetum males: N = 230; Table S2). To determine how body length differed as a function of latitude, species, and sex, we fit a linear model to the body size data, with latitude, species, sex, latitude x species, latitude x sex, and species x sex as fixed effects. We used a Type III ANOVA to assess the significance of model terms. Based on these results, we also fit individual geographic clines to each sex for each species.

Divergent selection on egg size: We hypothesize that variation in needle width among pine populations and species generates divergent selection on body size (egg size) in Neodiprion that use different pine hosts via a combination of constraints imposed by thin needles (favors smaller eggs) and selection on early larval survival (favors larger eggs). To test this hypothesis, we first verified that N. pinetum (thin-needled specialist) and N. lecontei (uses hosts with thicker needles) differ in egg size. To do so, we used females reared from wild-caught larvae collected from different sites in Kentucky and Tennessee between 2013 and 2015 (Table S3). In Neodiprion, all egg maturation occurs within the cocoon, from which females emerge with a full complement of eggs ready for oviposition (i.e., these species are pro-ovigenic). To measure egg length and width, we dissected eggs out of the abdomens of recently eclosed adult females and photographed a random subset of eggs at 10X total magnification using a Zeiss Discovery V8 stereomicroscope with an Axiocam 105 color camera and ZEN lite 2012 software (Carl Zeiss Microscopy, LLC Thornwood, NY). We then used the ZEN lite software to measure the length and width of 5 eggs from each of 5 females from each species (N = 25 eggs per species). We then calculated egg area after García-Barros (2000) using the following equation:

          egg size = (0.5236 x d² x h)^1/3,

where d = egg diameter (width) and h = egg length. To test for differences in egg size between N. pinetum and N. lecontei, we fit a linear mixed-effects model to the egg area data, with species as a fixed effect and female ID as a random effect (since multiple eggs were sampled from each female). To evaluate the significance of the species term, we used a Type II ANOVA.

N. lecontei and N. pinetum differ in many traits that affect hatching success on different types of pines, including ovipositor size and number of eggs per needle (Bendall et al. 2017; Fig. 3A, 3B). To evaluate the effect of egg size independent of other traits, we used a backcross design to generate recombinant hybrid progeny between N. lecontei and N. pinetum after Bendall et al. (2017). All females were lecontei-backcross females (one set of N. lecontei chromosomes and one set of chromosomes with approximately 50% N. lecontei and 50% N. pinetum ancestry). We used a lecontei-backcross design because preference for the thin-needled N. pinetum host is dominant (Bendall et al. 2017), and we wanted to ensure a good sample size of females that laid on both hosts. Backcross females were released individually into mesh cages (60 x 40 x 40 cm) containing two seedlings of Pinus banksiana (N. lecontei host) and two seedlings of P. strobus (N. pinetum host). Cages were checked daily for oviposition and once this occurred, live females (which tend to remain at the bottom of the egg-bearing branch until death) were preserved in 100% ethanol for subsequent morphological and molecular work. After preserving females, the number of eggs laid on each seedling was counted and monitored daily for hatching. Because females almost always cluster their eggs in a single branch tip, host choice was scored as a binary trait (P. strobus or P. banksiana). Once a hatchling was observed, eggs were given an additional 48 hours to provide sufficient time for hatching. All hatchlings were then removed by washing them off the pine seedling in ethanol. Hatchlings were counted and the hatching success for each female was calculated as the proportion of eggs that hatched (number of hatchlings/number of eggs).

After oviposition, females almost always have eggs left in the abdomen. We used these remaining eggs to quantify the average egg size for each female. To measure egg size, we first rehydrated preserved female abdomens by soaking them in five decreasing concentrations of ethanol (95%, 80%, 65%, 50%, and 25%) for 10 minutes each. After the lowest ethanol concentration, we soaked females in deionized water for 24 hours at room temperature. We then dissected the eggs from each female’s abdomen and placed them in 300 mL of a modified Ringer’s dissection solution (7.5 g/L NaCl, 0.35 g/L KCl, and 0.21 g/L MgCl2). We then imaged 1-10 eggs per female at 10X total magnification and measured eggs and calculated the area as described above. We then averaged these values to obtain a single average egg area per female. Although ethanol preservation and rehydration may cause preserved eggs to differ in size from fresh eggs, this approach should nevertheless reveal relative egg size differences among females (i.e., which females laid the largest or smallest eggs) because all eggs were treated in the same way. We are also assuming here that the size of leftover eggs correlates positively with the size of laid eggs, which will require future experiments to confirm.

In total, we measured hatching success, egg size, and host preference for N = 38 lecontei-backcross females. After removing two females that laid eggs on both pine species, our final sample size was N = 36 backcross females. To evaluate the effect of egg size (egg area) and host (P. banksiana or P. strobus) on hatching success, we used the glm function (lmerTest v3.1-3; Kuznetsova et al. 2017) to fit a logistic regression model to the hatching data (proportion hatched ~ egg area + host + area*host), followed by a Type III ANOVA to evaluate the significance of model terms. A significant egg area x host interaction would support our hypothesis that there is a trade-off between egg fit and egg provisioning.

Size-assortative mating and reproductive isolation: To determine whether Neodiprion adults mate assortatively by size and to quantify reproductive isolation between N. pinetum and N. lecontei, we first sampled larval colonies of both species from several locations throughout the eastern United States from June to August 2019 (Table S4 and Fig. S2). In both species, females tend to mate once and then lay a single clutch of eggs in one pine branch terminus. Thus, distinct larval clusters typically represent full-sib families. To maintain broad-scale geographic differences in behavior or morphology, we grouped field-collected larvae by state (Michigan, Kentucky, Indiana, and North Carolina) and then reared larvae to adults for mating assays and line propagation. Methods for line propagation are described elsewhere (Harper et al. 2016; Bendall et al. 2017). To minimize the evolutionary change in the lab, we used adults that were reared either from wild-caught colonies or first-generation lab colonies.

No-choice mating assays were performed from September to October 2019. We used no-choice assays because they are consistent with mating behaviors in the wild (Benjamin 1955; Wilson et al. 1992). Because reproductive isolation can vary across geographic space (Jiang et al. 2013; Rougemont et al. 2015), we assayed mating outcomes in three different crosses, each containing a different combination of N. lecontei and N. pinetum adults from different U.S. states: Cross 1 = North Carolina N. lecontei x Indiana N. pinetum; Cross 2 = Kentucky N. lecontei x Kentucky N. pinetum; and Cross 3 = Kentucky N. lecontei x Michigan N. pinetum. Each assay consisted of two arenas (plain white 8 x 11 pieces of printer paper) that were divided into six equally sized sections (dividing lines drawn with a black Sharpie pen; Fig. S3). Each arena was recorded by either a Logitech Carl Ziess Tessar HD 1080p or Microsoft LifeCam Cinema Model 1393 camera. We placed a small petri dish (5 cm x 1.5 cm) in each section of each arena, within which we placed a single male and female. Each assay (pair of arenas) consisted of three replicates of each of the four possible female-male pair types for a particular cross: (1) N. lecontei female x N. lecontei male, (2) N. lecontei female x N. pinetum male, (3) N. pinetum female x N. lecontei male, and (4) N. pinetum female x N. pinetum male (i.e., each assay contained 12 female-male pairs total). Based on the availability of males and females from each geographic region, we were able to complete 12 assays for Cross 1 (N = 36 replicates per pair type; 144 female-male pairs total), 5 assays for Cross 2 (N = 15 replicates per pair type; 60 female-male pairs total), and 5 assays for Cross 3 (N = 15 replicates per pair type; 60 female-male pairs total). 

During each assay, which lasted 2 hours, all mating events were recorded for each male-female pair. Each pair was assigned an arbitrary identifier and observed blind with respect to pair type. For an interaction to be considered a mating, the pair in question had to be properly aligned and physically connected (Fig. 1C) for at least 1 minute to ensure a sufficiently secure attachment for sperm transfer. Although we have not experimentally determined the minimum duration that a pair must be connected for sperm transfer to occur, previous observations have shown that sawfly pairs that mate for approximately 1 minute produce daughters (Author, personal observation), which indicates that fertilization has occurred as pine sawflies are haplodiploid (Knerer and Atwood 1973). At the end of each assay, the overall mating outcome of each pair was recorded as a “1” if the pair mated at least once or as a “0” if the pair never mated.  All males and all unmated females were immediately preserved in 100% ethanol. Mated females were allowed to lay eggs for further lab line propagation but were preserved live after egg laying (except for 17 females that we were unable to find). To confirm recorded mating outcomes, we reviewed all videos (again, blind with respect to each pairs’ identity) and scored each pair as described above. While reviewing videos, we also counted the number of mating attempts made by each male. For a behavior to be counted as a “mating attempt,” we required that the male curled his abdomen in a “U-shape” under the female’s abdomen in an attempt to initiate copulation. We then coded male behavior in each pair as a binary trait: “1” if the male made any attempt and “0” if the male made no attempts.

We used our mating assays to quantify the strength of prezygotic isolation between N. lecontei and N. pinetum both within each of the three crosses and globally (all crosses combined), following Sobel and Chen (2014). This method requires observed and expected mating frequencies for conspecific and heterospecific pairs. To obtain the observed frequency of conspecific matings, we divided the total number of conspecific pairs that mated (either N. lecontei female x N. lecontei male or N. pinetum female x N. pinetum male) by the total number of attempted conspecific crosses. To obtain the observed frequency of heterospecific matings, we divided the total number of heterospecific pairs that mated (either N. lecontei female x N. pinetum male or N. pinetum female x N. lecontei male) by the total number of attempted heterospecific crosses. For the expected conspecific and heterospecific mating frequencies under random mating, we used 0.5. We then calculated reproductive isolation using the following equation:

            RI = 1 – 2 x (H/H+C);

where H = heterospecific pairs and C = conspecific pairs. This equation yields reproductive isolation values ranging from 0 (no reproductive isolation) to 1 (complete reproductive isolation). To further evaluate evidence of prezygotic isolation, we used one-tailed Fisher’s exact tests in each cross and globally (fisher.test function) to determine whether heterospecific pairs were significantly less likely to mate than conspecific pairs. The Sobel and Chen (2014) equation and our statistical analysis ignores possible asymmetries within pair types. Therefore, to determine whether mating outcomes differed between the four different pair types (N. lecontei female x N. lecontei male; N. lecontei female x N. pinetum male; N. pinetum female x N. lecontei male; N. pinetum female x N. pinetum male), we also used the fisher.multcomp function from the RVAideMemoire v0.9-81-2 package (Hervé 2022), with FDR correction for multiple testing. For plotting proportion data, we used the DescTools v0.99.45 package (Signorell et al. 2022) to calculate 95% Clopper-Pearson confidence intervals for each cross and pair type.

To determine whether body size influenced mating outcomes, we measured body size for N = 511 preserved males and females from our mating assays. Because body length is likely to be especially important for proper alignment during mating (Fig. 1C and Video S1), we used body length as our measure of body size. Specifically, we used a Neiko Tools Digital Caliper (model 01407A) to measure each individual from the tip of the head to the tip of the abdomen to the nearest hundredth of a millimeter. To reduce measurement error, each sawfly was measured independently by two individuals, and the average body length (in mm) was used. To calculate a size differential for each pair, we subtracted the male body length from the female body length.

We then used logistic regression to evaluate how species, pair type, and size differences affect binary mating outcomes in each of the three crosses and globally. For each cross, we modeled mating outcome as a function of female species, male species, female species x male species interaction, and size differential. For the global model, cross was included as a random effect. Then, we used Type III ANOVAs to evaluate the significance of model terms. If the two species differ in willingness/motivation to mate, we expected significant male species or female species terms. If the type of male-female pair (conspecific versus heterospecific) affects mating outcomes independent of body size, we expected significant female species x male species interaction terms. And if there is size-based assortative mating independent of male-female pair type, we expected significant effects of size differentials in our models, with pairs that mated having a smaller size differential than those that did not. We also used logistic regression to evaluate how species, pair type, and size differences affect male mating attempts (coded as a binary trait) in each cross and globally. Then, we used Type III ANOVAs to evaluate the significance of model terms.

To further explore how variation in body size within and between species relates to variation in the strength of prezygotic isolation among the three crosses, we compared the size differences for all intraspecific pairs (N. lecontei female x N. lecontei male and N. pinetum female x N. pinetum male) to the size differences for all interspecific pairs (N. lecontei female x N. pinetum male and N. pinetum female x N. lecontei male) within each cross and globally. We then performed one-sided t-tests using the t.test function to evaluate whether the size differential for interspecific pairs was significantly greater than the size differential for intraspecific pairs.

All data files are .csv files and can be opened in R.