A statistical approach to understanding reproductive isolation in two sympatric species of tree crickets

Reproductive isolation is essential between sympatric species in order to avoid heterospecific matings, which typically result in either no or sterile offspring (Futuyma, 1998; Mayr, 1963). Various pre- and post-copulatory mechanisms exist to ensure reproductive isolation. The pre-copulatory mechanisms ensure that time and effort are not wasted in approaching inappropriate heterospecific mates. For example, in animals that signal their presence and mating status to their mates, the signals themselves may be species specific (Mayr, 1963). Thus it is also essential in signalling species that the receiver recognizes the conspecific signal as its own.

Signallers may employ various strategies to avoid heterospecific signal overlap such as partitioning themselves in time or space. For example, acoustically calling species may call at different seasons, or within a season may call at different times of the day to avoid overlap with heterospecifics (Brumm and Slabbekoorn, 2005; Gottsberger and Gruber, 2004; Greenfield, 1988). Signallers who overlap in their calling period can also partition themselves in space. Animals could distribute themselves in different vertical strata (Diwakar and Balakrishnan, 2007; Hödl, 1977; Sueur, 2002) or in different microhabitats to avoid signal overlap (Hödl, 1977).

Especially for sympatric species that are not partitioned in either time or space, receivers need to recognize the conspecific signal and respond specifically to it (Ryan, 1988). To attract mates, males of orthopteran and anuran species employ species-specific calls (Alexander, 1967). The correct recognition of these calls by potential mates is therefore of primary importance. Among the various temporal features, the micro-temporal feature pulse rate has been indicated to be most important in call recognition (Popov and Shuvalov, 1977; Schildberger, 1984; Thorson et al., 1982; Walker, 1957; Weber et al., 1981). That there was more complexity to the phenomenon was soon shown (Doherty, 1985; Huber et al., 1989). Macro-temporal features such as chirps were also subsequently demonstrated to play a role in call recognition (Gerhardt and Doherty, 1988). Various studies since have shown that both macro- and micro-temporal features play a role in call recognition (Clemens and Hennig, 2013; Hennig and Weber, 1997; Hennig et al., 2014; Pollack and Hoy, 1979). How call pattern recognition might actually take place neurologically has also been investigated (Meckenhäuser et al., 2013; Pollack, 2001). Among the various neurologically driven questions, whether an AND or OR gate exists between the pattern recognizers of macro- and micro-temporal features has been a subject of study (Grobe et al., 2012; Rothbart and Hennig, 2012).

From the characteristics of the features themselves, various hypotheses have been postulated for signal recognition. One such hypothesis, the invariant feature hypothesis (Emlen, 1972; Falls, 1963; Marler, 1960; Nelson and Marler, 1990) predicts that the features of the signal that have low intra-species variance would be the ones used for signal recognition. These features would be least likely to overlap with the signals of the heterospecifics. In contrast to this, the sound environment hypothesis (Bremond, 1976; Brémond, 1978; Dabelsteen and Pedersen, 1985; Emlen, 1972; Marler, 1960; Nelson and Marler, 1990) predicts that the features, which best separate the signals of one species from those of other species, would be the ones most likely to be used for signal recognition. The features may or may not show low variance. A characterization of the signal spaces of the conspecific and the heterospecifics would generate predictions as to which feature seems most likely to be used in species recognition. However, to determine this would require one to study the perception of the receivers.

How the receiver perceives and recognizes the conspecific call can be studied by examining the tuning of the receiver to the various temporal and spectral features. The tuning of the receiver to the signal can be elucidated by studying the preference functions. The preference function describes the variation in the response of the receiver to variation in the signal (Gerhardt et al., 2000; Ritchie, 1996; Wagner et al., 1995). Recognition of conspecifics is thus dictated by preference functions for the various features of the signal. Various studies have enumerated such individual preferences (Ryan and Keddy-Hector, 1992). Many studies have investigated the preference functions for individual features. Across several studies on insects and anurans, preference functions for dominant frequency and syllable period have been found to be stabilizing (Gerhardt, 1991; Gerhardt and Huber, 2002; Shaw and Herlihy, 2000), while those for chirp period and chirp duration have been found to be directional (Gerhardt, 1991; Ryan and Keddy-Hector, 1992). Usually a stabilizing preference indicates that the particular feature is used for species recognition, while a directional preference indicates that the feature plays a role in mate choice (Paterson, 1985).

