Divergence in female damselfly sensory structures is consistent with a species recognition function but shows no evidence of reproductive character displacement

Abstract Males and females transmit and receive signals prior to mating that convey information such as sex, species identity, or individual condition. In some animals, tactile signals relayed during physical contact between males and females before and during mating appear to be important for mate choice or reproductive isolation. This is common among odonates, when a male grasps a female's thorax with his terminal appendages prior to copulation, and the female subsequently controls whether copulation occurs by bending her abdomen to complete intromission. It has been hypothesized that mechanosensory sensilla on the female thoracic plates mediate mating decisions, but is has been difficult to test this idea. Here, we use North American damselflies in the genus Enallagma (Odonata: Coenagrionidae) to test the hypothesis that variation in female sensilla traits is important for species recognition. Enallagma anna and E. carunculatum hybridize in nature, but experience strong reproductive isolation as a consequence of divergence in male terminal appendage morphology. We quantified several mechanosensory sensilla phenotypes on the female thorax among multiple populations of both species and compared divergence in these traits in sympatry versus allopatry. Although these species differed in features of sensilla distribution within the thoracic plates, we found no strong evidence of reproductive character displacement among the sensilla traits we measured in regions of sympatry. Our results suggest that species‐specific placement of female mechanoreceptors may be sufficient for species recognition, although other female sensory phenotypes might have diverged in sympatry to reduce interspecific hybridization.

