Lack of assortative mating might explain reduced phenotypic differentiation where two grasshopper species meet.

Hybridization is an evolutionary process with wide-ranging potential outcomes, from providing populations with important genetic variation for adaptation to being a substantial fitness cost leading to extinction. Here, we focussed on putative hybridization between two morphologically distinct species of New Zealand grasshopper. We collected Phaulacridium marginale and Phaulacridium otagoense specimens from a region where mitochondrial introgression had been detected and where their habitat has been modified by introduced mammals eating the natural vegetation and by the colonization of many non-native plant species. In contrast to observations in the 1970s, our sampling of wild pairs of grasshoppers in copula provided no evidence of assortative mating with respect to species. Geometric morphometrics on pronotum shape of individuals from areas of sympatry detected phenotypically intermediate specimens (putative hybrids), and the distribution of phenotypes in most areas of sympatry was found to be unimodal. These results suggest that hybridization associated with anthropogenic habitat changes has led to these closely related species forming a hybrid swarm, with random mating. Without evidence of hybrid disadvantage, we suggest a novel hybrid lineage might eventually result from the merging of these two species.

INTRODUCTION

When related species occur in the same region, they are expected to be reproductively isolated, or their sympatry is unlikely to be stable (Coyne, 2001). Choosing the wrong mate will be selected against if the resulting offspring have low fitness (potentially zero if hybrids are infertile). Positive assortative mating is the nonrandom pair formation in which individuals with similar genotypes mate with one another and is widely accepted as a primary factor preventing hybridization and gene flow between species (Vaux et al., 2016). However, lineages that meet only occasionally may experience little selection pressure on mate choice. If populations are not reproductively isolated, hybridization will allow interspecific gene flow (Harrison & Larson, 2014) which is now considered by most biologists to be a fairly common and important evolutionary process (Mallet, 2005). Historical hybridization is revealed by introgressed mtDNA and cpDNA across species boundaries (cytonuclear discordance; Sloan et al., 2017). By increasing genetic variance, hybridization can contribute directly to diversification, adaptation to new environments (Lewontin & Birch, 1966), or even speciation (Abbott et al., 2013; Mallet, 2007). Nevertheless, hybridization is also a potentially rapid route to extinction when interbreeding results in a population that contains only individuals that are a mix of the original genetic lineages (Rhymer & Simberloff, 1996; Todesco et al., 2016).

Secondary contact of related species provides natural testing of the costs of imperfect assortative mating (Butlin & Ritchie, 2013). Sometimes narrow and stable hybrid zones emerge where hybrid individuals have low fitness, but parental species lack assortative mating (Barton & Gale, 1993). In other situations, hybrid swarms develop where there is neither assortative mating nor hybrid disadvantage, resulting in the expansion of the hybrid population and loss of pure parentals (Anderson, 1948; Harrison & Bogdanwicz, 1997; Langton‐Myers et al., 2019). Populations that contain a mixture of hybrids, backcrosses and later generations often exhibit a range of genotypic and phenotypic forms. The formation of hybrid swarms appears to be particularly influenced by human activity, for example among weedy species within disturbed habitats (Hasselman et al., 2014; Lowe & Abbott, 2015). More generally, rapid environmental change associated with anthropogenic climate and landscape modification appears to be responsible for many instances of contemporary hybridization, associated with shifting ranges (Larson et al., 2019; Mallet et al., 2011). It is even possible that novel genetic combinations generated within hybrid swarms could facilitate evolution within colonizing populations as they adapt to new environments (Morgan‐Richards et al., 2004; Seehausen, 2004).