Other studies have investigated the problem of species recognition using a combinatorial approach, incorporating receiver response based on various features of an acoustic signal. Studying preference functions for a combination of features of a song is necessary to understand the dynamics of tuning of a species to its conspecific call. When the responses to different component features of a signal are analysed combinatorially, the manner in which one feature modifies response to another feature and the amount of influence each feature has on the receiver become apparent. Although a combinatorial statistical approach can be applied in various ways [generalized additive model (GAM): Amezquita et al., 2011; Bentsen et al., 2006; canonical analysis: Brooks et al., 2005; Gerhardt and Brooks, 2009; second-order polynomial regression: Hennig et al., 2016; multidimensional scaling: Ryan et al., 2003], we chose in this study to examine it using a generalized linear mixed model (GLMM) framework. The advantage of a GLMM framework is that, like other linear modelling approaches, it affords us a predictive model. The model essentially gives an equation predicting the response of the animal to the various features incorporated. After obtaining such an equation, various values of the explanatory features can be incorporated to check how the response gets modified. Thus feature values of heterospecifics can be incorporated and response to feature values from different populations of the conspecific can be examined. The equation obtained from a GLMM analysis also informs about the relative weight associated with each feature in determining the overall response.

The pair of species we investigated in this paper are Oecanthus henryi and Oecanthus indicus. These are sympatric tree cricket species found in the same habitat, during the same season of the year and with similar activity periods (Metrani and Balakrishnan, 2005). Oecanthus henryi and O. indicus calling songs overlap significantly in their respective carrier frequencies (Metrani and Balakrishnan, 2005). The two species thus provide an ideal system to study pre-copulatory reproductive isolation through acoustic communication. The response of O. henryi females to male calling song based on the carrier frequency was found to be broadly tuned for the entire natural range of carrier frequencies (Mhatre et al., 2011). Such a broad tuning to a range of carrier frequencies that encompasses the heterospecific values (O. indicus calling song) implies that O. henryi females are unable to distinguish between the signals of O. henryi and O. indicus using spectral features. Temporal features of the signal may therefore be more important for species recognition and reproductive isolation.

The response space for a species can take various shapes when plotted along with the signals of the conspecific and heterospecific (Ryan and Rand, 2001). The response space may be (1) confined to the signal space of the conspecific, (2) encompass both the conspecific and heterospecific signal space, (3) show more preference to the conspecific and reduced but some response to the heterospecific, or (4) show more preference to the heterospecific and reduced response to the conspecific (Ryan and Rand, 2001).

In this paper we have investigated reproductive isolation between O. henryi and O. indicus based on their calling songs and receiver response spaces. We chose three temporal features, chirp period, chirp duration and syllable period (Fig. 1) for studying O. henryi female response to the calling song. Response functions were generated separately for the three temporal features using playback experiments. A response space was then generated for O. henryi females by statistical modelling, integrating the information from the individual response functions. By superimposing the signal spaces of O. henryi and O. indicus on the predicted response space of O. henryi we examined how the predicted response probability of O. henryi females changed between the signal spaces of the conspecific and the heterospecific. We argue that such quantitative models are important to study reproductive isolation based on song as they incorporate information on responses to several different acoustic features of the signal, and also allow for prediction of responses to novel combinations of features.

Oecanthus henryi Chopard 1936 and Oecanthus indicus Saussure 1878 are tree cricket species found in Southern India (Chopard, 1969; Metrani and Balakrishnan, 2005). The animals are found on bushes of the aromatic weed Hyptis suaveolens (Deb et al., 2012; Mhatre et al., 2011). The animals were collected from Hyptis suaveolens bushes in fields near the village of Ullodu, situated near Bangalore, India (13°38′27″N, 77°42′0″E).

Females were caught in the nymphal stage and reared in the laboratory to the final moult. After the final moult, each female was transferred to an individual plastic box (8 cm in diameter and 5 cm in height) with a perforated lid to ensure virginity. The females were used in experiments 15–21 days after the final moult. Virgin females were used in the experiments to ensure high motivation for phonotaxis. Animals were fed ad libitum on apple pieces and maintained on the natural night–day cycle at room temperature. Animal handling for behavioural experiments was done in accordance with national guidelines for the ethical treatment of animals.

