Patterns of morphological divergence across species’ ranges can provide insight into local adaptation and speciation. In this study, we compared phenotypic divergence among 4,221 crickets from 337 populations of two closely related species of field cricket, Gryllus firmus and G. pennsylvanicus, and their hybrids. We found that these species differ across their geographic range in key morphological traits, such as body size and ovipositor length, and we directly compared phenotype with genotype for a subset of crickets to demonstrate nuclear genetic introgression, phenotypic intermediacy of hybrids, and essentially unidirectional mitochondrial introgression. We discuss how these morphological traits relate to life history differences between the two species. Our comparisons across geographic areas support prior research suggesting that cryptic variation within G. firmus may represent different species. Our study highlights how variable morphology can be across wide-ranging species and the importance of studying reproductive barriers in more than one or two transects of a hybrid zone.

Gryllus, hybrid zone, local adaptation, morphology, Orthoptera, speciation

Phenotypic divergence can provide insight into evolutionary processes acting across different scales of biological organization. Within a single species, phenotypic divergence can reflect differences between environments, population histories, or a combination of these factors (Gavrilets et al. 2001, Uyeda et al. 2009, Runemark et al. 2010, Oneal and Knowles 2013, Jenck et al. 2020). Phenotypic divergence can signal the possible early stages of species differentiation (Wolf et al. 2008, González et al. 2011, Skoglund et al. 2015) and, in closely related species, can shed light on local adaptation and patterns of increasing divergence (Britch and Cain 2001, Shaw and Mullen 2011). Most studies of species divergence have limited replication across the ranges of a species pair, and the specific traits that maintain reproductive barriers between species are not always clear (Harrison and Larson 2016). Geographically comprehensive surveys of phenotypic divergence are much harder (Jiménez and Ornelas 2015, Wang et al. 2017, Polly and Wójcik 2019, Moran et al. 2020) but critical if we are to understand the origin and maintenance of species’ boundaries.

The relationship between divergent phenotypic characteristics and reproductive barriers is most easily studied in places where the ranges of closely related species overlap and heterospecific individuals mate and produce offspring (Barton and Hewitt 1985, Harrison 1990). In the resulting hybrid zone, as the different species co-exist, compete, and interbreed, phenotypic characteristics may be more variable among individuals when compared to the pure allopatric populations that lie outside the hybrid zone (Hollander et al. 2018, Sottas et al. 2018). By comparing the phenotypic variation between both conspecific allopatric and sympatric populations and between heterospecific populations, it becomes possible to examine the potential causes of phenotypic evolution, speciation, and how those mechanisms lead to the reproductive barriers that maintain species boundaries (Shaw and Mullen 2011).

In this study, we examined the phenotypic divergence between two closely related and geographically widespread species of North American field crickets, Gryllus pennsylvanicus Burmeister 1838 and G. firmus Scudder 1902, whose common ancestry dates to roughly 200,000 years ago (Willett et al. 1997, Maroja et al. 2009a). The more northern, inland species, G. pennsylvanicus, is broadly distributed throughout the United States, while the more southern, coastal species, G. firmus, is restricted to the east coast and west into Texas (Alexander 1968, Harrison and Arnold 1982, Weissman and Gray 2019). These species form a hybrid zone along the eastern front of the Appalachian Mountains (Harrison and Arnold 1982), and where they co-occur, they are isolated by multiple reproductive barriers. The most striking barrier is a one-way incompatibility: G. firmus females mated to G. pennsylvanicus males lay few eggs that do not hatch (Harrison 1983, Maroja et al. 2009b, Larson et al. 2012). These two species are also isolated by habitat: G. firmus is often found in sandy habitats and has lighter coloration and longer ovipositors that can presumably lay eggs deeper in sandy soils (Harrison 1986, Ross and Harrison 2006). Also, Gryllus firmus is a larger cricket, though size may vary with the length of the growing season (Masaki 1961). In some parts of the hybrid zone, G. firmus develops faster and emerges earlier in the season, leading to temporal isolation (Harrison 1985).

These morphological differences have been well characterized in a handful of locations within the hybrid zone (e.g., Connecticut), but whether these morphological traits are consistently different between G. firmus and G. pennsylvanicus remains an open question (Weissman and Gray 2019). When species differences are studied in only a few locations, it may be impossible to distinguish species-specific traits from within-species local adaptation. Morphological traits, such as lighter color and longer ovipositors, may have evolved in specific areas due to habitat selection. Likewise, body size may vary with climate and latitude. This paper presents the first geographically comprehensive comparison of G. firmus and G. pennsylvanicus by combining published and unpublished morphological datasets for these two species across their geographic ranges. Our dataset includes 4,221 crickets from 337 populations, spanning collections over four decades. We had three objectives. First, we quantified morphological divergence within and between species across their geographic ranges. Second, for populations near the hybrid zone, we tested whether traits that distinguish species correlate with ancestry. Finally, we examined the correlation between morphological traits and environmental variables across the ranges of these species. In doing so, we aimed to gain a greater understanding of how population variation and local adaptation contribute to divergence and speciation.