In New Zealand, recent anthropogenic habitat modification has resulted in some examples of inter‐species hybridisation (Morgan‐Richards et al., 2009), and detailed studies of grasshopper taxa have been especially revealing. Secondary contact and introgression of species within the genera Sigaus and Phaulacridium are associated with vegetation changes driven by agriculture (Dowle et al., 2014; Sivyer et al., 2018). Observations in the 1970s suggested that two New Zealand lowland grasshopper species were ecologically isolated in adjacent microhabitats that differed in aridity and vegetation (Westerman & Ritchie, 1984). The two Phaulacridium species were distinguished and formally described based on morphological differences in body size, angle of lateral carinae of pronotum, hind femur dimensions, and for males the frontal ridge between vertex and median ocellus on the head. Data from more than 120 adult specimens showed phenotypic clusters (Mallet, 1995) consistent with two species (Westerman & Ritchie, 1984). The parapatric distributions of Phaulacridium marginale and P. otagoense appeared to be stable over at least 12 years and interspecific mating was never observed in the wild nor in captivity (Westerman & Ritchie, 1984). Subsequently however, evidence for ongoing range expansion of one species (P. marginale, Figure 1) and putative interspecific hybridisation emerged (Goldberg et al., 2015; Sivyer et al., 2018), suggesting imperfect assortative mating with respect to species. The widespread replacement of native New Zealand forest with pastures appears to have facilitated the rapid expansion of the range of P. marginale to occupy the new open‐habitat. In contrast, due to the local semi‐arid mesoclimate much of the range of P. otagoense was open habitat prior to human arrival in New Zealand but is now grazed by sheep and rabbits promoting soil erosion and colonization by exotic plant species (Walker, 2000). Land modification including dam construction, irrigation of pastures and weed invasion is gradually transforming the limited natural open‐habitat in which P. otagoense exists (e.g. Duncan et al., 2001; Mark et al., 2009; Meurk, 2008; Walker, 2000), and potentially changing the interactions of grasshoppers previously shown to have distinct ecological requirements (Westerman & Ritchie, 1984).

FIGURE 1

Recorded spatial distribution of two species of New Zealand lowland grasshoppers, Phaulacridium marginale (red) and P. otagoense (blue) from Sivyer et al. (2018). An adult female specimen of each of these colour‐polymorphic species is shown

Here, we examined the extent to which the process of assortative mating operates to limit gene flow between these two endemic grasshopper species that were ecologically isolated from one another when studied 46 years ago (Westerman & Ritchie, 1984). We used the morphological traits first used to describe the two species in a clustering analysis to determine whether specimens fell into two groups as expected of two species. Gaussian mixture model‐based clustering analyses allowed us to determine that our sample was partitioned, and thus, we could assign each specimen to a morphological species. Despite finding evidence of hybridization, the key morphological traits used for initial species identification (hind femur length and shape and angle of pronotum edge) resulted in high assignment probabilities for the majority of our specimens. We identified locations where the two species met, and looked for evidence of positive assortative mating in wild populations from adult grasshoppers caught in copula. We compared the number of pairs we captured that comprised two individuals of the same species (conspecific) and the number of pairs involving one P. marginale and one P. otagoense (interspecific pair) with the number of each pair expected from a model of random mating with respect to species. We postulated that if assortative mating predominates (as suggested by Westerman & Ritchie, 1984), gene flow would be minimal, and the two species would remain morphologically differentiated and display a bimodal frequency of phenotypes. Alternatively, if prezygotic barriers break down, hybridisation would be possible in the absence of post‐zygotic barriers. Hybridization is expected to initially yield some individuals of intermediate phenotype, and many generations of hybridization would result in the loss of differentiation and a unimodal phenotype distribution. Therefore, we examined genetics and phenotype in more detail to look for evidence of gene flow. At one location, we documented the mitochondrial haplotype of males and females that had chosen to mate with one another. For all specimens, we used a geometric morphometric approach to characterize and quantify the fine‐scale phenotypic variation of pronotum shape. Geometric morphometrics allowed us to remove the confounding effect of size variation that could be influenced by local environment, thus enabling population comparisons to identify the degree of mixing of the two species. By combining evidence of current behaviour in the wild (conspecific or interspecific mate choice) with genetic and phenotypic data, we were able to consider both the process of mate choice and its outcome.