Experiments were conducted over a period of one year, from April 2013 to April 2014. Experiments were conducted indoors in an anechoic room (2.9 m×2.75 m×3.1 m) from 18:45 h to 21:15 h (the active period of animals in the field). The phonotaxis experimental set-up consisted of stripped down branches of Hyptis suaveolens arranged in the shape of a T. The T was placed vertically on the ground so that the horizontal part of the T was 60 cm above the ground. Two loudspeakers (X-mini Capsule Speaker V1.1, Xmi Pte Ltd, Singapore) were placed at the horizontal ends of the T (120 cm apart). The vertical branch of the T was placed so that it bisected the horizontal branch at its midpoint, thus the junction was 60 cm from the ends of the horizontal branch of the T [see (Mhatre et al., 2011) for details].

For each trial the animal was released at the bottom of the T. Single speaker playback experiments were conducted, where only one of the two speakers was active at a given time. The mute speaker at the other end was placed to check for directional bias. We performed experiments with a no-choice paradigm as previous data from the wild have indicated that female O. henryi perform sequential rather than simultaneous mate sampling (Deb and Balakrishnan, 2014). In such scenarios, to experimentally elucidate the response function, a no-choice paradigm is preferable (Gerhardt and Huber, 2002). Stimuli were played back at a sound pressure level (SPL) of 61 dB (re 2×10⁻⁵ Nm²) as measured from the junction of the T using a Brüel and Kjær SPL meter (type 2250, Brüel and Kjær ½ in microphone type 4189; Brüel and Kjær Sound and Vibration Measurement A/S, Nærum, Denmark). A trial was scored as positive if the animal reached the playback speaker within 180 s from release. The experiments were conducted in complete darkness and responses were recorded with an IR-sensitive video camera (HDR-CX730E, Sony Corporation, Tokyo, Japan). The anechoic room was maintained at a temperature of 25°C (±1°C), the mean ambient temperature in the field. Throughout the duration of the experiments the temperature in the room was monitored with a temperature meter (Testo 110, Testo AG, Lenzkirch, Germany); a room heater and an air-conditioner (LG, Seoul, Korea) were used to regulate the temperature as and when required. Animals were transferred into the temperature-controlled anechoic room at least an hour before the commencement of the experiments, to acclimatize them to the ambient temperature of the room.

Five sets of experiments were conducted to study response functions. While the temporal features of the calls were varied as per the treatment, the spectral feature (carrier frequency) was kept constant at the value for 25°C, which is 3000 Hz for all the treatments. All calls were synthesized using MATLAB 6.5 (The MathWorks Inc., Natick, MA, USA). A representative syllable of the natural song at 25°C was used to generate the envelope. The envelope was then used to synthesize a syllable at the chosen carrier frequency (3000 Hz). This syllable was repeated at different rates by modifying the intervals between the syllables to generate chirps containing different syllable periods (Fig. 1). Different chirp durations were achieved by modifying the number of syllables in a chirp. Modifying the intervals between successive chirps generated the required chirp periods.

For each of five sets an independent group of animals was tested. Within a set, a repeated measures design was followed, in which each animal was tested in all the treatments of a given set. Each set for an animal was completed on a single night. Between two trials on the same animal, a gap of at least 10 min was maintained. Stimuli were presented to the animals in random order.

The first set of experiments was designed to determine whether the macro- or micro-temporal features were more important for signal recognition in O. henryi. For this set of experiments, six stimuli were created. The first stimulus (Hen) of this set was a representative conspecific O. henryi call (Table S1), which served as the positive control. The second (Ind) was a representative O. indicus call (Table S1). For the next four stimuli the micro-temporal features (syllable period and syllable duration) and macro-temporal features (chirp period and chirp duration) of the calls of O. henryi and O. indicus were interchanged to create novel stimuli (Table S1, Fig. 2). The first three stimuli had the micro-temporal features of O. henryi, and the macro-temporal features of O. indicus. The first (HImed) of these three stimuli had chirp period (CP) and chirp duration (CD) close to the median values found for O. indicus at 25°C (median CP=720.3 ms, median CD=480 ms), the second (HIq3) had values that were close to the third quartile (3rd quartile CP=756.7 ms, 3rd quartile CD=620 ms) and the last (HImax) had values that were close to the maximum (maximum CP=1241 ms, maximum CD=1310 ms) found in O. indicus song. The syllable period and syllable duration had the mean values of an O. henryi call at 25°C. The last stimulus (IH) of this set had the micro-temporal features of an O. indicus call, and the macro-temporal features of an O. henryi call (Table S1).