Cricket collections.—We compiled a dataset of 4,221 crickets, the majority being G. pennsylvanicus but also G. firmus and their hybrids, from 337 collecting localities (Fig. 1). Crickets were sampled throughout the United States and Canada, with the largest collections coming from the northeastern United States and the hybrid zone. Sampling spanned 40 years (1983–2022), with collections performed by A.R. Byerly, E.L. Larson, L.S. Maroja, C.L. Ross, and R.G. Harrison. In addition to these previously unpublished morphological data, we included data from Ross and Harrison (2002), Larson et al. (2013), and Weissman and Gray (2019), with the latter being the most geographically widespread dataset. We also included morphological data from a newly described cricket species, G. thinos Weissman and Gray 2019, , which is closely related to G. pennsylvanicus and G. firmus (Gray et al. 2020). We included G. thinos to enable us to compare morphological variation within G. firmus to that of a closely related species that occupies the same habitat but is classified as a separate species.

Fig. 1. 

Map of North American cricket collecting locations. Allopatric populations of Gryllus firmus are in yellow, G. pennsylvanicus are in teal, G. thinos populations are in purple, and sympatric G. firmus and G. pennsylvanicus populations are in red. The size of the circle corresponds to the sample size for each location. A. Entire range of collection locations in the United States and Canada; B. Enlarged area of densely sampled locations in northeast, central, and southeast United States.

We categorized each collecting location as allopatric or sympatric based on past sampling of the field cricket hybrid zone (Harrison and Arnold 1982, Willett et al. 1997, Maroja et al. 2009a, Larson et al. 2013a, 2014). Populations in and near the hybrid zone often have individuals that are pure G. firmus or pure G. pennsylvanicus, but they also have many backcrosses and recent generation hybrids (Harrison and Bogdanowicz 1997, Maroja et al. 2009a, Larson et al. 2013a, 2014). Because of this, we considered any collecting locations that were near the hybrid zone to be “sympatric”. We also assigned each collecting location to a geographic region (labeled in Fig. 1). These regions, identified using climatological data (Karl and Koss 1984), were as follows: central (CTR: IL, IN, KY, MO, OH, TN, WV); east north central (ENC: IA, MI, MN, WI); northeast (NE: CT, DE, ME, MD, MA, NH, NJ, NY, PA, RI, VT); northwest (NW: ID, OR, WA); south (SO: AR, KS, LA, MS, OK, TX); southeast (SE: AL, FL, GA, NC, SC, VA); southwest (SW: AZ, CO, NM, UT); west (WE: CA, NV); and west north central (WNC: MT, NE, ND, SD, WY).

In all cases, crickets were collected by hand and maintained in plastic containers with food (cat and rabbit food), water vials, and shelter prior to freezing. Most samples were collected as adults, but in some cases, crickets were collected as late instar nymphs. Nymphs were allowed to mature to the adult stage in the laboratory before freezing. Most collections were done in August–September, but some crickets were collected in late July or early October.

Morphological measurements.—We focused only on traits that were measured using the same methods across different studies. Crickets were measured for body size, as gauged by either body length, femur length, and/or pronotum width. Body length was measured from the vertical surface of the face to the tip of the abdomen, straightening the body when necessary. Pronotum width was measured at the widest part of the pronotum. Femur length was measured from the proximal to distal end of the hind femur. Female ovipositor length was measured from the point of attachment on the abdomen to the distal end of the ovipositor. Because ovipositor length varies isometrically with body size (Suppl. material 1: fig. S1), we also calculated relative ovipositor length as the length of the ovipositor divided by pronotum width or femur length, depending on sample availability. We obtained all measurements using Vernier calipers and recorded values to the nearest 0.1 mm.

For a subset of samples where tegmina were available (31 allopatric crickets and 437 sympatric crickets), we measured their color using a USB4000 spectrophotometer with an Ocean Optics PX-2 pulsed xenon lamp and SpectraSuite v2.0 software. We mounted a probe on a metal stand at a 90° angle 0.7 mm from the surface of the tegmina. For each male, we recorded and averaged spectral reflectance for three points near the center of the tegmina. We recorded spectral measurements as the percentage of reflected light relative to a Spectralon white standard, restricted our analyses to wavelengths of 300–700 nm, and used a segmental classification method to estimate brightness, chroma, and hue using CLR v1.1 (Montgomery 2008). We calculated total brightness (B) as R300–700, which is the summed reflectance from 300 nm to 700 nm. We also divided our reflectance data into four bins of 100 nm each, calculated the total brightness for each bin (Br=600–700, By=500–600, Bg=400–500, and Bb=300–400), and then calculated chroma: √(BrBg)²+(ByBb)² and hue: arctan[(ByBb)/B]/[(BrBg)/B].