MATERIALS AND METHODS

Grasshopper sampling and identification

Shorthorn grasshoppers in the genus Phaulacridium (Key, 1992) occupy low elevation grassland habitats that were naturally scarce in New Zealand (Koot et al., 2020; Sivyer et al., 2018). Two endemic species differ in size, shape and geographic distribution, and have distinct genome sizes (18.15 pg cf 20.88 pg; Westerman & Ritchie, 1984). Phaulacridium marginale (Walker, 1870) is widespread throughout both major New Zealand islands and many near‐shore islands, whereas Phaulacridium otagoense (Westerman & Ritchie, 1984) is restricted to the semi‐arid central and southern South Island (Figures 1 and 2). In late summer (March 2018 and 2019), we searched areas where we expected to find the species together (Table S1; Sivyer et al., 2018). These grasshoppers are flightless so to avoid bias from the scale‐of‐choice effect (Rolán‐Alvarez et al., 2015) each location sample consists of grasshoppers collected from within 200m of each other. We hand‐collected adult males and females as individuals or as mating pairs and preserved these by freezing and storage in 95% ethanol. In addition, we remeasured a subset of grasshopper specimens collected in 2012–2015 (Table S1; Sivyer et al., 2018).

FIGURE 2

Collection of two lowland New Zealand grasshopper species for this study focussed on the region of overlap in their distributions in southern South Island. The number of Phaulacridium marginale (red) and P. otagoense (blue) specimens in population samples was determined using model‐based analysis of traditional morphology (Sivyer et al., 2018). Circle size is proportional to sample size at each location (see Table S2 for full details of location codes and species composition)

To identify each specimen, we applied a formalization of morphological characters used in the diagnostic key for these grasshopper species (Westerman & Ritchie, 1984) with digital measures of hind femur length, hind femur length/width ratio and pronotum angle at the anterior pronotum sulcus (Appendix S1; Figure S1). In combination with a parameterised Gaussian hierarchical clustering algorithm, these metrics give efficient diagnosis to identify two phenotypic clusters we recognized as species (Sivyer et al., 2018). Gaussian mixture model‐based clustering used the Mclust v5.4.2 package (Fraley & Raftery, 2003) in the R programming environment (R Core Team), analysing males and females separately as sexual dimorphism renders the phenotypic difference between sexes greater than species differences. Grasshoppers from 22 locations were included in this analysis (n = 535; Table S2). The optimal model of variance and number of clusters in the data were selected based on Bayesian information criteria (BIC), using the value of the maximized log likelihood with a penalty for the number of parameters in the model. The higher the BIC score, the lower the median classification uncertainty, the better the model fits the data (Appendix S1).

Because the size of an adult grasshopper will be determined by an interaction between genotype and the environment, femur length is likely to be a plastic trait; however, the two species differ in size substantially more than individuals of the same species differ among locations (Westerman & Ritchie, 1984). Although we classified specimens using the analysis of the full dataset, to study assortative mating we focussed on locations where the two species are sympatric, removing the environment as a source of size variation when comparing individuals within a location. We collected data on the composition of mating pairs and on the distribution of phenotypes only where both species were identified (see below). Evidence of assortative mating with respect to species could potentially be influenced by our ability to correctly identify species (which in turn might be obscured by the presence of hybrids). We used assignment probabilities from the Gaussian mixture model‐based clustering as an indication of potential error in our species identification. High assignment probability (1 – 0.95) could only be achieved if the variation within our sample contained two Gaussian distributions as expected of two morphological species. Low assignment probability to either species/cluster would suggest that some specimens had traits that were not unequivocally either P. marginale or P. otagoense.

Assortative mating