The next three experiments were conducted to ascertain the female response functions for CP, CD and syllable period (SP), respectively.

The first experiment contained a set of six stimuli where, while every other feature was held constant at the mean value of the O. henryi call (CD: 236.9 ms; SP: 17.3 ms), only the CP was varied. CP was varied both above and below the mean CP of an O. henryi call such that the minimum value was lower than the minimum found in the conspecific calling song at 25°C (minimum CP=500.4 ms), while the maximum value was close to the maximum found in the heterospecific calling song at 25°C (maximum CP=1214 ms). The CPs tested were: 300, 450, 600 (mean), 750, 900 and 1100 ms (Fig. 3).

The second experiment also contained a set of six stimuli, varying only in CD, while the others were held constant at the mean value of the O. henryi call (CP: 633.7 ms; SP: 17.3 ms). Different CDs were obtained by varying the number of syllables per chirp; from two syllables per chirp (CD: 29 ms) to 33 syllables per chirp (CD: 566 ms) (almost a trill), with steps at seven syllables (CD: 115.7 ms), 14 syllables (CD: 236.9 ms) (mean), 21 syllables (CD: 358.2 ms) and 28 syllables (CD: 479.5 ms). While the stimulus with minimum chirp duration was a bisyllabic chirp, the stimulus with maximum duration was close to the median of the chirp duration distribution found in the natural calling song of O. indicus, the heterospecific (480 ms).

A subset of this set included four stimuli to investigate response for chirp duration at finer resolution between seven and 14 syllables per chirp, and 14 and 21 syllables per chirp. The stimuli in this subset were chosen to include more points from within the natural calling song distribution of O. henryi, the conspecific. The stimuli in this subset included 10 (CD: 169 ms), 12 (203 ms), 16 (273 ms) and 18 syllables per chirp (307 ms). The group of animals for this subset was independent from the group used for the second experiment.

The third experiment included five stimuli. In this set only the syllable period was varied, with the chirp period held constant at the mean value of O. henryi call (CP: 633.7 ms). The chirp duration varied slightly (233–242.9 ms) around the mean chirp duration (236.9 ms) as the syllable period was varied and the number of syllables per chirp was held constant at 14. Syllable periods tested were: 14, 17 (mean of O. henryi song), 20, 23 and 25 ms. The minimum value of the syllable period used was 14 ms as the syllable duration at 25°C of O. henryi was 13 ms. The maximum value exceeded the maximum found in the natural calling song of O. indicus (19.22 ms).

In each of the above experiments the stimulus with the mean values of chirp period, chirp duration and syllable period for the respective sets served as the positive control for that set. Every individual tested was found to respond to the positive control, suggesting that females were motivated to respond to male calls.

In the previous sets one temporal feature was varied per set, keeping the others at their mean value. However, if the response depends not only on the additive effect of each individual feature but also on how one feature may modify the effect of another feature, then the previous sets would be unable to capture such an effect. In this set, therefore, two temporal features were co-varied and the response examined. In this experiment, two out of the three temporal features (chirp period, chirp duration and syllable period) were co-varied between the extremes of their already tested values (i.e. minimum and maximum values), while the third one was held constant at the mean value. Thus if chirp period was varied to its maximum at 1100 ms and chirp duration to its minimum at two syllables per chirp, then syllable period would be at its mean value of 17 ms. Eleven such stimuli were created (Fig. 4A). The twelfth stimulus in this set was the positive control, where all three temporal features had the mean values of their respective distributions at 25°C.

All the statistical analyses were carried out using R (version 3.0.3; R Core Team, 2014). In the phonotaxis experiments a female was scored as a positive responder if she went to the playback speaker within 180 s of release. Within a set, each female was tested once on each of the stimuli. The final score for each stimulus was the number of positive responders for that particular stimulus.