Although much is known about the importance of visual, auditory, and chemical signals and responses in sexual communication and species recognition, we know relatively little about other sensory modalities that may have strong effects on individual mating decisions. Tactile signals have been hypothesized to contribute to mating decisions (Mendelson & Shaw, 2012), but it is unclear whether tactile cues could represent a primary species recognition signal, given that visual, auditory, and chemical cues usually act earlier during the mating sequence. Research on the prevalence of tactile signals in mating decisions is limited (Coleman, 2008) partly because of the experimental challenge it poses. Whereas other sensory modalities present male signals to a focal female from a distance, studying female preference for tactile cues requires contact between males and females, which is not always easily achieved or quantified under controlled conditions. Despite this challenge, understanding the role of tactile signals along the continuum between intraspecific mate choice and interspecific RI is important because it broadens our understanding of the causes and consequences of a common pattern in nature-the rapid divergence of male genital morphology between species (Eberhard, 1985). It has been suggested that rapid genital differentiation can cause RI (Dufour, 1844), although mechanical incompatibilities between heterospecific male and female genitalia do not appear to be a common cause of RI (Masly, 2012;Shapiro & Porter, 2013;Simmons, 2018). However, observations both within Eberhard, 1994;Edvardsson & Göran, 2000;Frazee & Masly, 2015) and between species (Barnard et al., 2017;Coyne, 1993;Eberhard, 1992;Patterson & Thaeler, 1982;Robertson & Paterson, 1982) suggest that male reproductive structures may convey tactile information to females that affects their subsequent behavior and/or physiology. Although female genital structures often appear invariant among closely related species (Shapiro & Porter, 2013), subtle morphological differences (e.g., Kamimura & Mitsumoto, 2011;Yassin & Orgogozo, 2013) could enable females to detect variation among males' morphology. Female variation in detection ability could also occur in signal processing at the level of neurons, neural networks, and/or in the distribution and morphology of sensory structures that receive male tactile signals.
These sensory structures may exist not just in the female genitalia or reproductive tract, but in any region of the female that receives contact from male structures.
Concentrations of cuticular mechanoreceptors (sensilla) on the female thorax have been described in several coenagrionid damselfly genera Jurzitza, 1974Jurzitza, , 1975Robertson & Paterson, 1982;Tennessen, 1975). The morphology of these sensilla is consistent with a mechanosensory function and does not indicate that they conduct signals related to olfaction, hygroreception, or temperature reception (McIver, 1975). These sensilla reside in areas where males' grasping appendages contact the female thorax before and during mating, which has led to speculation that they allow females to evaluate male morphologies and discriminate conspecific from heterospecific males. This idea is based on demonstrated reductions in female receptivity when grasped by males with manipulated appendages (Loibl, 1958;Robertson & Paterson, 1982) or heterospecific or hybrid males (Barnard et al., 2017;Sánchez-Guillén et al., 2016;Sánchez-Guillén et al., 2016;Tennessen, 1975). In African Enallagma, species-specific placement of sensilla within the female mesostigmal plates appears to correspond to where they are grasped by the male (Robertson & Paterson, 1982), which further suggests that the sensilla receive tactile cues based on male morphology that aid in species recognition.
In insects, each cuticular mechanoreceptor is associated with a single sensory neuron (McIver, 1975;Keil, 1997). The thoracic sensilla thus represent a spatial matrix that can transmit signals to the female central nervous system based on the pattern in which the sensilla are stimulated. Greater numbers of these receptors are expected to enhance a female's sensory resolution by increasing the combinatorial complexity of tactile signals that she can perceive. For example, if a female possesses 25 sensilla, and each sensillum has two response states ("on" if contacted and "off" if not contacted), then the number of unique tactile patterns that the female could distinguish is 2 25 = 3.4 × 10 7 . A female that possesses just one additional sensillum would be able to distinguish among roughly twice as many tactile patterns (2 26 = 6.7 × 10 7 ). Should individual sensilla respond to quantitative variation in touch (rather than a binary response), this would dramatically increase the number of response states and therefore further enhance tactile acuity (e.g., Gaffin & Brayfield, 2017). Female damselfly thoracic sensilla thus present an external, quantifiable phenotype in which to investigate the mechanistic basis of tactile stimuli and female mating decisions.
The North American damselfly genus Enallagma includes several recently diverged species that often co-occur in the same habitats (Johnson & Crowley, 1980;McPeek, 1998), and do not engage in premating courtship (Barnard et al., 2017;Fincke, Fargevieille, & Schultz, 2007) or use chemical cues for mate selection (Rebora et al., 1998). A female's first opportunity to assess a potential mate occurs when the male uses his terminal appendages to grasp the mesostigmal plates on the female thorax to form "tandem," the premating position. The male superior grasping appendages (cerci) have speciesspecific morphologies, and differences in the morphology of these structures are the primary cause of RI in this genus (Paulson, 1974;Barnard et al., 2017). Two species, Enallagma anna and Enallagma carunculatum, have strikingly different cercus morphologies, yet occasionally hybridize in nature to produce males and females with reproductive structure morphologies that are intermediate to each of the pure species (Barnard et al., 2017;Donnelly, 2008;Johnson, 2009;Miller & Ivie, 1995). Females of both pure species discriminate strongly against both heterospecific and interspecific hybrid males that take them in tandem, which indicates that female E. anna and E. carunculatum can detect not only large differences in species-specific male stimulation, but also more subtle differences such as those that distinguish conspecific and hybrid males (Barnard et al., 2017).
Because it appears that mesostigmal sensilla mediate species recognition, they might be expected to show signs of reproductive character displacement (RCD): increased divergence of traits involved in RI in regions of sympatry between E. anna and E. carunculatum relative to regions of allopatry (Brown & Wilson, 1956;Howard, 1993;Pfennig & Pfennig, 2009). RCD can manifest phenotypically as divergence in either signaling traits or mate preferences in which sympatric females display stronger discrimination against heterospecific males than do allopatric females of the same species (e.g., Gerhardt, 1994;Gabor & Ryam, 2001;Albert & Schluter, 2004;Wheatcroft & Qvarnström, 2017). This strengthening of preference in sympatry may evolve via direct selection on adult prezygotic phenotypes, or via reinforcement, where selection against interspecific hybrids gives rise to selection for enhanced premating isolation between species (Dobzhansky, 1937). Enallagma anna and E. carunculatum can interbreed, but their hybrids experience significantly reduced fitness (Barnard et al., 2017). Female Enallagma experience frequent mating attempts from heterospecific males where both species co-occur (Paulson, 1974;Fincke et al., 2007;Barnard et al., 2017). These findings suggest that in sympatry, females may experience selection for stronger species discrimination ability. Studies of several Enallagma species (not including E. anna or E. carunculatum) have revealed that male cercus shape varies little among populations, even across large geographical regions (McPeek, Symes, Zong, & McPeek, 2011;Siepielski, McPeek, & McPeek, 2018). Enallagma anna and E. carunculatum appear to show similar patterns, at least in the western part of their distributions ( Figure S1, Supporting information). It is possible, however, that females in sympatry with other species are more sensitive to variation among males than are females of the same species in regions of allopatry, and this variation in sensitivity may be reflected in female sensilla traits.
Here, we use sensilla number, density, and location as proxies for female preference, to test the hypothesis that variation in female sensilla phenotypes supports a function in species recognition.
We tested this hypothesis by quantifying sensilla on the mesostigmal plates of a large set of E. anna and E. carunculatum females from multiple populations across the western United States and comparing phenotypes of each pure species from sympatric and allopatric populations to identify patterns consistent with RCD. We predicted that in sympatric populations, females would possess higher sensilla numbers, higher sensilla density, and/or different spatial distributions of sensilla within their mesostigmal plates when compared to females from allopatric populations.