Because interspecific mating had not previously been observed in the field or in captivity (Westerman & Ritchie, 1984), we housed the two species together to determine whether these two species could mate with each other in captivity. We collected grasshoppers from two locations where only a single species occurs (Phaulacridium otagoense from Quartz Reef, Otago and P. marginale from Turitea, North Island). Males of P. otagoense were held with female P. marginale and vice versa. This captive no‐choice arrangement used eight individuals of each sex of each species and ended when we observed interspecific mating.

To test for positive assortative mating with respect to species in the wild, we focussed on mating pairs captured in copula at five locations where both species were identified using morphology (total of 84 pairs; Table 1). We estimated a location‐specific species ratio from our five local samples to generate expected values under a random model for each location. We compared the number of each pair type observed (marginale/marginale, marginale/otagoense, otagoense/otagoense) with the number of each pair type expected under the null hypothesis of random mating with respect to species. Because our sample size in each pair category was relatively small (max = 29; Table 1), we pooled conspecific pairs as a unique category. To compare observed with expected values, we applied a binomial exact test implemented with the binom.test function from the native R package with correction for sample size.

TABLE 1

Lack of assortative mating with respect to species (morphologically identified) in the New Zealand grasshoppers Phaulacridium marginale ( ) and P. otagoense ( ) at five sympatric locations where mating pairs were observed in the wild

Location	Species					Mating pairs (Expected / Observed)			Exact p
	Individuals		In copula		Total	Conspecific	Interspecific	Total
	m	o	m	o
AxG Graveyard Gully Rd	22	6	3	3	34	1.83 / 2	1.17 / 1	3	1.0
AxL Little Valley Rd	3	8	1	27	39	11.43 / 13	2.58 / 1	14	0.49
AxM Marshall Rd	15	1	59	1	76	28.46 / 29	1.54 / 1	30	1.0
TeE Lake Tekapo East	4	29	4	38	75	16.99 / 17	4.00 / 4	21	1.0
TeW Lake Tekapo West	12	8	12	20	52	8.05 / 8	7.95 / 8	16	1.0

Individuals and mating pairs were collected to test for assortative mating. The expected proportion of conspecific pairs under a random model was computed as 2pq (with p and q as the local species proportion) and compared with the observed number of conspecific pairs with an exact binomial test (Exact p).

MtDNA haplotyping

We sought evidence of the interspecific mtDNA capture expected of hybridisation between the two species. Distinct mtDNA lineages are associated with each species, so we examined haplotypes of mating Phaulacridium pairs collected where both species occur. Sympatry was observed beside Lake Tekapo (Figure 2), and a total of 25 pairs in copula were caught for this analysis. We predicted that assortative mating would result in grasshopper pairs in which males and females had the same mitochondrial lineage. We extracted genomic DNA and sequenced a 755 bp fragment of the mitochondrial DNA gene cytochrome oxidase subunit I (COI) as previously described (Sivyer et al., 2018). Evolutionary relationships were inferred with a minimum spanning network (Bandelt et al., 1999) implemented in PopArt (Leigh & Bryant, 2015).

Phenotypic variation using pronotum shape

If positive assortative mating with respect to species predominates, we expected the two grasshopper species to be morphologically differentiated, and a sample containing both species would produce a bimodal frequency of phenotypes. To quantify phenotypic variation within and between population samples of the two species, we quantified the shape of their pronotum as this structure is not susceptible to arbitrary changes during preservation (Friedrich et al., 2014). The shape of the posterior margin of the pronotum provides diagnostic differences that distinguish many of the New Zealand grasshopper species and has previously been shown to be amenable to geometric morphometric approaches (Carmelet‐Rescan et al., 2021; Dowle et al., 2014; Ober & Connolly, 2015). We applied a geometric morphometric analysis to digital images of the pronotum of 535 grasshopper specimens allowing us to remove size variation, which although important for species diagnosis of Phaulacridium grasshoppers (Westerman & Ritchie, 1984) might be susceptible to environmental influence. Twelve geometric landmarks were identified and digitized around the perimeter of the dorsal surface of the pronotum, and measurement error was estimated by repeat mounting, photography and landmark selection of four grasshoppers ten times (Appendix S2; Figure S2). Procrustes analysis of variance was performed to compare variation in landmark location among iterations (Fruciano, 2016). A Procrustes fit aligned by principal axes was performed in MORPHOJ (Klingenberg, 2011). Nonshape variation was therefore mathematically removed as position, orientation and size superimposes landmark configuration using least‐squares estimates for translation and rotation. Linearly uncorrelated variables were obtained from a principal component analysis (PCA) across all individuals (n = 535) and all landmarks (n = 12) in MORPHOJ. Statistically significant principal components (PCs) were identified using the broken‐stick test on eigenvalues, implemented in the R package vegan 2.5‐2 (Oksanen & Blanchet, 2019).