For the macro- versus micro-temporal feature experiment, a McNemar chi-squared test was conducted to test whether the relative frequencies of positive and negative responses varied across the different stimuli. For the sets exploring the individual response functions, McNemar chi-squared tests were conducted to compare the relative frequencies of positive and negative responses between the stimulus with the highest response and that with the second highest response. McNemar chi-squared tests were employed because a repeated measures design was used in each set (Sokal and Rohlf, 2012).

For various values of chirp period (300–1100 ms), chirp duration (29–566 ms) and syllable period (14–25 ms), response probabilities were obtained from the model coefficients estimated through the GLMM. The result was a response space of O. henryi based on a combination of responses to chirp period, chirp duration and syllable period. On the generated response space of O. henryi individual calls of O. henryi and O. indicus were superimposed to examine the response probabilities predicted for each of these conspecific and heterospecific signals. The individual predicted response probabilities at each of the O. henryi and O. indicus signal values were also plotted to elucidate the difference in predicted response of O. henryi to the conspecific versus the heterospecific signals.

The signals of Oecanthus henryi that were superimposed on the generated response space were further evaluated. The coefficients of variation of the three features were determined and ranked to examine the relative variability in these features. A linear discriminant analysis was also performed on the signals of both the conspecific and heterospecific to elucidate which of the features (or combination of features) best distinguished the O. henryi signals from the O. indicus signals.

Fifteen animals were tested. All the animals showed positive response to the conspecific O. henryi song (Hen), while only one out of the fifteen responded to the heterospecific O. indicus song (Ind) (Fig. 2). The number of positive responses was low and progressively decreased for HImed, HIq3 and HImax, respectively, in which the stimuli consisted of the macro-temporal pattern of O. indicus, and the micro-temporal pattern of O. henryi. Response to the stimulus (IH) where the macro-temporal pattern was that of O. henryi and the micro-temporal pattern of O. indicus was as low as that for the heterospecific song.

There was a statistically significant difference in the relative frequencies of positive and negative responses between the O. henryi song and the treatment with the next highest number of responders (HImed) (Hen versus HImed, McNemar chi squared=7.1, P=0.007) indicating that macro-temporal features do play a significant role in determining female response. There was no significant difference between the negligible response to the O. indicus song and the song pattern with macro-temporal features of O. henryi and micro-temporal features of O. indicus (Ind versus IH, McNemar chi squared=0, P>0.9999). Thus changing only the micro-temporal features, syllable period and duration, to the values of that of the heterospecific song was sufficient to eliminate phonotactic response. There was no significant difference between the O. indicus song and HImed (Ind versus HImed, McNemar chi squared=3.2, P=0.07).

The variability in signal features of the conspecific Oecanthus henryi was analysed. Syllable period had the least coefficient of variation, followed by chirp duration and chirp period (Table 1). Between the signals of the conspecific (O. henryi) and heterospecific (O. indicus) a linear discriminant function analysis was also performed to assess which, out of the three features, contributed most to the discrimination between the signals. Syllable period had the highest loading, followed by chirp duration and chirp period (Table 1). These three features were sufficient to give rise to separate groups distinguishing O. henryi and O. indicus (Fig. 5).

Eighteen animals were tested. A closed response function was obtained with maximum response for the chirp period of 600 ms, which was closest to the mean chirp period of O. henryi (Fig. 3A), with all animals responding positively. The response to 600 ms differed significantly from the response to 450 and 750 ms, the stimuli which showed the next highest responses (600 versus 450 ms, McNemar chi squared=5.14, P=0.02; 600 versus 750 ms, McNemar chi squared=4.17, P=0.04). Minimum response was obtained for the stimulus with chirp period of 300 and 1100 ms. The response was thus highest for the mean and significantly decreased for values both higher and lower than the mean.