| Population sampling
We measured the sensilla traits of 29 E. anna females across 13 populations, and 74 E. carunculatum females across 28 populations ( Figure 1, Table 1). We classified each population as allopatric, locally allopatric, or sympatric. Sympatric populations are those where E. anna and E. carunculatum co-occur temporally as well as spatially.
Because E. anna's geographic range falls completely within E. carunculatum's range (Figure 1), only E. carunculatum has completely allopatric populations. We designated populations as locally allopatric at sites within the area of range overlap, but where only one species is known to occur based on occurrence data from OdonataCentral. org. Specimens were either dried or preserved in ethanol; neither preservation method alters the morphology of the hard cuticle that comprises the mesostigmal plates, nor the ability to visualize sensilla. Although some specimens were collected as early as 1945, the majority of samples (82 of 103) we studied were collected between 2012 and 2016.

| Trait imaging and quantification
We photographed each damselfly using a Nikon D5100 camera (16.2 MP; Nikon Corporation, Tokyo, Japan). We dissected the ventral thoracic cuticle from each female using forceps and imaged the mesostigmal plates using scanning electron microscopy ( Figure 2). Specimens were mounted on aluminum stubs with carbon tape, sputter-coated with gold-palladium, and imaged at, 200× magnification and 3 kV using a Zeiss NEON scanning electron microscope.
To avoid any potential bias during measurements, we blindcoded image files before measuring all sensilla traits. We measured abdomen length (abdominal segments 1-10, excluding terminal appendages) on the full-body photographs as an estimate for body size using the segmented line tool in ImageJ (Abramoff, Magalhaes, & Ram, 2004). We quantified sensilla traits on the right mesostigmal plate of each female damselfly unless the right plate was dirty or damaged, in which case we quantified the left plate (n = 20).
Sensilla counts on a subset of 57 females showed that left plate and right plate sensilla counts were highly correlated within individual females (r = 0.85). In cases where we quantified the left plate, we flipped the image horizontally, so it was in the same orientation as a right plate. We standardized the position of the mesostigmal plate in each image by cropping and rotating the image so that the lower medial corner of the plate was in line with the lower left corner of each image. We counted sensilla and obtained their x and y coordinates in ImageJ using the multipoint selection tool. When a sensillum breaks, it leaves behind distinctive round area with a central pore. In these cases, we recorded the area where the sensillum had been (Figure 2b-c).
We traced an outline around the plate image, excluding the lateral carina (Figure 2f), using a Wacom Cintiq 12WX tablet and stylus (Wacom, Saitama, Japan) and the freehand selection tool in ImageJ.
This procedure produced x and y coordinates that described the plate outline. We performed all measurements twice for each specimen.
Measurements across the two technical replicates were highly correlated (r abdomen = 0.95, n = 78; r count = 0.95, n = 103; r plate area = 0.99, n = 86), so we used the mean trait values of the two replicates in subsequent analyses. Seventeen samples were imaged at angles that allowed counting of the sensilla, but distorted the plate shape or distances between the sensilla. Those samples are included in analyses of sensilla number but were not included in the analyses of sensilla density or distribution.