Molecular markers.—A subset of the crickets in our dataset was previously genotyped for mitochondrial DNA haplotype (N = 1,132, Harrison et al. 1987, Harrison and Bogdanowicz 1997, Willett et al. 1997, Maroja et al. 2009a, Larson et al. 2013b) and/or 110 Single Nucleotide Polymorphisms (SNPs) from nuclear genes with elevated divergence between G. pennsylvanicus and G. firmus (N = 559, Larson et al. 2013a, 2014). Mitochondrial DNA haplotype was determined by sequencing cytochrome c oxidase I, the adjacent tRNA-Leu, and a portion of cytochrome c oxidase II (Harrison et al. 1987, Willett et al. 1997). SNPs were identified from transcriptomes of male accessory glands from two focal populations (Ithaca, NY and Guilford, CT; Andrés et al. 2013) were genotyped using the Sequenom MassARRAY platform (Larson et al. 2013a, 2014). We used these genotype data to recalculate the hybrid index while accounting for hemizygosity for male X-linked markers using the methods from Shastry et al. (2021). This was especially important because nearly half of these 110 SNPs are located on the X chromosome (Maroja et al. 2015, Gainey et al. 2018). We defined the hybrid index as the proportion of alleles that were inherited from G. firmus (hybrid index = 1; Guildford, CT (GUI); Tom’s River, NJ (TOM); and Parksley, MD (MET, a.k.a. PAR in Larson et al. 2013a, 2014) and G. pennsylvanicus (hybrid index = 0; Ithaca, NY (ITH); Scranton, PA (SCR); State College, PA (SCO)).

Analysis of morphological traits and molecular markers.—All analyses were conducted in R v4.1.2 (R Core Team 2020). To manipulate the data, we used the R packages dplyr v1.0.6 and tidyverse v1.3.1. To plot our data, we used the R packages ggplot2 v3.3.5 and ggpubr v0.4.0, and to make our maps, we used Maps v3.3.0. For statistical analyses, we used commands from the R packages MASS v7.3-54 and car v3.0-12. We used the R packages corrplot v0.92 and Hmisc v4.5-0 to determine environmental variable correlation. We used the R packages AICcmodavg v2.3-1 and MnMln 1.43.1 to rank models based on Akaike Information Criterion and test models.

To test for differences in morphological traits between species and regions, we used the Kruskal-Wallis test followed by a pairwise Wilcoxon rank sum test (PWRST) to determine differences between multiple groups. We chose these non-parametric tests because our dataset failed Levene’s test for homogeneity of variance. We quantified how well morphological traits could classify crickets using a linear discriminant analysis (LDA) on allopatric crickets. For all analyses, we present the unadjusted p-values and indicate in bold the values that were significant following FDR correction (Benjamini and Hochberg 1995).

Environmental predictors of species distributions.—We tested the relationships between phenotype and environmental variables that we predicted would be important in determining species range or local adaptation on two scales: 1) across species ranges and 2) at an intermediate scale in a well-characterized region of the hybrid zone (Connecticut). Across the species ranges, we used only allopatric crickets that were most clearly differentiated by morphology, and at the intermediate scale, we used both allopatric and sympatric crickets. We focused on the two phenotypes that best distinguished the two species and that were quantified in most of our samples: ovipositor length and pronotum width.