Fifteen animals were tested. Maximum response was obtained for the stimulus with mean O. henryi chirp duration of 14 syllables per chirp (∼237 ms), where all the fifteen animals responded positively (Fig. 3B). Chirp durations lower than the mean hardly elicited any response, while the response to values higher than the mean were significantly lower [14 syllables per chirp (237 ms) versus 21 syllables per chirp (358 ms), McNemar chi squared=4.16, P=0.04]. Responses similar to the mean chirp duration were also found for chirp durations with 10, 12, 16 and 18 syllables per chirp (tested in a different set) (Fig. S1). Response decreased progressively as chirp durations increased from more than 18 syllables per chirp. Thus for chirp duration also, a closed response function was obtained with maximum response at and around the mean chirp duration.

Twenty-one animals were tested. All the animals responded to the mean syllable period of O. henryi, 17 ms (Fig. 3C). Response decreased significantly at the lower value of syllable period, 14 ms (17 versus 14 ms, McNemar chi squared=5.1, P=0.02; 17 versus 20 ms, McNemar chi squared=8.1, P=0.004). As the mean syllable duration used was 13 ms, syllable periods less than 14 ms could not be tested. Female response decreased for syllable periods with values greater than the mean, showing a closed response function.

Twenty animals were tested. A stimulus with the characteristics of the mean call of O. henryi was used as a positive control, and all the animals responded positively to it (twelfth stimulus in Fig. 4B). When the syllable period (SP) was kept constant at its mean value and chirp period (CP) and chirp duration (CD) were at their maximum values (third stimulus), seven females out of 20 responded positively. Holding CD at mean value, 11 females out of 20 responded positively when SP was at its minimum and CP at maximum (sixth stimulus). When CP was kept at its mean value, the number of positive responses elicited was always less than five, as was for the rest of the stimuli (Fig. 4B).

The standard deviation for the random effects (effect of individual animals) in the mixed effects model was 0.88. All the fixed effects terms, including the interaction terms, were statistically significant in explaining the response (Table 2). When chirp period was reduced and increased by 25% of its median value, the changes in the predicted response probability were 0.121 and 0.001, respectively. The corresponding values for chirp duration were 0.056 and 0.010, while for syllable period they were 0.034 and 0.289. The effect size for syllable period was the most, indicating that syllable period had more influence than the other features in describing the response space.

The predicted response space of O. henryi for a combination of chirp period, chirp duration and syllable period is shown in Fig. 6. The response probabilities vary from 0.1 to 0.9 (shown by varying colours). Individual calls of O. henryi and O. indicus are also plotted on the response space, indicating the predicted response for the calls. Most of the O. henryi calls lie within the 0.8–0.7 predicted response probability, whereas the O. indicus calls lie mostly at the periphery (0.5–0.1 predicted response probability) (Fig. 6).

From this response space, individual response probabilities for each of the calling songs of both the species were also estimated. The calling songs were plotted as a function of their similarity (measured in terms of distance) from both the mean song of O. henryi and O. indicus. On these was superimposed the predicted probability value for each of these calls (varying colours) (Fig. 7). The predicted response probabilities to the various O. henryi songs varied from 0.69 to 0.86, with a mean of 0.8; the response probabilities to the various O. indicus songs varied from 0.001 to 0.6, with a mean of 0.32. The estimated response probability to the mean song of O. henryi (CP=629 ms, CD=241 ms, SP=17 ms) was 0.83 and that to the mean song of O. indicus (CP=720 ms, CD=480 ms, SP=21 ms) was 0.3. In general the estimated response probability of a song was higher if it was closer to the mean song of O. henryi and lower if songs were closer to the mean song of O. indicus. Songs of O. indicus that were further from both the mean songs had the lowest response probabilities (Fig. 7).

Stabilizing preference on features of a signal have traditionally been associated with aiding in mate recognition (Paterson, 1985). These features are usually also associated with low variability (Gerhardt, 1991). Of the temporal features, syllable period has often been found to be under stabilizing selection (Popov and Shuvalov, 1977; Stout and McGhee, 1988; Gerhardt, 1991; Wollerman, 1998; Shaw and Herlihy, 2000; Gerhardt and Huber, 2002; but see Hedrick and Weber, 1998). Chirp period and duration, however, have been found to be under directional selection (Hedrick, 1986; Gerhardt, 1991; Ryan and Keddy-Hector, 1992; Wagner et al., 1995; Gerhardt et al., 2000; Murphy and Gerhardt, 2000; Nandi and Balakrishnan, 2013; but see Hedrick and Weber, 1998). These features are also associated with having higher variability (Gerhardt, 1991). In studies of multivariate preference functions of females, Blankers et al. (2015) and Hennig et al. (2016) found stabilizing preference for the micro-temporal feature pulse rate and directional preference for the macro-temporal feature trill/chirp duration. We found closed response functions for both micro- and macro-temporal features.