| Sensilla trait analyses
We conducted all morphometric and statistical analyses using R v.
3.4.1 (R Core Team, 2015). We used the mesostigmal plate outline coordinates to calculate each plate's two-dimensional area. To calculate the area of the sensilla-covered region of each plate, we generated a polygon connecting the coordinates of the outermost sensilla and calculated the area within this outline. We determined the proportion of each plate that was covered by sensilla by dividing the sensilla area by total plate area. We calculated sensilla density in two ways. First, we divided sensilla number by the area of the sensilla-covered region. This measures the number of sensilla that occur in a particular area, but does not capture the relative arrangement of sensilla within that area.
Second, we computed the nearest neighbor distances among all sensilla within each plate based on their x and y coordinates and then calculated the mean and median nearest neighbor distances between the sensilla for each female. Nearest neighbor mean and median distances were highly correlated (r E. carunculatum = 0.83; r E. anna = 0.81), so we used the mean values for these measures in our analyses.
To determine whether larger females possess more sensilla, we regressed sensilla number against abdomen length. We found no significant relationship between these traits in either species (E. anna: We thus present the results that compare sensilla counts without correcting for differences in body size.

| Sensilla spatial analyses
To quantify sensilla distributions within each plate, we generated kernel density estimates (KDEs) for populations with at least four sampled individuals (six E. carunculatum populations, Table 3; and two E. anna populations, both sympatric) using the R package ks (Duong, 2016). First, we randomly selected one of the two replicate sets of sensilla and plate outline coordinates for each female.
To prepare the coordinate data for KDE analyses, we concatenated the sensilla and plate coordinates for each female and adjusted all F I G U R E 1 Sampling sites and species ranges. Enallagma anna's geographic range (red) occurs within Enallagma carunculatum's geographic range (orange). Names of sites associated with each number are described in Table 1. Symbol color indicates the species sampled and symbol shape indicates the population type. (Species ranges are adapted from Johnson, 2009;Paulson, 2009Paulson, , 2011 plate outlines to have an area of one. This standardized each set of sensilla coordinates for size, while maintaining their relative positions within each plate. Next, we translated each set of coordinates to place the origin of the coordinate system at the plate outline's centroid. We concatenated sensilla coordinates for all females sampled within each population to compute a representative KDE for each population. We compared sensilla distributions among E. carunculatum populations using pairwise KDE tests using the function kde.test with the default settings. This test returns a p-value that reflects the probability of generating the two respective KDEs from the same distribution of points. Because we performed multiple pairwise tests among E. carunculatum populations, we adjusted the resulting p-values using the false discovery rate (Benjamini & Hochberg, 1995 and the remaining 198 points were designated as sliding semilandmarks). We used the R package geomorph (Adams & Otarola-Castillo, 2013) to obtain an average two-dimensional plate shape for each population via general Procrustes analysis (Rohlf, 2011), which calculates a mean shape from the landmarks of a set of superimposed shapes.

| Statistical analyses
Some populations were well sampled, whereas others were represented by a single female (Table 1)
The sensilla also occurred in different locations on the mesostigmal plates of each species: They were more medially located in E. anna and more laterally located in E. carunculatum (Figures 3 and 4).

| E. carunculatum sensilla traits do not show a strong pattern of reproductive character displacement
We made several nonmutually exclusive predictions expected under RCD for the sensilla traits we measured in sympatric populations relative to allopatric populations. In particular, we predicted to observe at least one of the following phenotypic differences in sympatric females relative to allopatric females: (a) more numerous sensilla, (b) denser sensilla, and (c) sensilla concentrated in different regions of the mesostigmal plates. We did not find significant differences in any of these traits between sympatric and locally allopatric E. anna females (Table 2). However, because our E. anna samples included only four females from three locally allopatric populations, we could not perform a robust comparison of E. anna sensilla traits between populations that do or do not encounter E. carunculatum.
We thus focus our analysis on comparisons between sympatric and allopatric E. carunculatum populations, for which we had larger sample sizes.
Sympatric, locally allopatric, and fully allopatric E. carunculatum populations did not differ significantly from one another in sensilla number (Kruskal-Wallis 2 2 = 0.69, p = 0.71), proportion of the mesostigmal plate covered by sensilla (Kruskal-Wallis 2 2 = 2.16, p = 0.34), or sensilla density (overall density: Kruskal-Wallis 2 2 = 0.12, p = 0.94; mean distance between sensilla: Kruskal-Wallis 2 2 = 3.53, p = 0.17). In addition to divergence of mean trait values, RCD can also result in reduced trait variance in sympatry without affecting the mean (Pfennig & Pfennig, 2009). Sympatric E. carunculatum populations displayed less interpopulation variance than allopatric populations in both mean sensilla number (Figure 3a) and mean proportion of the plate covered by sensilla (Figure 3b). F I G U R E 4 Individual trait values for sensilla number, sensilla density, and proportion of plate containing sensilla. Each symbol represents a single female, separated by population along the y-axis. Horizontal lines indicate the mean value for each population type (completely allopatric, locally allopatric, or sympatric), calculated from population means. Populations are described in

TA B L E 2 Statistical comparison of sensilla traits in locally allopatric and sympatric Enallagma anna populations
However, these trends were not statistically significant (sensilla number: Bartlett's K 2 1 = 0.83, p = 0.36; proportion of plate covered by sensilla: Bartlett's K 2 1 = 1.86 p = 0.17). Interestingly, although mean trait values did not differ significantly between sympatric and allopatric populations, sensilla traits displayed considerable variation within the populations we sampled. For example, within a single population, a particular female might have twice as many sensilla than another female ( Figure 4). This pattern was also observed in the E. anna populations we studied.
Kernel density estimates comparisons did not reveal significant differences in sensilla distributions between sympatric and allopatric E. carunculatum populations (Table 3)

| D ISCUSS I ON
Enallagma anna and E. carunculatum females possess different numbers of sensilla in species-specific distributions on their mesostigmal plates. This result supports the idea that receptors that receive male stimuli will occur in patterns that correspond to the male organs during contact (Eberhard, 2010). An association between male morphology and female sensilla has been described for African Enallagma species (Robertson & Paterson, 1982), and our results show a similar pattern for two North American species. Enallagma anna male cerci are considerably larger than E. carunculatum cerci, and the observation that E. anna females had a larger number of sensilla compared to E. carunculatum females is consistent with the likelihood that E. anna male cerci make greater spatial contact with the mesostigmal plates.
We might, however, expect a stronger pattern of RCD in sympatric E. anna females because E. carunculatum males can take them in tandem relatively easily, whereas E. anna males are usually unsuccessful at taking E. carunculatum females in tandem (Barnard et al., 2017). This means that E. anna females may have more opportunities for mating mistakes than E. carunculatum females, which can result in stronger asymmetric RCD (Lemmon, 2009;Pfennig & Pfennig, 2009).
There are at least three potential explanations for the absence of RCD in the form of significant differences in the sensilla traits we  (Siepielski et al., 2018). Although RCD is most easily facilitated when the trait under selection already differs between species (Pfennig & Pfennig, 2009), these sensilla traits may have already diverged sufficiently enough to preclude strong selection on further divergence.
Second, it is possible that the external sensilla phenotypes we measured are not representative of proximate female sensory traits, and the variation that directs mating decisions occurs elsewhere within the female nervous system. For example, individual sensilla might differ in firing rate or sensitivity to pressure applied by the cerci. Any of these variables could differ between species or within the same species in allopatry and sympatry without noticeable differences in sensilla morphology.The direction of mechanosensor deflection is also important for stimulus detection (Keil, 1997), and different species' cercus morphologies may contact sensilla from different angles. Female mate preferences may also be influenced by the relative frequencies with which females encounter heterospecific and conspecific males and female sexual experience (e.g., Svensson, Runemark, Verzijden, & Wellenreuther, 1989).
Finally, although we did not detect a statistically significant difference between group means, the small differences we observed may still have biological relevance. If gaining just one additional mechanosensor can (at least) double a female's tactile discriminatory power (Gaffin & Brayfield, 2017), then females in a population with a seemingly minor upward shift in sensilla number could gain a substantial increase in their ability to detect and avoid mating with heterospecifics. Similarly, it is difficult to determine the features of sensilla density distributions that may influence female preference solely by conducting statistical tests between KDEs. Small spatial differences within largely similar patterns may not contribute a signal large enough to be captured in a statistical test, but still reflect salient variation in the way females receive tactile stimuli. This might include three-dimensional spatial differences that we were unable to measure here. Additionally, small sample sizes from some populations and high variation within populations may have limited our ability to detect evidence of RCD. We pooled geographically widespread populations because fragmenting the samples for regional comparisons would have compromised our statistical power to detect differences between sympatry and allopatry.
These possible explanations highlight the interesting avenues that female damselfly sensilla provide for investigating the mechanisms underlying how females evaluate male tactile signals to make mating decisions. The ability to quantify the number and locations of female mechanoreceptors in a region contacted by male reproductive structures complements our understanding of patterns of variation in male morphologies (Barnard et al., 2017;McPeek et al., 2011;McPeek, Shen, & Farid, 2009;McPeek, Shen, Torrey, & Farid, 2008). Females of both species display substantial intrapopulation variation in sensilla traits (Figure 4), and this variation may play a role in sexual selection and female preferences within species. Behavioral studies will be crucial to link mechanoreceptor phenotypes to female mating decisions and clarify whether and F I G U R E 5 Population kernel density estimates for E. carunculatum (a) and E. anna (b) sensilla. The shading indicates different regions of sensilla density: red represents the 75-99th percentile of sensilla density, orange represents the 50-74th percentile, and yellow represents the 25th-49th percentile. Each outline represents the average mesostigmal plate shape for the population. Asterisks indicate E. carunculatum populations whose KDEs are significantly different (*p < 0.05, ***p < 0.001). how sensilla traits influence both species recognition and sexual selection. For example, do females with more sensilla make fewer mating mistakes than females with fewer sensilla (Lemmon, 2009)?
Another outstanding question of this system is how the cerci stimulate individual sensilla during tandem. This might be determined by flash-freezing male-female tandem pairs and using micro-CT scanning to understand how the male and female structures interact, similar to a recent approach used in seed beetles (Dougherty & Simmons, 2017). Once we understand how cerci contact the sensilla, functional tests of sensilla electrophysiology could reveal how individual sensilla respond to stimulation and indicate whether certain sensilla make greater contributions to reproductive decisionmaking than others.
Female preference can drive sexual selection, promote trait divergence, and cause RI between species (Ritchie, 1996). A longstanding presumption in the literature on genital evolution and speciation has been that female reproductive morphologies are less variant or species-specific than male genitalia (reviewed in Shapiro & Porter, 2013). However, recent studies of variation in female reproductive structures suggest that variation does exist among individuals and species (Ah-King, Barron, & Herberstein, 2014), and our data highlight the importance of looking beyond the easily quantified external morphologies. When male reproductive structure morphologies are obviously divergent, but female morphologies are not, females may possess important variation at neurophysiological levels that affect how they evaluate male tactile signals, similar to the way females evaluate signals in other sensory modalities.

ACK N OWLED G M ENTS
We are grateful to P. Larson for SEM imaging help. We thank O. This work was supported by funds from the University of Oklahoma.
We are grateful to assistance with publication costs from OU's University Libraries Open Access fund.

CO N FLI C T O F I NTE R E S T
The authors declare that there is no conflict of interest regarding the publication of this article.

AUTH O R ' CO NTR I B UTI O N S
AAB and JPM planned the study; AAB collected the data and performed the analyses; AAB and JPM wrote the paper.

DATA ACCE SS I B I LIT Y
Data associated with this manuscript have been archived in Dryad Digital Repository (provisional https://doi.org/10.5061/dryad.vf807sd).