Research Article Journal of Orthoptera Research 2023, 32(2): 189-200 Geographic variation in phenotypic divergence between two hybridizing field cricket species Amy R. Byerty!, CLARA JENCK!, ALEXANDER R. B. Goetz!, David B. WEIssmAN2, Davip A. Gray, CHARLES L. Ross’, LUANA S. MAROJA®, ERICA L. LARSON! 1 Department of Biological Sciences, University of Denver, Denver, CO 80208, USA. 2 Department of Entomology, California Academy of Sciences, Golden Gate Park, San Francisco, CA 94118, USA. 3 Department of Biology, California State University, Northridge, CA 91330, USA. 4 School of Natural Science, Hampshire College, Amherst, MA 01002, USA. 5 Department of Biology, Williams College, Williamstown, MA 01267, USA. Corresponding authors: Amy R. Byerly (amy.byerly@du.edu); Erica L. Larson (erica.larson@du.edu) Academic editor: Hojun Song | Received 22 July 2022 | Accepted 23 March 2023 | Published 25 September 2023 https://zoobank. org/D3283705-EA19-4CB8-9ED6-E113287B2149 Citation: Byerly AR, Jenck C, Goetz ARB, Weissman DB, Gray DA, Ross CL, Maroja LS, Larson EL (2023) Geographic variation in phenotypic divergence between two hybridizing field cricket species. Journal of Orthoptera Research 32(2): 189-200. https://doi.org/10.3897/jor.32.90713 Abstract Patterns of morphological divergence across species’ ranges can pro- vide insight into local adaptation and speciation. In this study, we com- pared phenotypic divergence among 4,221 crickets from 337 populations of two closely related species of field cricket, Gryllus firmus and G. pennsyl- vanicus, and their hybrids. We found that these species differ across their geographic range in key morphological traits, such as body size and ovi- positor length, and we directly compared phenotype with genotype for a subset of crickets to demonstrate nuclear genetic introgression, pheno- typic intermediacy of hybrids, and essentially unidirectional mitochon- drial 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 impor- tance of studying reproductive barriers in more than one or two transects of a hybrid zone. Keywords Gryllus, hybrid zone, local adaptation, morphology, Orthoptera, speciation Introduction Phenotypic divergence can provide insight into evolutionary processes acting across different scales of biological organiza- tion. 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, Gonzalez 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 be- tween 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 Wojcik 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 in- dividuals mate and produce offspring (Barton and Hewitt 1985, Harrison 1990). In the resulting hybrid zone, as the different spe- cies 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 (Hol- lander et al. 2018, Sottas et al. 2018). By comparing the phenotypic variation between both conspecific allopatric and sympatric popu- lations and between heterospecific populations, it becomes possi- ble to examine the potential causes of phenotypic evolution, speci- ation, 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 be- tween 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 Copyright Amy R. Byerly et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. JOURNAL OF ORTHOPTERA RESEARCH 2023, 32(2) 190 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., Connecti- cut), but whether these morphological traits are consistently dif- ferent 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 distin- guish 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. Like- wise, body size may vary with climate and latitude. This paper pre- sents the first geographically comprehensive comparison of G. fir- mus and G. pennsylvanicus by combining published and unpub- lished 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 A 50 S > 40 ios} 23] @ 50 @ pennsylvanicus @ 100 =3=6©—_ firmus @ 150 @ thinos @ sympatric sd @ 200 ame Longitude A.R. BYERLY, C. JENCK, A.R.B. GOETZ, D.B. WEISSMAN, D.A. GRAY, C.L. ROSS, L.S. MAROJA AND E.L. LARSON 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. Materials and methods Cricket collections. —We compiled a dataset of 4,221 crickets, the majority being G. pennsylvanicus but also G. firmus and their hy- brids, from 337 collecting localities (Fig. 1). Crickets were sampled throughout the United States and Canada, with the largest collec- tions 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 mor- phological data, we included data from Ross and Harrison (2002), Larson et al. (2013), and Weissman and Gray (2019), with the lat- ter being the most geographically widespread dataset. We also in- cluded 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. We categorized each collecting location as allopatric or sympat- ric 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 lo- cations 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, 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. JOURNAL OF ORTHOPTERA RESEARCH 2023, 32(2) A.R. BYERLY, C. JENCK, A.R.B. GOETZ, D.B. WEISSMAN, D.A. GRAY, C.L. ROSS, L.S. MAROJA AND E.L. LARSON 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 laborato- ty before freezing. Most collections were done in August-Septem- ber, 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, straight- ening 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 abdo- men 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 ovi- positor 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 al- lopatric 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 mount- ed 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 re- flectance 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 es- timate 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: V(BrBg)?+(ByBb)? and hue: arctan|(ByBb)/B]/[(BrBg)/B]. Molecular markers.—A subset of the crickets in our dataset was pre- viously 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. fir- mus (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 identi- fied 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 191 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 AlCcmodavg v2.3-1 and MnMIn 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 pair- wise 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 vari- ance. 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 re- lationships 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 al- lopatric crickets that were most clearly differentiated by morphol- ogy, 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 tempera- ture, 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 conti- nental United States. For each site, we collected precipitation vari- ables 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 STATS- GO2 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 soiIDB v2.6.14. We used the following variables: average percent sand, average percent clay, and average percent organic matter. Due to the high JOURNAL OF ORTHOPTERA RESEARCH 2023, 32(2) 192 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 transporta- tion access at a 1-km resolution (Venter et al. 2016, 2018). We used model selection tests that included these 10 environ- mental 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 in- formation, including GPS coordinates and environmental data and scripts, are published in Dryad (doi:10.5061/dryad.jwstqjqdx). Results Estimates of body size.—In total, our dataset comprised 4,221 crick- ets, with > 1,100 crickets per sex for each morphological trait meas- ured, 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 de- pending 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 supple- mental datasets. Male and female individuals of both G. pennsylvan- icus and G. firmus had strong positive relationships between femur length and pronotum width (male G. pennsylvanicus: R? = 0.53, F, ,,, = 265, p < 2.2x10'° and male G. firmus: R’ = 0.76, F, ,, = 363.1, p < 2.2x10°'°, Suppl. material 1: fig. S1A; female G. pennsylvanicus: R? = 0.53, F, ,,, = 21, p < 2.2x10"° and female G. firmus: R’ = 0.74, F, ,, = 254.7, p < 2.2x10"°, Suppl. material 1: fig. S1B). Therefore, we used pronotum width as our estimate for overall body size to maxi- mize the number of individuals we could compare across datasets. A.R. BYERLY, C. JENCK, A.R.B. GOETZ, D.B. WEISSMAN, D.A. GRAY, C.L. ROSS, L.S. MAROJA AND E.L. LARSON In female individuals, pronotum width and ovipositor length were also positively related in both species (G. pennsylvanicus: R* = 0.44, F, 4. = 165.7, p < 2.2x10"° and G. firmus: R* = 0.26, F, ,, = 30.48, p = 3.44x10’, 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. thi- nos, and sympatric populations (e.g., G. firmus, G. pennsylvani- cus, and hybrids) in male body size (Kruskal-Wallis, 7? = 35.79, df = 3, p = 8.29x 10°), female body size (Kruskal-Wallis, 7? = 51.89, df = 3, p = 3.16x10"), female ovipositor length (Kruskal-Wallis, X = 1277.2, df = 3, p < 2.2x10"°), and relative ovipositor length (Kruskal-Wallis, 7? = 82.10, df = 3, p < 2.2x10"° ). When com- paring allopatric G. pennsylvanicus and G. firmus, male pronotum (p = 2.1x10°, Fig. 2A), female pronotum (p = 1.4x10", Fig. 2B), ovipositor length (p < 2.2x107°, Suppl. material 1: fig. S2A), and relative ovipositor length (p = 2.8x10°'°, Fig. 2C) were all signifi- cantly 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). 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. Table 1. Summary of sample sizes for morphological measurements by sex, population type, and region (See Fig. 1 for location information). Pronotum Width Femur Length Ovipositor Ovipositor Ovipositor Tegmina Females Males Females Males Length Pronotum Ratio Femur Ratio Color Totals 1203 1263 134 1213 4047 1174 1110 469 CTR 4 5 4 5 12 4 4 - allopatric 4 5 4 5 4 4 4 - sympatric - - - 8 - - - NE 993 1010 871 849 3739 969 851 449 allopatric 111 167 85 132 1480 108 82 23 sympatric 882 843 786 717 2259 861 769 426 NW 26 17 27 15 27 26 27 - allopatric 26 LZ 27 15 27 26 27, - SE 66 65 77 60 111 62 74 20 allopatric 66 65 Oe 60 89 62 74 8 sympatric - - - - 22 - - 18. SO 40 69 66 171 70 40 66 - allopatric 40 69 65 171 69 40 65 - sympatric - - 1 1 - 1 - SW 29 41 29 52 29 29 29 - allopatric 29 41 29 52 29 25 25 - WNC 45 56 60 61 59 44 59 - allopatric 45 56 60 61 59 44 ony - JOURNAL OF ORTHOPTERA RESEARCH 2023, 32(2) A.R. BYERLY, C. JENCK, A.R.B. GOETZ, D.B. WEISSMAN, D.A. GRAY, C.L. ROSS, L.S. MAROJA AND E.L. LARSON > Ww € a a b E E® =& Piet’ ) Female pronotum width (mm) (o>) 4 NE CTR SE SO SW Cc 4 e a de ab cd pl ex D 5 4 S 3 = fe) ° 2 oS > (eo) ® 3y%e = ai Is & er 2 NE CTR SE SO SW pennsylvanicus A.R. BYERLY, C. JENCK, A.R.B. GOETZ, D.B. WEISSMAN, D.A. GRAY, C.L. ROSS, L.S. MAROJA AND E.L. LARSON WNC NW NE SE SO be. be a b c NE SE SO WNC NW firmus 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. morphological traits were also correlated with mtDNA haplotypes: crickets that had G. pennsylvanicus mtDNA tended to be smaller (males: Kruskal-Wallis, x? = 43.14, df= 1, p = 5.1110"; females: Kruskal-Wallis, x* = 44.86, df= 1, p = 2.11x10-"), darker (Kruskal- Wallis, x? = 33.75, df = 1, p = 6.27x10°) crickets with shorter ovi- positors (Kruskal-Wallis, x* = 37.67, df= 1, p = 8.40x 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. Environmental predictors of morphology.—In allopatric populations throughout broad ranges, we found that latitude, elevation, average soil percent clay, and minimum and maximum temperatures cre- ated 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 char- acteristics of G. firmus (Suppl. material 1: fig. $3). In the subset of Connecticut sympatric and allopatric populations, minimum and maximum temperatures, as well as soil percent organic mat- ter, created the best model, with positive associations for all three variables and ovipositor length (Table 2, Suppl. material 1: fig. $3). JOURNAL OF ORTHOPTERA RESEARCH 2023, 32(2) A.R. BYERLY, C. JENCK, A.R.B. GOETZ, D.B. WEISSMAN, D.A. GRAY, C.L. ROSS, L.S. MAROJA AND E.L. LARSON 195 A B Cc 18 18 = ‘= = = 16 — 16 = = = & D . =) fe f 4 514 = r : a 5 iF - o 12 2 12 g Ss = 3 i r 10 : 10 1.0 NE Florida Texas thinos NE Florida Texas thinos NE Florida Texas thinos firmus firmus firmus 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. Box- plots 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. A B Hybrid index 8 8 1.00 0.75 7 = 7 a mi € e 0.50 & e E e 3 ~ * < ‘i e6 s % s 0.25 g 6 e = 6 : @ 6 € 0.00 E Bh z = e g ®@ 2 e 5 ° a ° e =n © é ro @ ov) s 5 rm e 4 4 3 3 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 Hybrid Index Hybrid Index Cc 4.0 D e 2.5 S 3.5 2 = ® jo) e = c 4 E x) e aa — = 4) is} e Soo @@ 3 9 2 on Gs 2 > © - e (oy E gv ® ee 7 e 3 ° ee e % pa 2 -2.5 r) c S e ® -5.0 2.0 -7.5 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 Hybrid Index Hybrid Index 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. JOURNAL OF ORTHOPTERA RESEARCH 2023, 32(2) A B KKKK pane. KKKK € = © 8 5 5 ‘= 5 35 6 5 6 x) ec € fe) S ot roe o4 a4 o Ss 5 LL firmus pennsylvanicus firmus pennsylvanicus A.R. BYERLY, C. JENCK, A.R.B. GOETZ, D.B. WEISSMAN, D.A. GRAY, C.L. ROSS, L.S. MAROJA AND E.L. LARSON Cc D KKKK KKKK 4 al a, D 5 S = ‘Gnd fo) (©) = (40) (7p) Cc a? = Ss D0 fo) 2 oO $2 2 is = ke a) or pennsylvanicus firmus firmus pennsylvanicus Fig. 6. Morphological traits tended to correspond to mtDNA haplotypes. A. Male pronotum width; B. Female pronotum width; C. Rela- tive 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. Table 2. Results of linear regression and AIC to test the relationship between environmental variables and morphological traits in female crickets of both species. | indicates variables where values are based on the year the samples were collected. A. Ovipositor length Df Sum of Sq RSS AIC Coefficient St. Error t-value p-value (Intercept) - 887.86 321.96 16.213 0.132 122.417 < 2.00E-16 Latitude 1 9.283 897.14 322.33 0.545 0.358 1523 0.129 Precipitation!’ 1 1.054 886.81 323.69 - - - - Longitude 1 0.389 887.47 323.86 - - - - Human Footprint 1 0.132 887.73 323.93 - - - - Avg Soil % Sand 1 0.015 887.85 323.96 - - - - Avg Soil % Organic Matter 1 0.004 887.86 323,96 - - - - Elevation 1 26.035 913.90 326.55 -0.638 0.250 -2.551 0.011 Avg Soil % Clay 1 26.562 914.42 326.68 0.365 0.142 2.577 0.011 Minimum Temperature’ 1 29.311 SUASLT, 327.36 -1.244 0.459 -2.707 0.007 Maximum Temperature’ 1 124.629 1012.49 349.91 2.281 0.409 5.582 6.89E-08 B. Pronotum width (Intercept) - 56.539 25309 5.83503 0.036 162.636 < 2.00E-16 Maximum Temperature’ 1 0.381 36.158 -253.69 - - - - Precipitation! 1 0.176 36.363 20 2;¢3 - - - - Human Footprint 1 0.147 36.392 -252.59 - - - - Elevation 1 0.103 36.436 -252.38 - - - - Avg Soil % Clay 1 0.088 36.451 -252.31 - - - - Avg Soil % Organic Matter 1 0.013 36.526 -251.96 - - - - Minimum Temperature’ 1 2.964 39.503 -242.56 -0.233 0.064 -3.670 3.27E-04 Avg Soil % Sand 1 3953 40.492 -238.33 -0.180 0.043 -4.238 3.73E-05 Longitude 1 4.618 41.157 -235.55 -0.187 0.041 -4.580 9.07E-06 Latitude 1 12.890 49.429 -204.23 -0.536 0.070 -7.652 1.53E-12 Discussion 2009a), can be largely differentiated by a combination of body 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 far- ther into the southeast. Divergence in morphology, nuclear and mitochondrial DNA, and reproductive barriers have been careful- ly 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. size, male tegmina color, and female ovipositor length (Fig. 2). At the same time, there is regional variation in these traits with- in each species (Fig. 3). These differences may be due to local adaptation of life history traits such as egg diapause and devel- opment time (discussed below) or phenotypic plasticity. How- ever, 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. Impor- tantly, 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 larg- er G. pennsylvanicus group (Weissman and Gray 2019, Gray et al. 2020). Unfortunately, we do not have a phylogeny that includes JOURNAL OF ORTHOPTERA RESEARCH 2023, 32(2) A.R. BYERLY, C. JENCK, A.R.B. GOETZ, D.B. WEISSMAN, D.A. GRAY, C.L. ROSS, L.S. MAROJA AND E.L. LARSON 200 - ah ol Oo mtDNA haplotype | Northern firmus ) Southern firmus Mi Northern pennsylvanicus I Southern pennsylvanicus Number of crickets ° oO 50 - 0.50 0.75 1.00 Hybrid index 0.25 0.00 Fig. 7. Mitochondrial DNA introgression is largely asymmetric. Crickets with G. firmus ancestry at nuclear markers (hybrid in- dex = 1) often had G. pennsylvanicus mtDNA haplotypes. genes from both Texas and Florida G. firmus and northeastern G. firmus, so the relationships among these groups are still un- clear. However, the combination of distinct morphology and phy- logenetic relationships suggests that at least one cryptic species of Gryllus exists, a situation that will not be resolved without fur- ther genotyping and/or evaluations of reproductive compatibility among these populations. Intermediate phenotypes in hybrid zone crickets. —The morphologi- cal 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 allopat- ric 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 oviposi- tor length had clinal variations from G. pennsylvanicus-like to G. firmus-like, with highly admixed individuals having inter- mediate 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 mark- ers often had G. pennsylvanicus mtDNA (Fig. 7). This pattern fits with what we expect based on the one-way prezygotic in- compatibility 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). 197 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 re- late 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 over- winter as eggs, as opposed to most field crickets that overwinter as early instar nymphs (Alexander 1968, Harrison and Bogdano- wicz 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 compli- cated. Ovipositor length does not necessarily determine egg-laying depth; instead, females may wield long ovipositors at different an- gles (Réale and Roff 2002). It is also not clear exactly how the as- sociation 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 Harri- son 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 (Kivela 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 JOURNAL OF ORTHOPTERA RESEARCH 2023, 32(2) 198 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 mul- tiple generations per year. Insects can shift from continuous de- velopment 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 ac- commodate 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. Howev- er, there is some evidence that development time in G. firmus var- ies 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 simultane- ously (Harrison 1985). In Florida, G. firmus is reported to have multiple generations per year (Walker, personal observation; re- ported 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. Conclusions In studies of speciation and to understand the effects of local selection, it is critical to quantify morphological and genetic varia- tions across the geographic range of widespread species. The field cricket hybrid zone is an example of how important the larger geo- graphic 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 Har- rison 1989, Ross and Harrison 2002). However, the strong soil as- sociation 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 un- derstanding how species distributions and hybrid zones shift over time and in a changing climate (Britch and Cain 2001, Taylor et al. 2015). Author Contributions ARB, CJ, DBW, DAG, CLR, LSM and ELL collected the data; ARB, CJ, and ELL combined the datasets; AG obtained the envi- ronmental data and advised on analyses. ARB conducted all analy- ses. ARB and ELL wrote the manuscript with contributions from all authors. A.R. BYERLY, C. JENCK, A.R.B. GOETZ, D.B. WEISSMAN, D.A. GRAY, C.L. ROSS, L.S. MAROJA AND E.L. LARSON Acknowledgements 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 Den- ver 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). References Alexander RD (1968) Life cycle origins, speciation, and related phenom- ena in crickets. The Quarterly Review of Biology 43: 1-41. https://doi. org/10.1086/405628 Andrés JA, Larson EL, Bogdanowicz SM, Harrison RG (2013) Patterns of transcriptome divergence in the male accessory gland of two closely related species of field crickets. Genetics 193: 501-513. https://doi. org/10.1534/genetics.112.142299 Barton NH, Hewitt GM (1985) Analysis of hybrid zones. Annual Review of Ecology and Systematics 16: 113-148. https://doi.org/10.1146/an- nurev.es. 16.110185.000553 Benjamini Y, Hochberg Y (1995) Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society. Series B, Statistical Methodology 57: 289-300. https://doi.org/10.1111/j.2517-6161.1995.tb02031.x Britch SC, Cain ML (2001) Spatio-temporal dynamics of the Allonemo- bius fasciatus— A. socius mosaic hybrid zone: a 14-year perspective. Molecular Ecology 10: 627-638. https://doi.org/10.1046/j.1365- 294x.2001.01215.x Broughton RE, Harrison RG (2003) Nuclear gene genealogies reveal his- torical, demographic and selective factors associated with speciation in field crickets. Genetics 163: 1389-1401. https://doi.org/10.1093/ genetics/163.4.1389 Gainey DP, Kim JY, Maroja LS (2018) Mapping reduced introgression loci to the X chromosome of the hybridizing field crickets, Gryllus firmus and G. pennsylvanicus. PLoS ONE 13: e0208498. https://doi. org/10.1371/journal.pone.0208498 Gavrilets S, Arnqvist G, Friberg U (2001) The evolution of female mate choice by sexual conflict. Proceedings of the Royal Society of London. Series B: Biological Sciences 268: 531-539. https://doi.org/10.1098/ rspb.2000.1382 Gompert Z, Lucas LK, Fordyce JA, Forister ML, Nice CC (2010) Sec- ondary contact between Lycaeides idas and L. melissa in the Rocky Mountains: extensive admixture and a patchy hybrid zone. Mo- lecular Ecology 19: 3171-3192. https://doi.org/10.1111/j.1365- 294X.2010.04727.x Gonzalez C, Ornelas JE Gutiérrez-Rodriguez C (2011) Selection and geographic isolation influence hummingbird speciation: genetic, acoustic and morphological divergence in the wedge-tailed sabrew- ing (Campylopterus curvipennis). BMC Evolutionary Biology 11: 38. https://doi.org/10.1186/1471-2148-11-38 Good JM, Hird S, Reid N, Demboski JR, Steppan SJ, Martin-Nims TR, Sul- livan J (2008) Ancient hybridization and mitochondrial capture be- tween two species of chipmunks. Molecular Ecology 17: 1313-1327. https://doi.org/10.1111/j.1365-294X.2007.03640.x Gray DA, Weissman DB, Cole JA, Lemmon EM, Lemmon AR (2020) Mul- tilocus phylogeny of Gryllus field crickets (Orthoptera: Gryllidae: Gryllinae) utilizing anchored hybrid enrichment. Zootaxa 4750: 328-348. https://doi.org/10.11646/zootaxa.4750.3.2 Harrison RG (1983) Barriers to gene exchange between closely relat- ed cricket species. I. Laboratory hybridization studies. Evolution; International Journal of Organic Evolution 37: 245-251. https://doi. org/10.2307/2408333 JOURNAL OF ORTHOPTERA RESEARCH 2023, 32(2) A.R. BYERLY, C. JENCK, A.R.B. GOETZ, D.B. WEISSMAN, D.A. GRAY, C.L. ROSS, L.S. MAROJA AND E.L. LARSON Harrison RG (1985) Barriers to gene exchange between closely related cricket species. II. Life cycle variation and temporal isolation. Evo- lution; International Journal of Organic Evolution 39: 244-259. https://doi.org/10.2307/2408360 Harrison RG (1986) Pattern and process in a narrow hybrid zone. Heredity 56: 337-349. https://doi.org/10.1038/hdy.1986.55 Harrison RG (1990) Hybrid zones: windows on evolutionary process. Oxford Surveys in Evolutionary Biology 7: 69-128. Harrison RG, Arnold J (1982) A narrow hybrid zone between closely related cricket species. Evolution 36: 535. https://doi.org/10.2307/2408099 Harrison RG, Rand DM (1989) Mosaic hybrid zones and the nature of spe- cies boundaries. Speciation and its Consequences: 111-133. Harrison RG, Bogdanowicz SM (1995) Mitochondrial DNA phylogeny of North American field crickets: perspectives on the evolution of life cy- cles, songs, and habitat associations. Journal of Evolutionary Biology 8: 209-232. https://doi.org/10.1046/j.1420-9101.1995.8020209.x Harrison RG, Bogdanowicz SM (1997) Patterns of variation and link- age disequilibrium in a field cricket hybrid zone. Evolution 51: 493. https://doi.org/10.2307/2411122 Harrison RG, Larson EL (2014) Hybridization, introgression, and the na- ture of species boundaries. The Journal of Heredity 105: 795-809. https://doi.org/10.1093/jhered/esu033 Harrison RG, Larson EL (2016) Heterogeneous genome divergence, dif- ferential introgression, and the origin and structure of hybrid zones. Molecular Ecology 25: 2454-2466. https://doi.org/10.1111/mec.13582 Harrison RG, Rand DM, Wheeler WC (1987) Mitochondrial DNA varia- tion in field crickets across a narrow hybrid zone. Molecular Biology and Evolution 4: 144-144. Hauffe HC, Searle JB (1993) Extreme karyotypic variation in a Mus muscu- lus domesticus hybrid zone: the tobacco mouse story revisited. Evo- lution; International Journal of Organic Evolution 47: 1374-1395. https://doi.org/10.2307/2410154 Hollander J, Montano-Rend6n M, Bianco G, Yang X, Westram AM, Duvaux L, Reid DG, Butlin RK (2018) Are assortative mating and genital di- vergence driven by reinforcement? Evolution Letters 2: 557-566. https://doi.org/10.1002/ev13.85 Jenck CS, Lehto WR, Ketterman BT, Sloan LE, Sexton AN, Tinghitella RM (2020) Phenotypic divergence among threespine stickleback that differ in nuptial coloration. Ecology and Evolution 10: 2900-2916. https://doi.org/10.1002/ece3.6105 Jiménez RA, Ornelas JF (2015) Historical and current introgression in a Mesoamerican hummingbird species complex: a biogeographic per- spective. Peer) 4: e1556. https://doi.org/10.7717/peerj.1556 Karl T, Koss WJ (1984) Regional and national monthly, seasonal, and an- nual temperature weighted by area, 1895-1983. https://repository. library.noaa.gov/view/noaa/ 10238 [April 12, 2021] Kivela SM, Valimaki P, Carrasco D, Maenpéa MI, Oksanen J (2011) Lati- tudinal insect body size clines revisited: a critical evaluation of the saw-tooth model. The Journal of Animal Ecology 80: 1184-1195. https://doi.org/10.1111/j.1365-2656.2011.01864.x Larson EL, Hume GL, Andrés JA, Harrison RG (2012) Post-mating prezy- gotic barriers to gene exchange between hybridizing field crickets. Journal of Evolutionary Biology 25: 174-186. https://doi.org/10.1111/ j.1420-9101.2011.02415.x Larson EL, Andrés JA, Bogdanowicz SM, Harrison RG (2013a) Differential introgression in a mosaic hybrid zone reveals candidate barrier genes. Evolution; International Journal of Organic Evolution 67: 3653- 3661. https://doi.org/10.1111/evo.12205 Larson EL, Guilherme Becker C, Bondra ER, Harrison RG (2013b) Structure of a mosaic hybrid zone between the field crickets Gryllus firmus and G. pennsylvanicus. Ecology and Evolution 3: 985-1002. https://doi. org/10.1002/ece3.514 Larson EL, White TA, Ross CL, Harrison RG (2014) Gene flow and the maintenance of species boundaries. Molecular Ecology 23: 1668- 1678. https://doi.org/10.1111/mec.12601 Maroja LS, Andrés JA, Harrison RG (2009a) Genealogical discordance and patterns of introgression and selection across a cricket hybrid 199 zone. Evolution; International Journal of Organic Evolution 63: 2999-3015. https://doi.org/10.1111/j.1558-5646.2009.00767.x Maroja LS, Andrés JA, Walters JR, Harrison RG (2009b) Multiple barriers to gene exchange in a field cricket hybrid zone. Biological Journal of the Linnean Society 97: 390-402. https://doi.org/10.1111/j.1095- 8312.2009.01201.x Maroja LS, Larson EL, Bogdanowicz SM, Harrison RG (2015) Genes with restricted introgression in a field cricket (Gryllus firmus/Gryllus penn- sylvanicus) hybrid zone are concentrated on the X chromosome and a single autosome. G3 Genes|Genomes|Genetics 5: 2219-2227. https://doi.org/10.1534/g3.115.021246 Masaki S (1961) Geographic variation of diapause in insects. Bulletin of the Faculty of Agriculture, Hirosaki University 7: 66-98. Masaki S (1967) Geographic variation and climatic adaptation in a field cricket (Orthoptera: Gryllidae). Evolution; International Journal of Organic Evolution 21: 725-741. https://doi.org/10.2307/2406770 Masaki S (1979) Climatic adaptation and species status in the lawn ground cricket. Oecologia 43: 207-219. https://doi.org/10.1007/BF00344771 Melo-Ferreira J, Boursot P, Suchentrunk FE, Ferrand N, Alves PC (2005) Invasion from the cold past: extensive introgression of mountain hare (Lepus timidus) mitochondrial DNA into three other hare spe- cies in northern Iberia. Molecular Ecology 14: 2459-2464. https:// doi.org/10.1111/j.1365-294X.2005.02599.x Moran PA, Hunt J, Mitchell C, Ritchie MG, Bailey NW (2020) Sexual se- lection and population divergence III: Interspecific and intraspecific variation in mating signals. Journal of Evolutionary Biology 33: 990- 1005. https://doi.org/10.1111/jeb.13631 Mousseau TA (1997) Ecotherms follow the converse to Bergmann’s rule. Evolution; International Journal of Organic Evolution 51: 630-632. https://doi.org/10.2307/2411138 Oneal E, Knowles LL (2013) Ecological selection as the cause and sexual differentiation as the consequence of species divergence? Proceedings of the Royal Society B, Biological Sciences 280: 20122236. https://doi. org/10.1098/rspb.2012.2236 Parsons SMA, Joern A (2014) Life history traits associated with body size covary along a latitudinal gradient in a generalist grasshopper. Oeco- logia 174: 379-391. https://doi.org/10.1007/s00442-013-2785-6 Poggio L, de Sousa LM, Batjes NH, Heuvelink GBM, Kempen B, Ribeiro E, Rossiter D (2021) SoilGrids 2.0: producing soil information for the globe with quantified spatial uncertainty. Soil 7: 217-240. https://doi. org/10.5194/soil-7-217-2021 Polly PD, Wojcik JM (2019) Geometric morphometric tests for phenotypic divergence between chromosomal races. In: Shrews, Chromosomes and Speciation. Cambridge University Press, 336-364. https://doi. org/10.1017/9780511895531.011 Rand DM, Harrison RG (1989) Ecological genetics of a mosaic hybrid zone: mitochondrial, nuclear and reproductive differentiation of crickets by soil type. Evolution; International Journal of Organic Evo- lution 43: 432-449. https://doi.org/10.2307/2409218 R Core Team (2020) R: A language and environment for statistical com- puting. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/ Réale D, Roff DA (2002) Quantitative genetics of oviposition behaviour and interactions among oviposition traits in the sand cricket. Animal Behaviour 64: 397-406. https://doi.org/10.1006/anbe.2002.3084 Ross CL, Harrison RG (2002) A fine-scale spatial analysis of the mosaic hybrid zone between Gryllus firmus and Gryllus pennsylvanicus. Evo- lution; International Journal of Organic Evolution 56: 2296-2312. https://doi.org/10.1111/j.0014-3820.2002.tb00153.x Ross CL, Harrison RG (2006) Viability selection on overwintering eggs in a field cricket mosaic hybrid zone. Oikos 115: 53-68. https://doi. org/10.1111/j.2006.0030-1299.15054.x Runemark A, Hansson B, Pafilis P, Valakos ED, Svensson EI (2010) Is- land biology and morphological divergence of the Skyros wall liz- ard Podarcis gaigeae: a combined role for local selection and ge- netic drift on color morph frequency divergence? BMC Evolutionary Biology 10: 269. JOURNAL OF ORTHOPTERA RESEARCH 2023, 32(2) 200 Shastry V, Adams PE, Lindtke D, Mandeville EG, Parchman TL, Gompert Z, Buerkle CA (2021) Model-based genotype and ancestry estimation for potential hybrids with mixed-ploidy. Molecular Ecology Resources 21: 1434-1451. https://doi.org/10.1111/1755-0998.13330 Shaw KL, Mullen SP (2011) Genes versus phenotypes in the study of spe- ciation. Genetica 139: 649-661. https://doi.org/10.1007/s10709-011- 9562-4 Skoglund S, Siwertsson A, Amundsen P-A, Knudsen R (2015) Morphologi- cal divergence between three Arctic charr morphs - the significance of the deep-water environment. Ecology and Evolution 5: 3114-3129. https://doi.org/10.1002/ece3.1573 Sottas C, Reif J, Kuczynski L, Reifova R (2018) Interspecific competition promotes habitat and morphological divergence in a secondary con- tact zone between two hybridizing songbirds. Journal of Evolutionary Biology 31: 914-923. https://doi.org/10.1111/jeb.13275 Tauber CA, Tauber MJ (1981) Insect seasonal cycles. Annual Review of Ecology and Systematics 12: 281-308. https://doi.org/10.1146/an- nurev.es. 12.110181.001433 Taylor SA, Larson EL, Harrison RG (2015) Hybrid zones: windows on cli- mate change. Trends in Ecology & Evolution 30: 398-406. https://doi. org/10.1016/j.tree.2015.04.010 Uyeda JC, Arnold SJ, Hohenlohe PA, Mead LS (2009) Drift promotes speciation by sexual selection. Evolution; International Journal of Organic Evolution 63: 583-594. https://doi.org/10.1111/j.1558- 5646.2008.00589.x Venter O, Sanderson EW, Magrach A, Allan JR, Beher J, Jones KR, Possin- gham HP, Laurance WE, Wood P, Fekete BM, Levy MA, Watson JEM (2016) Sixteen years of change in the global terrestrial human foot- print and implications for biodiversity conservation. Nature Commu- nications 7: 12558. https://doi.org/10.1038/ncomms12558 Venter O, Sanderson EW, Magrach A, Allan JR, Beher J, Jones KR, Poss- ingham HP, Laurance WE, Wood P, Fekete BM, Levy MA, Watson JE (2018) Last of the Wild Project, Version 3 (LWP-3): 2009 Human Footprint, 2018 Release. Palisades, New York: NASA Socioeconomic Data and Applications Center (SEDAC). https://doi.org/10.7927/ H46TOJQ4 [Accessed October 18, 2021] Walker TJ (1980) Mixed oviposition in individual females of Gryllus firmus: Graded proportions of fast-developing and diapause eggs. Oecologia 47: 291-298. https://doi.org/10.1007/BF00398519 Wang G-H, Li H, Zhao H-W, Zhang W-K (2017) Detecting climatically driven phylogenetic and morphological divergence among spruce (Picea) species worldwide. Biogeosciences 14: 2307-2319. https:// doi.org/10.5194/bg-14-2307-2017 Weissman DB, Gray DA (2019) Crickets of the genus Gryllus in the United States (Orthoptera: Gryllidae: Gryllinae). Zootaxa 4705: 1- 277. htt- ps://doi.org/10.11646/zootaxa.4705.1.1 Willett CS, Ford MJ, Harrison RG (1997) Inferences about the origin of a field cricket hybrid zone from a mitochondrial DNA phylogeny. He- redity 79: 484-494. https://doi.org/10.1038/hdy.1997.188 Wolf JBW, Harrod C, Brunner S, Salazar S, Trillmich F, Tautz D (2008) Tracing early stages of species differentiation: ecological, morpho- logical and genetic divergence of Galapagos sea lion populations. BMC Evolutionary Biology 8: 150. https://doi.org/10.1186/1471- 2148-8-150 A.R. BYERLY, C. JENCK, A.R.B. GOETZ, D.B. WEISSMAN, D.A. GRAY, C.L. ROSS, L.S. MAROJA AND E.L. LARSON Supplementary material 1 Author: Amy R. Byerly, Clara Jenck, Alexander R. B. Goetz, Da- vid B. Weissman, David A. Gray, Charles L. Ross, Luana S. Maroja, Erica L. Larson Data type: docx Explanation note: table $1. P-values for PWRST posthoc contrasts of allopatric G. pennsylvanicus populations by region (See Fig. 3). P- values marked with *** lost significance after correction. table S2. P-values for PWRST posthoc contrasts of allopatric G. firmus populations by region (See Fig. 3). P-values marked with *** lost significance after correction. table $3. P-values for PWRST post- hoc contrasts of G. thinos populations and G. firmus populations in the Northeast, Florida, and Texas. See Fig. 4. figure S1. Rela- tionship among phenotypic characteristics for allopatric popu- lations. (a) Male Gryllus firmus (yellow) and G. pennsylvanicus (teal) linear models show positive relationships in femur length and pronotum width (G. firmus: R2 = 0.76, F1,117 = 363.1, p-val- ue < 2.2e-16; G. pennsylvanicus: R2 = 0.53, F1,233 = 265.0, p-value < 2.2e-16). Female linear models show positive relationships in both (b) femur length and pronotum width (G. firmus: R2 = 0.74, F1,89 = 254.7,p-value < 2.2e-16; G. pennsylvanicus: R2 = 0.53, F1,192 = 217.0, p-value < 2.2e-16) and (c) ovipositor length and pronotum width (G. firmus: R2 = 0.26, F1,87 = 30.48, p-value = 3.44e-07; G. pennsylvanicus: R2 = 0.44, F1,214 = 165.7, p-value < 2.2e-16). figure S2. Ovipositor length differences between species and among populations of each species. A. Ovipositor length dif- ferences between G. firmus, G. pennsylvanicus and sympatric popu- lations (G. firmus vs G. pennsylvanicus p < 2.0E-16; G. firmus vs sympatric p <2.0E-16; G. pennsylvanicus vs sympatric p <2.0E-16). B. Ovipositor length differences among populations of G. penn- sylvanicus. Posthoc p-values are presented in Suppl. material 1: table S1. C. Ovipositor length differences among populations of G. firmus and G. thinos. Posthoc p-values are presented in Suppl. material 1: table $3. Boxplots indicate the mean values of each trait, quartiles and the range of the data (whiskers). Individual data points are overlaid as scatterplots. Letters indicate the signifi- cant differences among groups (PWRST with corrected p-values < 0.05). figure S3. Scatterplots of significant AIC environmental variables. For all female allopatric populations: ovipositor length vs. latitude (a.), elevation (b.), average soil percent clay (c.), mini- mum temperature (d.), and maximum temperature (e.). For all female allopatric populations: pronotum width vs. latitude (f.), longitude (g.), average soil % sand (h.), and minimum tempera- ture (i.). For female allopatric and sympatric Connecticut popula- tions: ovipositor length vs. minimum temperature (j.), maximum temperature (k.), and average soil % organic matter (I.). Copyright notice: This dataset is made available under the Open Database License __(http://opendatacommons.org/licenses/ odbl/1.0/). The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited. Link: https://doi.org/10.3897/jor.32.90713.suppl1 JOURNAL OF ORTHOPTERA RESEARCH 2023, 32(2)