However, recent studies of multivariate selection using the fitness landscape concept of genotypic space (Wright, 1931) have shown that features previously thought to be under directional selection can also contribute to stabilizing selection. Temporal features of the calls of the cricket Teleogryllus commodus such as trill number (number of trills), chirp pulse number and inter-call duration were found to be under multivariate stabilizing selection both in laboratory (Brooks et al., 2005) and field experiments (Bentsen et al., 2006). The coefficient of variation of these temporal features was around 12% for chirp pulse number (Bentsen et al., 2006; Brooks et al., 2005), around 40% for trill number (Bentsen et al., 2006; Brooks et al., 2005) and for inter-call duration varied between 12% (Brooks et al., 2005) and 86% (Bentsen et al., 2006). A similar study on the frog Hyla versicolor found directional selection on call duration and call rate, but found stabilizing selection on the combination of pulse rate, call duration and call rate (Gerhardt and Brooks, 2009). The coefficient of variation on pulse rate was lowest at 8%, while call duration was 24% and call rate was 17%. The coefficient of variation we obtained for chirp duration and chirp period was comparable to the coefficient of variation in similar features of previous studies (Table 1). We also obtained closed response functions for chirp duration and chirp period, and find all the three temporal features contributing significantly to the multivariate response space for O. henryi. Call features with relatively high variability may also contribute to stabilizing selection when analysed in a multivariate framework. The dichotomy of invariant features for mate recognition and variable features for sexual selection may not be necessarily true, though more studies of such kind are needed to explore this further. In our study, the closed response functions on chirp period and duration may indicate selection pressure on these features not to approach values typical of the same features of the sympatric heterospecific songs, thus achieving pre-copulatory reproductive isolation.

The first phonotactic experiment showed that both the macro- and micro-temporal patterns are important for conspecific song recognition. However, the difference in O. henryi female response to conspecific and heterospecific song was more pronounced when the signals differed in their micro-temporal pattern, thus we infer that the micro-temporal pattern plays a more important role in signal recognition than the macro-temporal patterns. The importance of micro-temporal features was also confirmed in the mixed model analysis, as the syllable period had the highest influence of all the temporal features. Syllable period or pulse rate has traditionally been regarded as the most important feature in song recognition (Popov and Shuvalov, 1977; Walker, 1957; Weber et al., 1981), to the extent that it was once regarded as the sole temporal feature necessary for song recognition (Schildberger, 1984; Thorson et al., 1982). However, studies have subsequently shown other features to be necessary as well (Doherty, 1985; Huber et al., 1989). In this study we found that acoustic species recognition in O. henryi is dictated by all the three temporal features individually and also by their combined effect, with syllable period being the most important temporal feature.

Our findings characterizing the features according to intra-species variance and how much they contributed in distinguishing the conspecific from the heterospecific did not help us to distinguish between the invariant feature hypothesis and sound environment hypothesis. Our results actually provided evidence for both the hypotheses. Syllable period was both the least variable and the one that had the highest loading in the linear discriminant analysis (Table 1). Interestingly, syllable period also had the highest influence in describing the response space.