We identified 10 environmental variables that might be good predictors of species’ distributions based on the natural history of these species and prior studies of the field cricket hybrid zone (longitude, latitude, elevation, precipitation, minimum temperature, maximum temperature, human footprint, and three soil characteristics; see Larson et al. 2013b). Elevation, precipitation, and temperature data were collected from the PRISM Climate Group website (https://prism.oregonstate.edu/). Elevation was calculated using an 800-m digital elevation model of the continental United States. For each site, we collected precipitation variables and minimum and maximum temperatures for the year in which each cricket was collected. PRISM data were not available for sites in Canada. Soil data were collected from the USDA STATSGO2 soil survey (US sites: https://www.nrcs.usda.gov/wps/portal/nrcs/detail/soils/survey/geo/?cid=nrcs142p2_053629) and the Soil Landscapes of Canada database (Canada sites: https://sis.agr.gc.ca/cansis/nsdb/slc/v3.2/index.html). For a subset of sites in the northeastern United States, we used soil data from ISRIC SoilGrids (Poggio et al. 2021) due to the smaller spatial scale. These data were accessed and compiled using the R package soilDB v2.6.14. We used the following variables: average percent sand, average percent clay, and average percent organic matter. Due to the high intercorrelation of soil variables confirmed through correlation matrix, we excluded average soil percent silt from further analyses. We also obtained spatial data from the Last of the Wild Global Human Footprint dataset (version 3), consisting of anthropogenic impact measured by population density, land use, and transportation access at a 1-km resolution (Venter et al. 2016, 2018).

We used model selection tests that included these 10 environmental variables to find the combination of variables that best explains morphological variation. We ranked competing models using Akaike Information Criterion (AIC), and we reported the models with the highest goodness-of-fit.

Data accessibility.—All morphological data and collection site information, including GPS coordinates and environmental data and scripts, are published in Dryad (doi: 10.5061/dryad.jwstqjqdx).

Estimates of body size.—In total, our dataset comprised 4,221 crickets, with > 1,100 crickets per sex for each morphological trait measured, except for male tegmina color (Table 1). We first evaluated the relationship between three morphological traits that reflect overall body size in crickets: body length, femur length, and pronotum width. We found that body length measurements could vary depending on how crickets responded to being frozen in the lab or other factors such as number of eggs or last meal (see also Weissman and Gray 2019). Consequently, we chose to exclude body length measurements from our analyses but include them in our supplemental datasets. Male and female individuals of both G. pennsylvanicus and G. firmus had strong positive relationships between femur length and pronotum width (male G. pennsylvanicus: R² = 0.53, F_1,233 = 265, p < 2.2×10^-16 and male G. firmus: R² = 0.76, F_1,117 = 363.1, p < 2.2×10^-16, Suppl. material 1: fig. S1A; female G. pennsylvanicus: R² = 0.53, F_1,192 = 21, p < 2.2×10^-16 and female G. firmus: R² = 0.74, F_1,89 = 254.7, p < 2.2×10^-16, Suppl. material 1: fig. S1B). Therefore, we used pronotum width as our estimate for overall body size to maximize the number of individuals we could compare across datasets. In female individuals, pronotum width and ovipositor length were also positively related in both species (G. pennsylvanicus: R² = 0.44, F_1,214 = 165.7, p < 2.2×10^-16 and G. firmus: R² = 0.26, F_1,87 = 30.48, p = 3.44×10^-7, Suppl. material 1: fig. S1C). In comparisons with G. thinos, we used femur length to estimate body size to maximize the number of individuals in those comparisons.

Morphological differences between species.—There were significant differences among allopatric G. pennsylvanicus, G. firmus, G. thinos, and sympatric populations (e.g., G. firmus, G. pennsylvanicus, and hybrids) in male body size (Kruskal-Wallis, 𝝌² = 35.79, df = 3, p = 8.29×10^-8), female body size (Kruskal-Wallis, 𝝌² = 51.89, df = 3, p = 3.16×10^-11), female ovipositor length (Kruskal-Wallis, 𝝌² = 1277.2, df = 3, p < 2.2×10^-16), and relative ovipositor length (Kruskal-Wallis, 𝝌² = 82.10, df = 3, p < 2.2×10^-16) . When comparing allopatric G. pennsylvanicus and G. firmus, male pronotum (p = 2.1×10^-5, Fig. 2A), female pronotum (p = 1.4×10^-11, Fig. 2B), ovipositor length (p < 2.2×10^-16, Suppl. material 1: fig. S2A), and relative ovipositor length (p = 2.8×10^-16, Fig. 2C) were all significantly different. However, for each of these traits, there was still considerable overlap between allopatric species. Ovipositor length had the most striking differences between species (Suppl. material 1: fig. S2A), even when controlling for body size (Fig. 2C).

Fig. 2. 

Allopatric populations of G. firmus and G. pennsylvanicus differ in overall body size and ovipositor length. A. Male pronotum width by species; B. Female pronotum width by species; C. Relative ovipositor length (ovipositor length/pronotum width). Boxplots indicate the mean values of each trait, quartiles, the range of the data (whiskers), and outliers. Individual data points are overlaid as scatterplots. Letters indicate the significant differences among groups (PWRST with corrected p-values < 0.05).

For males, tegmina color alone classified most individuals from allopatric populations as either G. pennsylvanicus or G. firmus (LDA, misclassification rate of 3%). One of the 24 G. pennsylvanicus males was misclassified as G. firmus, and zero of the 7 G. firmus males were misclassified as G. pennsylvanicus. When looking at male body size alone, the misclassification rate was much higher at 23%, with 56 of the 268 G. pennsylvanicus males misclassified and 27 of the 90 G. firmus males misclassified. There was not enough overlap in body size and tegmina color data to perform these analyses using both variables. For females, body size and relative ovipositor length classified most individuals from allopatric populations as either G. pennsylvanicus or G. firmus (LDA, misclassification rate 12%). Fifteen of the 189 G. pennsylvanicus were misclassified as G. firmus. and 17 of the 90 G. firmus were misclassified as G. pennsylvanicus.

Crickets from areas near the hybrid zone, which we refer to as sympatric, had considerable overlap with those from allopatric populations. Sympatric crickets were not different from G. firmus for male body size, but they were, on average, larger than G. pennsylvanicus (G. pennsylvanicus: p = 6.0×10^-6, G. firmus: p = 0.16, Fig. 2A) but were still different from both allopatric species for female body size (G. pennsylvanicus: p = 9.4×10^-7, G. firmus: p = 0.00032, Fig. 2B), female ovipositor length (G. pennsylvanicus: p < 2.0×10^-16, G. firmus: p < 2.0×10^-16, Suppl. material 1: fig. S2A), and female relative ovipositor length (G. pennsylvanicus: p = 4.6×10^-8, G. firmus: p = 1.0×10^-9, Fig. 2C). This suggests that while these sympatric populations may have individuals that are more G. firmus-like or G. pennsylvanicus-like, they still have intermediate morphology compared to allopatric populations.

Intraspecific variation in key morphological traits.—We then tested how these traits varied across the different geographic regions of each species. We found differences among regions of G. pennsylvanicus for male pronotum (Kruskal-Wallis, 𝝌² = 56.11, df = 6, p = 2.76×10^-10), female pronotum (Kruskal-Wallis, 𝝌² = 63.44, df = 6, p = 8.9×10^-12), ovipositor length (Kruskal-Wallis, 𝝌² = 185.72, df = 6, p < 2.2×10^-16), and relative ovipositor length (Kruskal-Wallis, 𝝌² = 33.6, df = 6, p = 8.03×10^-6). Male and female G. pennsylvanicus were largest in the southern and midcentral US (SE, SO, SW, CTR, Fig. 3A, B), and they had the smallest body size in the northern west (WNC, NW). There were differences among regions in G. firmus pronotum width (Kruskal–Wallis, males, 𝝌² = 9.27, df = 2, p = 0.01; females, 𝝌² = 9.15, df = 2, p = 0.01), in ovipositor length (Kruskal–Wallis, 𝝌² = 78.65, df = 2, p < 2.2×10^-16), and relative ovipositor lengths (Kruskal-Wallis, 𝝌² = 54.49, df = 2, p = 1.47×10^-12). Male and female G. firmus were larger in the south than in the northeast, while G. firmus in the south were not significantly different from crickets in either the northeast or the southeast (Fig. 3A, B). In G. pennsylvanicus, ovipositor length varied by region. Eastern populations (NE, SE) had the shortest ovipositors, and those from the central US (CTR) had the longest ovipositors, although there was a very limited sample size for this region (Suppl. material 1: figs S2B, S3C). There was considerable variation in ovipositor length among G. firmus populations; southern G. firmus females had significantly shorter relative ovipositors than G. firmus in the southeast, who in turn had shorter relative ovipositors than G. firmus in the northeast (Suppl. material 1: figs S2C, S3C). However, G. firmus in the southeast had very similar absolute ovipositor lengths to northeastern G. firmus but had larger body sizes, whereas southern G. firmus simply had shorter ovipositors (Suppl. material 1: fig. S2C).

Fig. 3. 

Cricket body size and relative ovipositor length varies by geographic region. A. Male pronotum width by species and region; B. Female pronotum width by species and region; C. Relative ovipositor length by species and region. Boxplots indicate the mean values of each trait, quartiles, the range of the data (whiskers), and outliers. Individual data points are overlaid as scatterplots. Letters indicate the significant differences among groups within each species (PWRST with corrected p-values < 0.05), and exact p-values are presented in Suppl. material 1: tables S1, S2. See Fig. 1 for location information.

Recent work by Weissman and Gray (2019) documented cryptic variation in southern USA G. firmus, so we took a closer look at these populations, separating crickets collected in Florida from those collected in Texas. We also included the recently described closely related species G. thinos, which is sympatric with Texas G. firmus (Weissman and Gray 2019). We found that male (Kruskal–Wallis, 𝝌² = 29.26, df = 3, p = 1.98×10^-6) and female (Kruskal-Wallis, 𝝌² = 24.88, DF = 3, p = 1.63×10^-5) body size and ovipositor length (ovipositor length: Kruskal-Wallis, 𝝌² = 101.39, df = 3, p < 2.2×10^-16; relative ovipositor: Kruskal-Wallis, 𝝌² = 89.57, df = 3, p < 2.2×10^-16) differ among these groups (Fig. 4, Suppl. material 1: fig. S2). Compared to northeastern G. firmus, Florida G. firmus were much larger (Fig. 4A, B) but had only slightly larger ovipositor lengths (Suppl. material 1: fig. S2C), giving them shorter relative ovipositors (Fig. 4C). Texas G. firmus did not differ in overall body size from northeastern G. firmus but had even shorter relative ovipositor lengths (Fig. 4C, Suppl. material 1: fig. S2C). The magnitude of the morphological differences among Florida, Texas, and northeastern G. firmus is similar to that of the differences between G. firmus and the recently described G. thinos. Gray et al. (2020) found that G. firmus in Texas and Florida are genetically distinct groups, with Texas G. firmus sister to G. pennsylvanicus and Florida G. firmus sister to both G. pennsylvanicus and Texas G. firmus. Altogether, the morphological differences and phylogenetic relationships support the findings by Weissman and Gray (2019) that Texas G. firmus may be an undescribed cryptic species.

Fig. 4. 

Morphological variation in G. firmus consistent with proposed cryptic species. A. Male femur length; B. Female femur length; C. Relative ovipositor length (ovipositor length/femur length). There is considerable morphological variation among northeastern, Florida, and Texas G. firmus, which is similar to the magnitude of morphological divergence observed in the closely related species G. thinos. This combined with genetic divergence suggests there may be cryptic species in what is currently considered G. firmus. Boxplots indicate the mean values of each trait, quartiles, the range of the data (whiskers), and outliers. Individual data points are overlaid as scatterplots. Letters indicate the significant differences among groups (PWRST with corrected p-values < 0.05), and exact p-values are presented in Suppl. material 1: table S3.

Morphology in sympatric populations.—For the subset of crickets that were from the hybrid zone or nearby (sympatric populations) and were also genotyped with molecular markers, we looked at the relationship between admixture and morphological traits. We found that each trait had a similar transition from G. pennsylvanicus to G. firmus, with highly admixed individuals having intermediate phenotypes (Fig. 5). We found that male pronotum (R² = 0.19, F_1,279 = 63.35, p = 4.38×10^-14), male tegmina color (R² = 0.31, F_1,133 = 60.82, p = 1.62×10^-12), female pronotum (R² = 0.28, F_1,275 = 107.3, p < 2.2×10^-16), and relative ovipositor length (R² = 0.47, F_1,270 = 243.1, p < 2.2×10^-16) all had strong correlation with the hybrid index. Because the SNPs used to calculate the hybrid index are concentrated on the X chromosome (54 out of 110 (Maroja et al. 2015, Gainey et al. 2018)), females (XX) were more likely to be classified with an intermediate hybrid index than males (XO). Overall, morphological traits were also correlated with mtDNA haplotypes: crickets that had G. pennsylvanicus mtDNA tended to be smaller (males: Kruskal-Wallis, 𝝌² = 43.14, df = 1, p = 5.11×10^-11; females: Kruskal-Wallis, 𝝌² = 44.86, df = 1, p = 2.11×10^-11), darker (Kruskal-Wallis, 𝝌² = 33.75, df = 1, p = 6.27×10^-9) crickets with shorter ovipositors (Kruskal-Wallis, 𝝌² = 37.67, df = 1, p = 8.40×10^-10) (Fig. 6). We found that crickets with G. firmus ancestry at nuclear markers (hybrid index = 1) often had G. pennsylvanicus mtDNA haplotypes (Fig. 7), indicating asymmetric introgression of the mtDNA.

Fig. 5. 

Crickets with more hybrid background have intermediate morphological traits. The relationship between the hybrid index (an estimate of ancestry proportions, G. pennsylvanicus = 0 and G. firmus = 1) and A. Male pronotum width; B. Female pronotum width; C. Relative ovipositor length, and D. Male tegmina color.

Fig. 6. 

Morphological traits tended to correspond to mtDNA haplotypes. A. Male pronotum width; B. Female pronotum width; C. Relative ovipositor length; and D. Male tegmina color. Boxplots indicate the mean values of each trait, quartiles, the range of the data (whiskers), and outliers. Individual data points are overlaid as scatterplots.

Fig. 7. 

Mitochondrial DNA introgression is largely asymmetric. Crickets with G. firmus ancestry at nuclear markers (hybrid index = 1) often had G. pennsylvanicus mtDNA haplotypes.

Environmental predictors of morphology.—In allopatric populations throughout broad ranges, we found that latitude, elevation, average soil percent clay, and minimum and maximum temperatures created the best model for ovipositor length. Latitude, longitude, soil percent sand, and minimum temperature created the best model for pronotum width (Table 2). Average soil percent clay and higher minimum and maximum temperatures were positively associated with longer ovipositor lengths, and higher minimum temperatures were positively associated with larger body size, which are characteristics of G. firmus (Suppl. material 1: fig. S3). In the subset of Connecticut sympatric and allopatric populations, minimum and maximum temperatures, as well as soil percent organic matter, created the best model, with positive associations for all three variables and ovipositor length (Table 2, Suppl. material 1: fig. S3).

Cryptic diversity in a wide-ranging species.—The hybrid zone between the field crickets G. firmus and G. pennsylvanicus has been a model for understanding speciation (Harrison and Rand 1989, Harrison and Larson 2014). The field cricket hybrid zone stretches from the northeastern United States as far south as Virginia and likely farther into the southeast. Divergence in morphology, nuclear and mitochondrial DNA, and reproductive barriers have been carefully studied in several major regions of the hybrid zone (Harrison 1985, Rand and Harrison 1989, Ross and Harrison 2002, Maroja et al. 2009a, 2009b, Larson et al. 2012, 2014). Yet, even in this well-studied system, there is geographic diversity across the ranges of these species that complicates their relationships.

Our results confirm that allopatric populations of these two species, defined by genetic markers (Harrison and Arnold 1982, Willett et al. 1997, Broughton and Harrison 2003, Maroja et al. 2009a), can be largely differentiated by a combination of body size, male tegmina color, and female ovipositor length (Fig. 2). At the same time, there is regional variation in these traits within each species (Fig. 3). These differences may be due to local adaptation of life history traits such as egg diapause and development time (discussed below) or phenotypic plasticity. However, in some cases, they may also indicate cryptic diversity in field crickets.

In their revision of North American field crickets, Weissman and Gray (2019) proposed that there was cryptic diversity in the southern populations of G. firmus, particularly in Texas. Importantly, our phenotypic comparisons confirmed that Texas and Florida G. firmus are morphologically distinct from northeastern G. firmus (Fig. 4). In a recent nuclear phylogeny, Texas and Florida G. firmus-like crickets also formed distinct clusters within the larger G. pennsylvanicus group (Weissman and Gray 2019, Gray et al. 2020). Unfortunately, we do not have a phylogeny that includes genes from both Texas and Florida G. firmus and northeastern G. firmus, so the relationships among these groups are still unclear. However, the combination of distinct morphology and phylogenetic relationships suggests that at least one cryptic species of Gryllus exists, a situation that will not be resolved without further genotyping and/or evaluations of reproductive compatibility among these populations.

Intermediate phenotypes in hybrid zone crickets.—The morphological traits that best distinguish species in allopatry can also be used to distinguish these species in or near the hybrid zone. In this study, we took a conservative approach to defining allopatric and sympatric populations. Allopatric populations were those well outside of where the two species co-occur and are typically populations that have been genotyped with species-diagnostic markers. We found that in sympatry, crickets that were mostly G. firmus or mostly G. pennsylvanicus at nuclear markers (Larson et al. 2013a, 2014) had morphological traits that are also G. firmus-like or G. pennsylvanicus-like. Both male and female body size, male tegmina color, and relative ovipositor length had clinal variations from G. pennsylvanicus-like to G. firmus-like, with highly admixed individuals having intermediate phenotypes (Fig. 5). Male tegmina color stood out as having the fewest individuals with intermediate hybrid index values (Fig. 5D), but this is most likely because the SNPs used to calculate the hybrid index were predominately X-linked, so male XO crickets were rarely heterozygous at those SNPs and had overall lower hybrid indices (Larson et al. 2014, Maroja et al. 2015, Gainey et al. 2018).

The relationship between morphology and mitochondrial haplotype was less clear for populations near or in the hybrid zone. Crickets that were mostly G. firmus at the nuclear markers often had G. pennsylvanicus mtDNA (Fig. 7). This pattern fits with what we expect based on the one-way prezygotic incompatibility between G. firmus females and G. pennsylvanicus males (Harrison 1983, Maroja et al. 2009b, Larson et al. 2012). All F1 hybrids are produced from crosses with G. pennsylvanicus mothers; thus, G. pennsylvanicus mtDNA will be more likely to introgress into G. firmus. Even rare instances of hybridization might lead to mtDNA introgression, such as the mtDNA capture observed in many mammal species (Melo-Ferreira et al. 2005, Good et al. 2008).

Adaptations to soil type.—Ovipositor length is one of the most striking morphological differences between G. firmus and G. pennsylvanicus. Female crickets use their ovipositors to lay their eggs in the soil, and ovipositor length has been hypothesized to relate to the soil type and/or the depth of egg laying (Masaki 1979). The depth of egg laying may be a particularly critical life-history trait in G. pennsylvanicus and G. firmus because these species overwinter as eggs, as opposed to most field crickets that overwinter as early instar nymphs (Alexander 1968, Harrison and Bogdanowicz 1995). For eggs to be viable, they must withstand low winter temperatures and freeze/thaw cycles (Ross and Harrison 2006). Throughout its range, G. firmus is most often found on sandy coastal soils (Harrison and Arnold 1982, Weissman and Gray 2019) and tends to have a longer ovipositor than G. pennsylvanicus (Fig. 3, Suppl. material 1: fig. S2). This may be an adaptation to laying eggs deeper in sandy substrates in response to intermittent rainfall and the risk of eggs drying out (Walker 1980). In some parts of the hybrid zone, such as Connecticut, the association with different soil types is striking. The two species have been found on micro habitat patches of loam (G. pennsylvanicus) and sandy (G. firmus) soils in Connecticut (Harrison 1986, Harrison and Rand 1989, Rand and Harrison 1989), and interactions between the two species occur across these habitat patch boundaries on a scale of only hundreds of meters (Ross and Harrison 2002, Larson et al. 2014).

Despite what appears to be strong habitat associations, the relationship between soil type and ovipositor length is complicated. Ovipositor length does not necessarily determine egg-laying depth; instead, females may wield long ovipositors at different angles (Réale and Roff 2002). It is also not clear exactly how the association between ovipositor length and soil type is maintained. Females of both species prefer to lay eggs in loamy soil, and there is no difference in overwintering egg viability in different soil types (Ross and Harrison 2006). Finally, these associations are clearly established only in a small part of the species’ ranges, i.e., Connecticut (Rand and Harrison 1989, Ross and Harrison 2002, Larson et al. 2013b). Even where soil associations appear to be the strongest, the transition from sandy to loamy soils is more gradual and less distinct than we might expect based on the patchiness of G. firmus and G. pennsylvanicus populations (Ross and Harrison 2002, Larson et al. 2014). Here we find that both across the broad ranges of these species and at an intermediate scale in the Connecticut hybrid zone, there is no strong association between ovipositor length and sandy soils. In fact, we tend to see crickets with longer ovipositors on clay soils (Table 2). This might be due to the different methods used to quantify soil type (soil survey data versus on-site soil sampling), but altogether, this suggests that habitat associations in these species are variable and should be investigated further.

Body size, climate, and life cycle.—In insects, seasonality and the length of the growing season are critical to the rate of development and adult body size (Masaki 1961, Tauber and Tauber 1981). This is particularly true for hemimetabolous insects, which often go through many nymphal stages and have long development times before reaching their full size and sexual maturity (Kivelä et al. 2011). Insects at higher latitudes have shorter growing seasons and, as a result, may develop more quickly or reach an overall smaller body size (Masaki 1967, Parsons and Joern 2014). This pattern of smaller body sizes at higher latitudes is sometimes referred to as the converse of Bergman’s rule, which states that individuals have larger body sizes in colder climes (Masaki 1967, Mousseau 1997). We see this pattern most clearly in G. pennsylvanicus, where we found that populations with the smallest body sizes tended to be farther north (WNC and NW, Fig. 2). Indeed, we found that crickets at higher latitudes had, on average, smaller body sizes and that there was a significant relationship between body size and latitude (Table 2).

We may not expect a direct relationship between body size and latitude if the length of the growing season allows for multiple generations per year. Insects can shift from continuous development in the south to univoltine (one generation per year) in the north (Masaki 1961, 1967). As a result, there may be regions where body size is smaller than expected based on latitude to accommodate multiple generations per year. We did not find this pattern in our results, but we may not have had the resolution of latitudinal samples to see a sawtooth pattern in body size. However, there is some evidence that development time in G. firmus varies with latitude. In Virginia, G. firmus emerge earlier in the season than G. pennsylvanicus, leading to temporal isolation in that part of the hybrid zone, but in Connecticut, the two emerge simultaneously (Harrison 1985). In Florida, G. firmus is reported to have multiple generations per year (Walker, personal observation; reported in Weissman and Gray 2019), where throughout its range, it otherwise appears to have a single generation per year (Walker 1980). Notably, despite having many generations per year, Florida G. firmus are considerably larger than northern populations. It is unclear whether there is a continuous shift in life cycle across the range of G. firmus or if Florida G. firmus have a distinct life history from other G. firmus.

In studies of speciation and to understand the effects of local selection, it is critical to quantify morphological and genetic variations across the geographic range of widespread species. The field cricket hybrid zone is an example of how important the larger geographic context can be. In some regions of the field cricket hybrid zone, G. pennsylvanicus and G. firmus have a patchy distribution, and G. firmus crickets are found on sandy soils (Rand and Harrison 1989, Ross and Harrison 2002). However, the strong soil association breaks down in other regions of the hybrid zone (Larson et al. 2013b) and across their geographic range, suggesting that the soil association may be a result of local adaptation or colonization history (Hauffe and Searle 1993, Gompert et al. 2010). Our results provide a foundation for future geographically expansive studies that compare genetic divergence and the role of specific traits in reproductive barriers to better understand local adaptation and speciation in this system. More broadly, this is an example of how critical it is to move studies of speciation beyond the comparison of a few focal populations. Geographically expansive studies of phenotypic and genetic divergence will also be important for understanding how species distributions and hybrid zones shift over time and in a changing climate (Britch and Cain 2001, Taylor et al. 2015).

ARB, CJ, DBW, DAG, CLR, LSM and ELL collected the data; ARB, CJ, and ELL combined the datasets; AG obtained the environmental data and advised on analyses. ARB conducted all analyses. ARB and ELL wrote the manuscript with contributions from all authors.

Rick Harrison collected many of the samples in this dataset. He was a wonderful mentor and friend to many of the authors and an inspiration to all of us. We thank the University of Denver Ecology and Evolution group, especially Shannon Murphy and Jonathan Velotta, for their feedback on this manuscript. This work was supported by grants to ARB from the Explorer’s Club, The Orthopterists’ Society, and the Society of Systematic Biologists and National Science Foundation grants to ELL (DEB 2012041) and LSM (DEB 2012060).