In each population sample, we expected a unimodal Gaussian distribution of variation on the PC1 axis where a single species was present, and a bimodal distribution where two reproductively segregated species occurred, so we tested for deviation from a unimodal Gaussian distribution using d’Agostino omnibus normality test. We applied this test to our nine largest population samples (n > 20; Table S2). We visualized trait variation using stacked histograms of species assignment based on morphometric analysis, and density distributions for all individuals per site. Data underlying all our analyses are available in the Dryad Digital Repository (https://datadryad.org/stash/share/P2STPRfZafMtGovKx6GUwXmlZ9YIB9Lgwxk3SEPq1X8.)

RESULTS

Sampling and taxonomic identification

Grasshoppers were collected in lowland grassland habitats at 22 locations mostly in the region of sympatry, and data gathered from 535 individuals (217 males and 318 females), including 218 specimens caught in copula (109 pairs; Table S2). Despite potential hybridisation, Gaussian mixture modelling with three traditional morphological traits distinguished two clusters for females and three clusters for males. The best fitting model for females was ellipsoidal, equal volume, shape and orientation with two clusters (EEE2; BIC = −2,550.92) corresponding to the two grasshopper species. For males, the optimal model had three clusters when all specimens were analysed (K = 3, diagonal, equal volume and shape (EEI3) model; BIC = −1,548.211). The two‐cluster model for males (EEE2) was only slightly less well fitting (BIC = −1,561.006; Figure S3) and was used for species identification of male specimens. Morphometric assignment of specimens to two morphological species (clusters) without a priori identification revealed that nine of our population samples had only P. marginale individuals, six had only P. otagoense and seven contained both species (Figure 2). The geographic range of Phaulacridium otagoense is nested inside that of P. marginale (Figure 1). The majority of assignment probabilities were more than 0.95 (410/535) but for 9.53% of our sample species assignment probability was less than 0.8 (see Dryad Digital Repository https://datadryad.org/stash/share/P2STPRfZafMtGovKx6GUwXmlZ9YIB9Lgwxk3SEPq1X8).

That interspecific mating can occur was confirmed by males and females of each species accepting the other as mates in captivity. In both species, males are smaller than females but P. otagoense are smaller than P. marginale so male P. marginale are a substantial physical burden for female P. otagoense (Figure S4).

In the wild, fifteen interspecific mating pairs of grasshopper were identified among 109 in copula samples (Table 1, Table S2). In ten, the male was P. marginale and five the male was P. otagoense. We used the species frequency at each site to generate expected values for interspecific pairs under a random mating hypothesis. None of our five sympatric population samples where grasshoppers were caught in copula differed significantly from expectations of random mating with respect to species (Table 1). Assignment probability of species identification for the specimens collected in copula from locations where both species occur was <0.8 for 15.12% of our sample. To determine whether this level of uncertainty would prevent us detecting positive assortative mating with respect to species, we analysed our data again, after rescoring ~15.12% of each sample to more closely resemble assortative mating (e.g. reclassifying approximately 15% of pairs as conspecific). In each sample, the new distribution of mating pairs (conspecific or interspecific) did not differ from expectations of random assorting. It would require 30% of our grasshopper pairs to be misclassified, for detection of potential assortative mating to fail due to difficulty assigning species.

mtDNA haplotypes

We obtained COI sequences from 95 of our new grasshopper specimens, including representatives of populations not previously sequenced (e.g. Hawkesburn Road, Table S1), and combined with existing data (Sivyer et al., 2018). Analysis of the resulting 755bp alignment of 243 sequences revealed five main haplotype clusters including the identification of a new lineage (V) from Phaulacridium otagoense collected in the southwest of this species’ range (Figure S5; Bannockburn Sluicings (DB), Hawkesburn Road (DH)). Lineage I is associated with P. marginale throughout its range in both South and North Island (Goldberg et al., 2015). Four distinct, geographically limited mtDNA lineages (II, III, IV and V) are now known within P. otagoense (Figure S5).

Of 25 random Phaulacridium pairs collected in copula at two locations beside Lake Tekapo (east + west), only six (24%) comprised male and female with the same mtDNA lineage (Table 2). All other pairs were combinations of lineage I (P. marginale) and lineage IV (P. otagoense). Only one of 25 pairs comprised male and female with both P. otagoense morphology and mtDNA (Table 2).

TABLE 2

Lack of assortative mating has resulted in mitochondrial introgression between two grasshopper species beside Lake Tekapo, New Zealand

Morphotype of pairs	mtDNA lineage of pairs				Total
	o♀ o♂	o♀ m♂	m♀ o♂	m♀ m♂
o♀ o♂	1	10	2	3	16
o♀ m♂		1	3	1	5
m♀ o♂			3		3
m♀ m♂	1				1
Total	2	11	8	4	25

Combinations of morphotype and mtDNA lineage encountered among 50 Phaulacridium otagoense (o) and P. marginale (m) grasshoppers collected from the wild in copula. We have combined collections from the east (TeE) and west (TeW) shore of Lake Tekapo.

Phenotypic shape variation

Interspecific copulation could lead to the presence of hybrids in the wild, and we observed intermediate sized individuals in sympatry. Grasshopper females are larger than males but size may also differ due to environmental effects; therefore, we chose to examine pronotum shape of adult grasshoppers in detail. To allow us to combine male and female data and focus on shape variation, we used a geometric analysis to eliminate isometric size variation by superimposition via Procrustes transformation, while capturing fine‐scale shape variation among individuals (Bookstein, 1991; Rohlf & Marcus, 1993). Our estimates of measurement error associated with photographing and digitizing the images were assessed by the calculation of an interclass correlation coefficient that expresses the repeatability of this analysis. The closer the ratio is to one, the smaller the measurement error variation is compared with the biological variation between individuals. Our interclass correlation coefficient (r = .9773) suggests that less than 4% of the variation observed among individuals is due to measurement error.

The first two PCA axes encompassed statistically significant pronotum shape variation (PC1 = 50.85%, PC2 = 27.12%; Figure S6) and a naïve Gaussian mixture model algorithm (Mclust v5.4.2) revealed four pronotum shape clusters within our dataset (EEE4 model; BIC = 3,655.958) spanning variation in the two species and two sexes. The first principal component of variation distinguished the two species and PC2 distinguished the two sexes, with greater shape variation detected within P. otagoense than within P. marginale (Figure S6). Overall, P. otagoense had a wider pronotum than P. marginale and had more acute external angles on the lateral carinae (Figure S6). We used a stacked barplot of the assignment probabilities (membership) of each specimen to the pronotum shape clusters (Figure 3, Figure S6) and found 60 grasshoppers (11% of sample) displayed pronotum shape with low assignment probability to alternative clusters (between 0.2 and 0.8). Most of these came from locations where traditional morphometrics indicated the two species were sympatric (Figure 2). To examine further the relative abundance of these putative products of hybridisation, and therefore stringency of barriers to reproduction, we analysed the prevalence of alternative pronotum shapes in each population sample containing ≥20 individuals. Histograms show the distribution of specimens along the PC1 axis that summarizes species differences (Figure 3). Considering the overall density of the PC1 score within each population sample, we observed variance consistent with convergence of the means where there was evidence of two taxa being present (Figure 3). Population samples from areas of sympatry had mean PC1 close to zero and unimodal distributions (d’Agostino omnibus normality test, overall .21 < p < .98) except in Little Valley Road (AxL) where a bimodal distribution was confirmed (d’Agostino omnibus normality test, p = .007).

FIGURE 3

Grasshopper phenotype suggests a hybrid swarm. Phaulacridium marginale and P. otagoense pronotum shape converge to unimodal morphotype distribution in areas of sympatry. (a) Stacked barplot of assignment probabilities from Gaussian mixture model assignment of pronotum geometric landmark data (PC1 and PC2). (b) Histogram stack plot of morphotype assignments from population samples with more than 20 individuals. Assignment to P. marginale cluster >0.8 (orange), assignment to P. otagoense cluster (turquoise), and putative hybrids with low assignment probability (0.2–0.8; yellow). (c) Density plot of PC1 scores of Phaulacridium grasshoppers at each site with colour gradient by PC1 score. Nine population samples with more than 20 specimens are shown, among which five locations where the two species are present in sympatry (AxM, TeW, AxG, TeE and AxL). Full sampling details are provided in Table S2

DISCUSSION

Hybridization between recently diverged taxa is common, for example due to climate change‐induced range shifts in already modified landscapes, resulting in secondary contact. Although hybridisation has the potential to beneficially increase genetic diversity in small populations or introduce alleles that increase fitness (e.g. adaptive introgression), there is also concern that hybridization may result in biodiversity loss (Sánchez‐Guillén et al., 2013). If assortative mating with respect to species predominates, reproductive isolation can persist despite sympatry (Larson et al., 2013). If, however, mate choice strategies result in weak assortative mating, but fitness of hybrids is low, a stable hybrid zone could emerge because a small reduction in hybrid fitness produces a narrower hybrid zone than a strong but imperfect mating preference for conspecific partners (Irwin, 2020). If premating barriers break down, whether for environmental or evolutionary reasons, numerous hybrids will emerge in the form of a hybrid swarm, potentially replacing both parental species.

Despite an absence of hybridisation less than 50 years ago (Westerman & Ritchie, 1984), our sampling of Phaulacridium grasshoppers in copula shows that heterospecific mating does occur naturally in the wild (Table 1). At locations with both grasshopper species, we found mating was random with respect to morphological cluster (species). Even if we consider that 15% of our specimens might be incorrectly identified our observations were so similar to expectations generated by random mating the same conclusion is drawn; error in species identification could not explain our results. However, our results suggest that these two morphologically distinct species may be less distinct when sympatric. We detect not only a discordance between phenotype and mtDNA haplotype but also the deficit of intermediates in sympatry, (as required by the genotypic cluster definition; Mallet, 1995) . Westerman and Ritchie (1984) inferred that P. marginale and P. otagoense were restricted to parapatric microhabitats but less than 50 years later we found them together. This is consistent with population genetic evidence for recent range expansion of P. marginale associated with ongoing habitat modification and mitochondrial introgression in areas of sympatry (Figure 1; Table 2). Although ecological conditions might influence the rate of interspecific hybridisation, a lack of fitness disadvantage appears to explain our observation of mtDNA introgression/cytonuclear discordance (Rosenthal, 2013; Sivyer et al., 2018).

Some grasshopper specimens collected from sympatry had low assignment probabilities for species identification, which is consistent with hybridisation producing phenotypically intermediate individuals. Morphology can be influenced by environment (e.g. diet) confounding comparisons among populations, and bias from correlated measurement error can also suggest assortative mating is occurring (Class et al., 2017). To accommodate this, we used the shape of the grasshopper pronotum as a proxy for individual phenotype so individuals of different sizes could be compared using a geometric morphometric approach with Procrustes transformation. We note too that measurement error can produce skewed distributions with the potential to yield erroneous detection of bimodality (Freeman & Dale, 2013). We therefore obtained an estimate of measurement error by re‐measuring a subset of individuals multiple times and found it to be negligible (<4% of detected variation). Analysis of samples from single locations allowed us to control for environmental variation within populations where we detected a wide range of size and shape variation suggestive of high genetic diversity. Although the two grasshopper species have distinct pronotum shapes when allopatric, we identified unimodal distributions of pronotum shape in most samples where the two species were sympatric. These unimodal distributions are consistent with a situation where assortative mating is weak (Jiggins & Mallet, 2000) and correspond with our direct observation of adults of the two species in copula in the wild and in captivity. From most sympatric locations, the abundance of intermediate individuals in our samples yielded a continuous range of phenotypes (unimodal) suggesting that F1 hybrids are fertile and are producing backcross offspring. The exception was at Little Valley Road Alexandra (AxL) where there was a signal of two phenotypes (bimodal; Figure 3). The presence of two distinct phenotypes despite catching an interspecific mating pair could reflect local difference in frequency of P. marginale and P. otagoense at Little Valley Road (AxL), and/or more recent contact between the two species here. A single trait such as size or mtDNA lineage could be misleading when investigating assortative mating as many traits have location effects or are easily introgressed (Harrison & Larson, 2014). By combining approaches that examined the process (mate choice) and the outcome of hybridisation (gene flow and phenotype distributions), we revealed a compelling signal that these grasshoppers were randomly mating with respect to species.

The consequence of random mating and backcrossing without selective disadvantage is a localized hybrid swarm characterized by high genetic and phenotypic variability when compared with pure parental species (Allendorf et al., 2001; Langton‐Myers et al., 2019) and novel genetic diversity (Seehausen, 2004). This is consistent with evidence from mtDNA for gene flow between P. marginale and P. otagoense (Table 2) and individuals with intermediate phenotype in sympatry. An increase in the abundance of one species in the range of the other species (or observation at an early stage of the contact) would lead to an excess of parental phenotypes (Jiggins & Mallet, 2000) and a bimodal distribution suggestive of assortative mating. This may explain our observations at Little Valley Road (AxL). However, we found that the distribution of Phaulacridium pronotum shape variation was unimodal in most samples from sympatric locations, meaning that the products of hybridisation (intermediate phenotypes) are overrepresented with respect to parental forms. We assume that at equilibrium the centre of the final Gaussian will depend on the initial local proportion of the two species and selection (Harrison & Bogdanowicz, 1997). The current distribution of morphological shapes is likely influenced by ongoing dispersal of individual grasshoppers into the areas where hybrids are common. Phaulacridium otagoense is surrounded by the widespread P. marginale (Figure 1) so that there are multiple locations of contact between the two. Although most Phaulacridium are flightless, the change detected in less than 50 years suggests rapidly shifting distributions coinciding with agricultural intensification in the region. In the light of intense environmental modification, the lack of reproductive isolation between P. marginale and P. otagoense has led to the formation of hybrid swarms in several regions of sympatry. This could ultimately lead either to local extinction of pure forms of both parental species or to replacement of one lineage by a parental species with introgressed alleles.

CONFLICT OF INTEREST

The authors have no conflict of interest to declare.

AUTHOR CONTRIBUTIONS

MMR and SAT conceived and designed the study and collected grasshoppers. MV photographed all specimens and carried out the geometric morphometric and statistical analyses and extracted DNA. SAT sequenced and analysed mtDNA. MV and SAT created figures. All authors worked on the manuscript, gave final approval for publication and agree to be held accountable for the work presented therein.

PEER REVIEW

The peer review history for this article is available at https://publons.com/publon/10.1111/jeb.13879.

Supplementary Material

Click here for additional data file.