We used statistical modelling to generate a response space for O. henryi. Response functions based on phonotactic experiments for the three temporal features as well as their interactions were integrated using a generalized linear mixed models approach. While relatively rare, other studies have similarly used a statistical approach to describe response spaces. Amezquita et al. (2011) used a generalized additive model framework with one spectral and one temporal feature, and the response space they described for each anuran species was restricted to a probability threshold of 0.9. Ryan et al. (2003) considered four temporal and three spectral features, and reduced it to two dimensions by multi-dimensional scaling. They then used quadratic equations to arrive at a generalization gradient for response. To show how song recognition takes place in three species of field crickets, Blankers et al. (2015) represented the results of individual preference functions of pulse rate and trill duration in a bivariate response space. They used linear modelling for constructing the bivariate response space. In studies of multivariate selection on the songs of the field cricket Teleogryllus commodus (Bentsen et al., 2006; Brooks et al., 2005) and of the anuran species Hyla versicolor (Gerhardt and Brooks, 2009) artificial songs were created where various song features were varied randomly from a normal distribution of the same features. The response to these songs was then used to describe matrices of directional, quadratic and correlational selection using a multiple regression framework. Further, canonical analysis was applied and the matrices were rotated so that along each major axis (number of major axes determined by the number of features varied) the value for directional selection as well as eigenvectors for non-linear selection was obtained. In a recent paper Hennig et al. (2016) explored the multivariate response of three allopatric field cricket species. Similar to our approach, from experiments elucidating the individual preference functions for carrier frequency, micro- and macro-temporal features, a multivariate response space for all the three species was created. Second-order polynomial regression was used. From the model, the differential weighting of each feature and interactions, with similarities and differences between the three species, was expounded (Hennig et al., 2016).

We attempted to expand the statistical approach taken by previous studies in several ways. We did not restrict the response probability to any particular threshold. We used GLMM, which allows for greater flexibility while working with non-normal binomially distributed data. The GLMM approach also allowed the construction of a predictive model to estimate the response of O. henryi females to different combinations of values of chirp period, syllable period and chirp duration, as well as to heterospecific calls. It also allowed us to factor in interactions between the features.

In general, the model predicted high response towards conspecific calls and much lower response to calls in the heterospecific response space. The stochastic nature of our model, however, gave rise to a few constraints. The response probability for the mean conspecific calling song as predicted by our model was found to be less than the value observed in our experiments (0.83 versus 1; Fig. 1). Similarly, the response probability for the mean heterospecific song predicted by our model was higher than that obtained in phonotaxis experiments (0.3 versus 0; Fig. 1). These caveats could be due to the nature of the quadratic equation utilized in the mixed model. As a quadratic equation is symmetrical, it approximates the humped shape accurately but is limited in accommodating the extreme values. Future studies could further evaluate the details of the shape of the response functions to better capture female response to call features.

The statistical model of response space constructed for O. henryi indicates maintenance of reproductive isolation between O. henryi and O. indicus through the temporal features of their calling songs. The predicted probability of response of the conspecific (O. henryi) calls was consistently higher than those of the heterospecific (O. indicus) calls. Thus the predicted response is not entirely limited to the conspecific signal space but also extends to the heterospecific signal space but with lower probabilities. This corresponds to the third scenario described in Ryan and Rand (2001) where the response space is neither completely limited to the conspecific, nor does it extend equally over both conspecific and heterospecific signals. Various other studies have also shown response of females to heterospecific songs (Backwell and Jennions, 1993; Doherty and Howard, 1996; Oldham and Gerhardt, 1975).

However, as indicated above, the predicted response probability to the mean song of O. indicus is high when compared with the behaviour. The model is, however, useful because it indicates which combinations of the temporal features found in the heterospecific songs elicit relatively high response probabilities and hence need to be tested. Testing females on combinations of features, which were not included for the model generation, would also serve to externally validate this model and this is planned to be done in future experiments. This model, moreover, can be used more generally and is not restricted to the reported pair of sympatric species. The response of O. henryi to various combinations of chirp period, chirp duration and syllable period can be evaluated, in effect testing the response to other heterospecific signals. This approach could also be used to test and construct similar multivariate response spaces of other cricket species. It would be very interesting indeed to construct the same for the females of the sympatric heterospecific O. indicus, and compare the response spaces between the two sympatrics. In O. indicus where the signal variability is more than that in O. henryi, one could expect the space for the high response probabilities to be broader.

The predicted response space, by delimiting the combinations giving rise to high response probabilities, has implications for the conspecific signal. It places constraints on the variation of the signal allowed within the species. The response to signal spaces of conspecifics from different populations can be evaluated using this model. It can then indicate how the signal has changed, if it has, between different populations and the reasons behind the change can be further explored.

We thank Manjunatha Reddy for helping with collection and maintenance of animals. M.B. is also thankful to Sandeep Pulla for help with R, Dipanjan Bhattacharya for help with MATLAB, and Diptarup Nandi and Rittik Deb for insightful discussions.