Physical linkage and mate preference generate linkage disequilibrium for behavioral isolation in two parapatric crickets.

Behavioral isolation is a potent barrier to gene flow and a source of striking diversity in the animal kingdom. However, it remains unclear if the linkage disequilibrium (LD) between sex-specific traits required for behavioral isolation results mostly from physical linkage between signal and preference loci or from directional mate preferences. Here, we test this in the field crickets Gryllus rubens and G. texensis. These closely related species diverged with gene flow and have strongly differentiated songs and preference functions for the mate calling song rhythm. We map quantitative trait loci for signal and preference traits (pQTL) as well as for gene expression associated with these traits (eQTL). We find strong, positive genetic covariance between song traits and between song and preference. Our results show that this is in part explained by incomplete physical linkage: although both linked pQTL and eQTL couple male and female traits, major effect loci for different traits were never on the same chromosome. We suggest that the finely tuned, highly divergent preference functions are likely an additional source of LD between male and female traits in this system. Furthermore, pleiotropy of gene expression presents an underappreciated mechanism to link sexually dimorphic phenotypes.

Behavioral isolation is often one of the first and most potent forms of reproductive isolation to arise (Mayr 1963; Coyne and Orr 2004). This is somewhat paradoxical given that gene flow is often ongoing early in speciation (Kirkpatrick and Ravigne 2002; Bolnick and Fitzpatrick 2007; Nosil 2008) and behavioral isolation typically requires linkage disequilibrium (LD) between at least two loci to be maintained. Gene flow and subsequent recombination threaten to break down this LD eroding isolation (Pinho and Hey 2010). Accordingly, the genetic architecture of behavioral isolation is a key feature that may predict the likelihood of speciation.

Often, LD needs to be maintained among multiple loci: Females may select males (or vice versa) based on multiple traits, each with different genetic underpinnings (Candolin 2003; Bro‐Jorgensen 2010) and LD between these traits will increase divergence. Furthermore, both signal and preference are often polygenic (Bakker and Pomiankowski 1995; Ritchie and Phillips 1998; Gleason et al. 2002; Chenoweth and Blows 2006; Chenoweth and McGuigan 2010). The magnitude of LD required to maintain behavioral isolation in the face of gene flow is directly related to (1) the number of loci underlying signaling phenotypes (as well as the number of signaling phenotypes) and preferences, (2) their physical location in the genome, and (3) their effect sizes (Via and Hawthorne 1998; Coyne and Orr 2004; Arbuthnott 2009).

The theoretical literature has shown that speciation proceeds more readily if signal and preference loci are physically linked (Kirkpatrick and Hall 2004; Servedio 2009). Physical linkage may be attained via a single locus with pleiotropic effects, close genomic proximity of separate signal and preference loci, or self‐referencing, that is, trait‐preference matching (Servedio and Boughman 2017; Kopp et al. 2018). A further alternative is parental imprinting of mating preferences. In that case, signal‐preference inheritance essentially becomes a form of self‐referencing (Verzijden et al. 2005; Servedio et al. 2009). However, we have limited empirical insights into the genetic architecture of behavioral isolation. Current data provide evidence for linked signal and preference loci (pleiotropy or close linkage) in certain systems such as for color morphs and preferences in Heliconius butterflies (Kronforst et al. 2006; Merrill et al. 2011, 2018) and Medaka fish (Fukamachi et al. 2009), acoustic communication in Laupala crickets (Shaw and Lesnick 2009), pheromone signals and discrimination in Drosophila melanogaster (Marcillac et al. 2005), and morph and morph preference in Erythrura finches (Pryke 2010).

However, other systems show reduced linkage between trait and preference loci, such as cosmopolitan/Zimbabwe D. melanogaster (Ting et al. 2001) or completely unlinked signal and preference genes such as in moths (Dopman et al. 2004; Groot et al. 2009; Gould et al. 2010; Koutroumpa et al. 2016). When signal and preference loci are unlinked, gene flow will quickly break down LD. Theoretical results show that speciation is unlikely in that case (Servedio and Boughman 2017; Kopp et al. 2018) unless signals are “magic traits” that are locally adaptive (Dieckmann and Doebeli 1999; Kirkpatrick and Ravigne 2002; Doebeli 2005) or if preferences are “magic preferences” that are under divergent ecological selection (van Doorn et al. 2004; Weissing et al. 2011). Even after initial divergence, variation in unlinked or partially linked signals and preferences is constantly threatened by the homogenizing effects of gene flow (Servedio and Bürger 2018).

A key determinant of whether sexual selection aids or hinders divergence and thus the likelihood of speciation is the shape of the preference function (Servedio and Boughman 2017; Kopp et al. 2018). Unimodal preference functions centered around the same mean value across diverging populations (Fig. 1A) can lead to stabilizing sexual selection, which impedes divergence (van Doorn et al. 2004; Weissing et al. 2011; Kopp et al. 2018) or even lead to homogenization of preference loci across populations (Servedio and Burger 2014). If females from different populations prefer the same male traits, this leads to an increase in gene flow, which then homogenizes the two gene pools at the trait and preference loci. Strong physical linkage of preference alleles with signal alleles would mitigate these counterproductive effects during divergence and increase the likelihood of speciation (Servedio and Burger 2014). However, with open‐ended (Fig. 1B), relative preferences, or with strongly divergent unimodal preferences (Fig. 1C), sexual selection can facilitate divergence and speciation will proceed more readily, especially when signaling traits are also ecologically relevant (Kondrashov and Kondrashov 1999; Doebeli 2005). With these types of preference functions, strong physical linkage between preference and signal loci may not be necessary for divergence (Lande 1982). More empirical examples where we have knowledge of both the shape of preference functions as well as information about the genetic architecture of signals and preferences are needed to empirically test these theoretical results.

Figure 1

Schematic of male trait distribution and female preference function. (A) Unimodal, nondivergent preferences lead to stabilizing selection. (B) Open‐ended or (C) strongly divergent preference functions exert strong directional selection on the male trait, thereby creating genetic covariance even if loci reside on different (parts of) chromosomes.

Additional sources of variation in the genetic architecture may result from the fact that linkage will be detected more readily in systems where loci of relatively large phenotypic effect are linked among suites of traits or among signals and preferences. When more subtle aspects of the genetic architecture of coevolving traits are linked, important associations might be missed using standard quantitative trait locus (QTL) methods. Behavioral variation between closely related species in general, and variation in sexually dimorphic traits contributing to behavioral barriers specifically, is strongly influenced by regulatory variation (Williams and Carroll 2009; Etges 2014; Mack and Nachman 2017), but the genetic architecture of expression and its role in speciation are widely underappreciated empirically and theoretically (Mack and Nachman 2017). One intriguing but unexplored example is that gene expression variation, due to the ubiquitous pleiotropy of regulatory variants (Chesler et al. 2005; Gibson and Weir 2005), could be a powerful means of generating additional LD between signal and preference if these traits have shared regulatory pathways or otherwise coexpressed loci. Although traditional QTL analysis may uncover these variants, integrating pQTL (i.e., with a behavioral, morphological, or physiological trait as the response variable) and eQTL (i.e., with the expression level of a gene or transcript as the response variable) presents a potentially powerful approach to uncover additional loci of small effect if eQTL harbor mutations with weak effect on the phenotype but with sufficiently strong effects on gene expression related to that phenotype.

Here, we use QTL mapping to identify the number, distribution, and effect size of loci associated with variation in the multivariate acoustic mate signal and with a major dimension of sexual selection resulting from female preference for the song rhythm in the field crickets Gryllus rubens and G. texensis. These sibling species are widely distributed across the Eastern and Southern United States (Alexander 1962; 0057alker 2017). Acoustic mate choice is a major driver of reproductive isolation (Walker 1998; Gray and Cade 2000; Gray 2005; Blankers et al. 2015a), which is strong. Evidence shows that no natural hybrids have been documented and no females inseminated with heterospecific sperm have been collected (Gray and Cade 2000). Demographic analyses show that gene flow ceased roughly 18,000 years ago after initial divergence commenced 0.5 million years ago (Blankers et al. 2018b). Two male song traits that have diverged strongly between the species, pulse rate (i.e., the repetition rate of sound pulses) and carrier frequency (the pitch of the song), are both associated with unimodal preference functions. However, pulse rate preferences are finely tuned to the male song and strongly divergent among species, whereas carrier frequency preferences are broadly overlapping across species (Blankers et al. 2015a,b).

We hybridized wild‐caught parental lines in the lab to obtain segregating mapping populations and looked for associations between transcriptome‐wide single nucleotide polymorphism (SNP) markers and variation in pulse rate and carrier frequency as well as for pulse rate preferences (pQTL scan). We then correlated phenotypic variation in the mapping population with variation in gene expression across more than 27,312 transcripts and performed an eQTL scan for all trait‐associated transcripts. Our results on patterns of linkage among pQTL and eQTL significantly advance our understanding of the genetic architecture of behavioral isolation and provide important new insights into mechanisms of trait‐preference coevolution and divergence.

Material and Methods

Crickets were collected from allopatric locations (G. texensis: 84 females from Austin [TX], Lancaster [TX], and Round Rock [TX]; G. rubens: 76 females from Gainesville [FL], Lake City [FL], and Live Oak [FL]) but patterns of reproductive isolation are similar across zones of sympatry and allopatry indicating that reinforcement is absent in this system (Izzo and Gray 2004). We generated eight mapping families encompassing all four possible types of backcrosses to pure G. rubens (Fig. S1 and Supporting Information Methods) using parental individuals selected to maximize the potential phenotypic space for the hybrid offspring. Selected pairs were kept in the breeding boxes with water and food ad libitum and oviposition substrate for one week after the day the first eggs were recorded, after which both individuals were sacrificed and processed for RNA sequencing.

Males were recorded individually overnight in a dark, anechoic chamber. Digitization of and parameter estimation from individual male song recordings were done using custom software written by R.M.H in LabVIEW 2009 (National Instruments, Austin, TX). We retained pulse rate and carrier frequency for further analyses as these are the traits that differ most strongly between G. texensis and G. rubens (Blankers et al. 2015a,b). Recording temperature (25.1 ± 1.05°C [mean ± SD]) was used to standardize the measurements. Female preferences were tested under dark, anechoic conditions using a trackball system (Hennig et al. 2016): a Styrofoam sphere floating on pressurized air that can be easily moved by the cricket, while infrared sensors underneath record the sphere's movement in lateral and longitudinal directions. Custom software (written in LabVIEW 2009) was used to present stimuli (Table S1) as well as negative (silence and pure frequency tones) and positive (highly attractive stimulus) controls and to analyze the feedback from the optical sensors. The lateral movement of a female during signal presentation was averaged between the consecutive playbacks from the two speakers (order of active and silent speaker was randomized across trials) and normalized with respect to the response to the attractive control signal.

Individual preference functions can be deconstructed into separate components (Bailey 2008; Fowler‐Finn and Rodríguez 2012; Kilmer et al. 2017). The stimulus with the highest phonotactic response is the peak preference and quantifies the strongest preferred trait value. The degree to which males with nonpreferred trait values are discriminated against is the preference strength. Finally, the amount by which a male can differ from the most preferred trait value and still attract a female is measured by the tolerance, typically calculated as the width of the preference function at 67% of the maximum response. However, G. rubens and G. texensis differ strongly in peak preference, but only marginally in preference strength (Blankers et al. 2015b) and tolerance (unpublished results); only peak preference is sufficiently divergent for QTL mapping; In an alternative approach, we quantified preference functions more broadly by projecting individual responses of all backcrosses to all eight stimuli onto a linear discriminant function (obtained using “lda” in the R‐package “MASS”) (Venables and Ripley 2002), which had been trained on parental data (N = 73 G. rubens and N = 44 G. texensis females). This linear discriminant score will be referred to as pulse rate preference function from hereon, because it describes multiple aspects of interspecific variation in female preference through the variable correlation of all test patterns with the linear discriminant function (Table S2).

After phenotyping and/or crossing, each individual was played back its control stimulus for 10 min, preserved in RNAlater following the manufacturer's recommendations, and transferred to –80°C. All libraries were sequenced on a HiSeq 2000 (Illumina, San Diego, CA) at a depth of 13 libraries per lane with paired‐end 100 bp reads. Reads were processed using Flexbar (Dodt et al. 2012) and transcript‐level information was obtained by mapping the reads against the G. rubens reference transcriptome (Berdan et al. 2016) using Bowtie2 (Langmead and Salzberg 2012). SNPs were called using the Genome Analysis Toolkit (DePristo et al. 2011; Van der Auwera et al. 2013) and filtered using GATK and VCFtools (Danecek et al. 2011). Additional details on library preparation, SNP calling and filtering can be found in the Supporting Information Methods.

LINKAGE MAPPING

We conducted a chi‐square test for every SNP to determine if the segregation of alleles fit an autosomal or a sex‐linked model (false discovery rate [FDR] corrected P < 0.1). We removed loci if more than two families had missing genotype data and retained one SNP per transcript. Exceptions were made for eight loci that were of special interest. Because crickets have XX‐XO sex determination, only families with F1 hybrid dams have recombining sex chromosomes. We were able to recover sufficient X‐linked markers only in a single family of backcross type D (Fig. S1). All linkage and QTL mapping information for the X chromosome is thus based on that single family of 40.

Linkage maps were generated in JoinMap 4.1 (van Ooijen 2006) for each family individually. The total sample size was 288 (143 females and 145 males) and family sizes varied between 25 and 43. Linkage groups (LGs) were created by removing duplicate markers and then using a log‐of‐odds (LOD) threshold equal to 4.0 or 5.0 to generate LGs. The Kosambi mapping function was used to convert recombination frequencies to centi‐Morgans (cM). A consensus map was constructed using the map integration tool. Linkage groups from individual families were joined if they shared two or more markers. This map was used for the QTL analyses below.

We expanded this map to create a “dense” linkage map, which contained markers for as many transcripts as possible. This allowed us to find post hoc locations for loci that may be picked up in past or future scans. To do this, we merged genotypes across all families, combining both sexes, excluding markers with high levels of segregation distortion, and imported the genotype data in a new JoinMap file. With the original map as scaffolding, we used the regression algorithm and Kosambi mapping function to place as many markers as possible on the new map. This approach can bias QTL mapping for which the dense map was therefore not used, but gives information about the linkage of many (potentially interesting) genes to markers present in the map used for QTL mapping.

HERITABILITY AND GENETIC COVARIANCE

To estimate narrow‐sense heritability of and genetic covariance among male signal traits and female preferences, we used phenotypic data from grandparental and parental lines and their backcross offspring. We first fitted mixed models in lme4 (Bates et al. 2014) and estimated heritability using REML. We then fitted Bayesian Animal models in MCMCglmm (Hadfield 2010) using an inverse Wishart prior (Gelman 2006) and checked for autocorrelation, effective sample size, and chain convergence following the MCMCglmm course notes (Hadfield 2012). We then fitted multiresponse models with male pulse rate, male carrier frequency, and female pulse rate preference as response variables and ran 1 million iterations. The median and 95% Honest Posterior Density (HPD) interval of the heritability of each trait and genetic covariance (and correlation) between each trait pair were estimated from the posterior distribution, discarding the first 100,000 samples as burn‐in.

We used similar models to estimate the heritability for each of the 27,312 transcripts, except for these models we only ran 100,000 iterations to accommodate computational resources. Due to the asymptotic patterns of some of the posterior distributions (approximating but not overlapping zero), we considered all transcripts with the lower tail of the 95% HPD interval higher than 0.01 to have nonzero heritability.

pQTL MAPPING

The goal here was to establish the number and distribution of genetic loci contributing to variation in the main divergent phenotypes used in intersexual acoustic communication. We used R/QTL (Broman et al. 2003) in R (R Development Core Team 2016) to detect QTL for pulse rate, carrier frequency, and pulse rate peak preference and preference function (linear discriminant scores) at FDR <5% (“significant”) or FDR <63% (“suggestive”), following recommendations by the Complex Trait Consortium (Abiola et al. 2003). We excluded two males for which song recordings did not meet minimal quality standards leaving 143 females and 142 males for pQTL mapping. We first used “scanone” with Haley–Knott regression (Haley and Knott 1992) to identify the single strongest QTL for each trait, followed by 1000 permutations to establish a significance threshold at the 5% and 63% level. We then used the multiple‐QTL model approach (Broman and Sen 2009) to scan for additional QTL, refining QTL positions and establishing whether the model LOD score increased beyond the penalized LOD score threshold. The thresholds for FDR equal to 5% and 63% were obtained using 1000 permutations of the “scantwo” function. Cross type was included as a covariate in the models initially but removed if not significant. The magnitude of the additive effects and the 95% Bayesian credible interval was estimated from the model. To estimate the true number of loci underlying the phenotypic traits, we used a custom code, based on the study of Otto and Jones (2000), to estimate the QTL detection threshold, the true number of loci, and the amount of missing variation given the results of our experiment. The code is available at github.com/thomasblankers/statistics/QTL_power_detect.r.

eQTL MAPPING

The goals here were (1) to identify transcripts for which expression covaries with the main phenotypic traits used by males (pulse rate, carrier frequency) and females (pulse rate preference) in intersexual acoustic communication and (2) to unravel the genetic architecture of the expression of these transcripts. Reads from all backcross individuals were separately aligned to the reference transcriptome using Bowtie (Langmead et al. 2009) and transcript abundances were calculated for each stage using RSEM (Li and Dewey 2011). We imported the read abundance data into R using “tximport” (Soneson et al. 2015). We performed a differential expression analysis using a continuous model with both pulse rate and carrier frequency as fixed effects and cross as a covariate (expression ∼ cross + pulse rate + carrier frequency) to account for cross effects and for the correlation between traits (see Results). We fit these models in DESeq2 (Love et al. 2014) with Wald's test for significance. For pulse rate preference, we similarly fit a continuous DESeq model with cross as a covariate. For all models, we considered transcripts with a Benjamini–Hochberg (Benjamini and Hochberg 1995) corrected P‐value <0.01 to be significantly associated with the trait of interest.

We then performed a similar analysis using robust regression models in the R package limma (Ritchie et al. 2015). Here, we used log2‐TMM normalization to normalize the count data using the “calcNormFactors” and “cpm” function in edgeR. We then fitted linear models with robust regression and estimated empirical Bayes statistics for differential expression. All loci with adjusted P‐value <0.01 were considered significant.

For eQTL mapping, we kept only those loci that were significant in both the DESeq2 and the limma analysis and that had nonzero heritability. We retained two sets: one including all the above transcripts (“permissive” set) and the other containing those that have a relatively strong relationship with the trait (“conservative” set). The latter set consisted of transcripts that had partial η 2 values (“eta.square” function in the R package heplots; Fox et al. 2018) for their association with variation in the trait in the top 25% (for pulse rate η 2 > 0.13, for carrier frequency η 2 > 0.07, and for pulse rate preference η 2 > 0.13).

We used the “mqmscanall” function on a trait‐by‐trait basis to perform a multiple eQTL scan for each transcript. We included 16 cofactors in the analysis, one for each LG at the median marker. We obtained LOD thresholds corresponding to FDR <5% using 1000 permutations of the “mqmscanall” function to establish significance of eQTL. For all significant eQTL, we checked if they were cis (eQTL and corresponding transcript on the same LG) or trans (eQTL and corresponding transcript on different LGs) by comparing the eQTL location with the position of the transcript on the dense map.

To examine whether the transcripts that covaried in expression with male and female traits were also differentially expressed between species, we performed a differential expression analysis using the grandparents used to create the QTL mapping families as well as individuals previously sequenced (using similar methods as described above) for a population genetic study (Blankers et al. 2018b). Transcripts were considered differentially expressed if the adjusted P‐value of Wald's significance test was <0.01 and read count differed at least twofold (log2‐fold difference ≥1).

Results

PHENOTYPES

Phenotypes were unimodally distributed within species or cross line (Fig. 2, Table 1). Values for first generation interspecific hybrids were intermediate but biased toward the maternal parent for pulse rate (G. rubens dam: 55.27 pulses s−1, G. texensis dam: 61.96 pulses s−1; t 22 = –6.72; P < 0.0001) and carrier frequency (G. rubens dam: 4.82 kHz, G. texensis dam: 5.01 kHz; t 22 = –3.9146; P = 0.0004). Backcross distributions were also unimodal and intermediate between interspecific hybrids and G. rubens. All traits follow expectations for polygenic, additive inheritance. In addition to the preference measurements used in the downstream analyses (i.e., the peak preference and the discriminant function score), the preference function scores for pulse rate were unimodal in both species, F1 hybrids, and first‐generation backcrosses (Fig. S2).

Figure 2

Phenotypic distributions. (A) Schematic crossing design. Diploid Gryllus rubens (blue) and G. texensis (red) were crossed to obtain heterozygote first generation hybrid offspring in both cross directions. All possible combinations of hybrid‐Gryllus rubens were paired to create segregating backcross offspring. (B–E) Phenotypic distributions of parental (top panels), hybrid (middle panels), and backcross (bottom panels) offspring. Male pulse rate and carrier frequency are shown in (B) and (C); female preference is shown in (D) (pulse rate preference, i.e., linear discriminant scores representing composite phonotactic scores on all eight pulse rate test stimuli) and (E) (peak preference). The inset map in (E) shows the approximate geographic distribution of the parental species and their zone of overlap in the United States based on the study of Walker (2017).

Table 1

Phenotypic distributions of parental and hybrid generations

	Males	Pulse rate		Carrier frequency		Females	Peak pulse rate preference		Pulse rate preference function
	n	Mean	SD	Mean	SD	n	Mean	SD	Mean	SD
G. rubens	73	45.34	3.86	4.73	0.27	24	50.00	0.00	–5.77	0.71
G. texensis	44	66.88	5.40	5.18	0.22	17	68.79	7.99	5.77	1.31
F1 rubtex	22	55.27	3.96	4.82	0.27	14	61.60	7.87	1.93	4.17
F1 texrub	28	61.96	2.79	5.08	0.17	12	61.42	8.55	1.37	4.49
Backcross	142	51.45	6.71	4.91	0.33	143	60.12	8.93	0.56	4.53

Notes: Pulse rate in pulses per second; carrier frequency in kilo Hertz; peak preference in pulses per second; and pulse rate preference function in dimensionless units of correlation with the first discriminant function. Rubtex are F1 individuals with a Gryllus rubens dam and texrub are F1 individuals with a Gryllus texensis dam. Backcross is the mean of all four possible backcrosses to Gryllus rubens.

We placed a total of 330 markers on our genetic map (Table S3). The markers were grouped in 15 autosomal LGs, one more than the number of autosomes for G. rubens (Yoshimura 2005), and an X‐linked group with a total map distance of 254.4 cM, an average marker spacing of 0.81 cM, and a maximum marker spacing of 14.40 cM. Linkage groups varied in length from 0.99 cM to 41.5 cM (mean 17.3 cM). On our dense map, we were able to place a total of 1611 markers. Of these, we were able to determine the position of 1349 markers, the remaining 262 markers have a LG assigned but not a position.

REML heritabilities estimated from sire variance were 0.91, 0.51, and 0.61 for pulse rate, carrier frequency, and pulse rate preference, respectively. The Bayesian Animal models gave similar results with 95% HPD between 0.48 and 0.99, 0.49 and 0.74, and 0.27 and 0.99, respectively. Correlations among the traits were high: median correlations were 0.49 for corr(pulse rate, carrier frequency), 0.92 for corr(pulse rate, pulse rate preference), and 0.46 for corr(carrier frequency, pulse rate preference); 95% HPD intervals did not overlap with zero: 0.24–0.70; 0.58–0.99; 0.16–0.71. All genetic covariances were positive, indicating that an increase in one trait was associated with an increase in the other trait.

Using single interval mapping, we detected only a single significant QTL for each trait, except for peak preference for which both an autosomal and an X‐linked QTL were significant at α = 0.05 (Fig. S3). Because single QTL scans have limited power to detect small effect QTL, we added the significant QTL identified in single interval mapping to a multiple QTL model (MQM) and proceeded to scan for additional QTL at 5% (i.e., significant QTL) and 63% (i.e., suggestive QTL) FDR. In the final MQM (Fig. 3; Table 2), we identified one significant and four suggestive pQTL for pulse rate (all on autosomes), two significant autosomal and one suggestive X‐linked pQTL for carrier frequency, two significant and two suggestive autosomal pQTL for the pulse rate preference function, and one significant autosomal and one significant X‐linked as well as two suggestive autosomal QTL for peak pulse rate preference.

Figure 3

pQTL scan. For each of the four traits, the log‐of‐odds scores along the 16 linkage groups are shown by the intensity of blue hues. The scale is shown on the top right. 95% Bayesian confidence intervals for significant (solid) and suggestive (dashed) are shown as boxes projected onto the heatmap. Red arrows indicate pQTL explaining >10% of the backcross variance. pr, pulse rate. For single quantitative trait locus interval mapping, see Figure S3.

Table 2

pQTL effects

LG	pQTL location (cM)	Nearest marker	Log‐of‐odds	AA	AB	QTL effect (trait mean ± SE)	% Species difference	% Backcross variance
Pulse rate
1	39.4	c214087_g2_i1	3.21	48.72	51.83	3.62 ± 0.94‡	17.24	14.8
3	6	c215368_g2_i3	1.04	49.85	50.70	2.11 ± 0.98*	10.03	4.6
5	25	c218669_g2_i3	0.97	49.14	51.32	2.15 ± 1.05*	10.23	4.3
10	0	c186619_g1_i1	1.14	49.31	51.39	2.33 ± 1.06*	11.09	5.1
12	0	c204487_g2_i1	1.22	49.11	51.23	2.18 ± 0.94*	10.37	5.5
cross			2.99			1.47 ± 0.39‡	37.92	13.8
Carrier frequency
1	39	c214087_g2_i1	1.93	4.84	4.98	0.16 ± 0.05†	34.89	8.9
3	1	c203593_g1_i1	2.89	4.83	5.00	0.19 ± 0.05‡	41.76	13.4
X	9.9	c205832_g2_i1	0.90	4.86	4.98	0.14 ± 0.06*	30.20	4.1
Pulse rate preference function
2	4.4	c142606_g1_i1	1.63	−0.38	1.32	1.84 ± 0.68†	15.94	7.4
3	0.1	c218168_g2_i2	2.28	−0.65	1.31	2.19 ± 0.67†	18.96	10.5
5	2	c217193_g1_i1	2.35	−0.49	1.37	2.20 ± 0.67†	19.07	10.8
12	2	c203868_g1_i1	1.91	−0.6	1.69	2.08 ± 0.70†	18.03	8.7
cross			2.7			0.53 ± 0.15‡	4.57	12.4
Peak pulse rate preference
2	5.1	c214277_g1_i1	1.94	58.29	61.56	3.99 ± 1.35†	21.23	8.8
3	0.1	c218168_g2_i2	1.35	58.25	61.20	3.30 ± 1.36*	17.55	6.1
5	0	c212100_g1_i1	3.00	57.52	62.17	5.07 ± 1.36‡	26.96	13.9
X	9.8	c205832_g2_i1	1.08	57.35	62.56	3.94 ± 1.61*	20.99	4.8
cross			2.17			0.94 ± 0.30‡	5.02	9.9

Notes: For each trait, the linkage group, the location, the nearest marker, the log‐of‐odds (LOD) score, the genotypic effects, and the pQTL effects (expressed in trait mean change and number of standard deviations in Gryllus rubens) are given, and percentage of species difference and backcross variance explained of each of the pQTL effects is shown. Significant pQTL (<5% false discovery rate based on penalized LOD score improvement of the multiple QTL model) are in bold. All pQTL effects are significantly larger than zero: * P < 0.05; † P < 0.01; ‡ P < 0.0.01.

The QTL for peak preference and preference function were largely overlapping, excepting an X‐linked QTL for peak preference and a suggestive QTL on LG 12 for preference function (Fig. 3, Fig S3). Carrier frequency and pulse rate also mapped to similar regions, with significant QTL overlapping on LG 1, and a suggestive pulse rate QTL on LG 3 overlapping with a significant QTL for carrier frequency. There was also QTL colocalization between male and female traits on LG 3 (all traits, QTL for pulse rate is suggestive), LG 5 (suggestive QTL for pulse rate, significant QTL for preference), LG 12 (suggestive QTL for both pulse rate and preference), and the X chromosome (carrier frequency and peak pulse rate preference; Fig. 3, Fig. S3).

The effect sizes for each QTL are shown in Table 2. For pulse rate, haploid allelic effects from five loci explained a total of 12.39 pulses s−1, or 34.3% of the backcross variance. For carrier frequency, this was 0.44 kHz or 26.4% of the variance across three loci, and for peak preference and pulse rate preference function, the total of four QTL effects was 16.07 pulses s−1 and 7.81 or 37.4% and 33.6% of the backcross variance, respectively. The combined effect size expressed as percentage of the difference between parental mean phenotypic values is much larger, but we note that these estimates are biased upward due to our selective breeding of individuals from the extremes of the distributions. All QTL effects were significant (P‐value for one sample t‐test <0.05) and of the same sign (i.e., G. texensis alleles always increased the trait values). Cross type effects were significant for all traits except carrier frequency but are not included in the sum of haploid allelic effects.

Using equation 6 in Otto and Jones (2000), we estimated the true number of loci (95% confidence interval) to be 23.30 (8.36–50.1) for pulse rate, 7.25 (2.61–15.59) for carrier frequency, 9.08 (2.83–21.11) for pulse rate peak preference, and 18.22 (5.66–42.32) for pulse rate preference function.

OVERLAP OF OUTLIERS WITH pQTL

Because we have access to the full (annotated) transcriptome of each individual, we asked whether pQTL regions contain candidate loci for behavioral isolation between G. rubens and G. texensis. Candidate loci should (1) be located within a QTL region and (2) show signatures of divergence. We scanned for loci that met these criteria using previously published data. In Blankers et al (2018b), we identified 231 loci that were potentially under selection using a combination of population genetic statistics (D and Tajima's D). Of these 231 loci, we were able to place 122 on our dense map (Table S4). Nine of these loci (six for carrier frequency, five for pulse rate, and four for female preference; note that some of the QTLs overlap) were within 2.5 cM of a QTL and are promising candidate genes (see Discussion). A cutoff at 2.5 cM from the peak means that LD between QTL peaks and any potentially interesting genes decays at less than ∼2.5% per generation. We acknowledge that other thresholds could be similarly rationalized and would affect the results somewhat.

We identified 430, 35, and 26 transcripts for which expression covaried with pulse rate, carrier frequency, and pulse rate preference, respectively (Tables S5–S7). The average narrow‐sense heritability, , of these transcripts was 0.43 ( = 0.37 for all 27,312 transcripts) with a minimum of 0.00. After removing 69, two, and one transcript(s) with very low heritability (lower 95% HPD interval <0.01), = 0.49 with a minimum heritability of 0.07. The partial η 2 for trait variation explained varied between 0.05 and 0.23 when considering all transcripts with nonzero heritability (“permissive” set) and between 0.12 and 0.23 when considering only the transcripts in the top 25% for the magnitude of trait association (“conservative” set; 109, 10, and eight transcripts for pulse rate, carrier frequency, and pulse rate preference, respectively). Some of these transcripts (43 out of 430 transcripts for pulse rate, 12 out of 35 for carrier frequency, and five out of 26 for preference) were also differentially expressed between the pure species (Tables S5–S7, Fig. S4).

eQTL were significant between LOD >3.0 and LOD >2.0 depending on the trait and set of transcripts. We detected a total of 56 significant eQTL, 15 of which were from the conservative set of trait‐associated transcripts, the remaining 41 from the permissive set. Of these, 39 from the permissive (six from the conservative) transcripts covaried with pulse rate, six (four) with carrier frequency, and 11 (five) with pulse rate preference (Fig. 4, Table S8). The majority of the transcripts for which an eQTL was identified did not have a LG assigned and therefore information about cis versus trans regulation is limited; however, all seven eQTL for transcripts of known location were cis regulated (Table S8).

Figure 4

eQTL scan. The log‐of‐odds (LOD) score traces from the multiple quantitative trait locus (QTL) model (one cofactor per linkage group) are shown for all transcripts that significantly covaried with pulse rate (top panels), carrier frequency (middle panels), and pulse rate preference (bottom panels). The horizontal black solid line shows the LOD threshold above which the false discovery rate is below 5%. All significant eQTL are shown in red. The heatmap insets below each panel show the pQTL results (Figure 3) for comparison. See also Figures S5–S7 for transcript specific eQTL scan results.

We find significant eQTL on most LGs, except LG 6, 10, and 13 (Fig. 4, Table S8). There is a trend for the full set of eQTL to colocalize as 40/56 eQTL colocalized with at least one other eQTL (Fig. 4, Table S8). This trend was apparent in the conservative set when comparing trait‐specific eQTL but not when comparing eQTL among traits (Fig. 4, Table S8). In the permissive set of trait‐associated transcripts, we observe much more extensive colocalization both within and between traits as well as between male pulse rate and female pulse rate preference (most notably on LG 5). Some of the eQTL locations also correspond to or are closely linked to pQTL locations discussed in the previous section (pulse rate: LG 1, LG 3, LG 5, and LG 12; carrier frequency: LG 3; pulse rate preference: LG 3 and 5; Fig. 4). However, there are also LGs that have eQTL for transcripts associated with a trait for which there was no pQTL on the LG.

Discussion

Behavioral barriers to gene flow arise through differentiation in male and female mating communication traits and are a powerful mechanism to promote divergence in the earliest stages of speciation (Coyne and Orr 2004). To better understand this process, we need to examine how the genetic architecture of signaling and preference traits as well as the shape of preference functions and their effect on signal distributions contribute to generating LD between coevolving male and female traits. In this study, we jointly examined the genetic architecture of male signal traits (pulse rate and carrier frequency of the song) and female preferences in two species of North American field crickets G. rubens and G. texensis. These species diverged ∼0.5 million years ago followed by a long period of bidirectional gene flow that lasted until ∼18,000 years ago (Blankers et al. 2018b). Preference functions for pulse rate closely track male song distributions, and both male and female traits have diverged conspicuously between the species (Gray and Cade 2000; Blankers et al. 2015a,b).

Our results reveal physical linkage between the two coevolving song traits, pulse rate and carrier frequency, as well as between coevolving male pulse rate and female pulse rate preference. However, the pQTL of largest effect was never shared between any two traits. We extended our analysis of the genetic architecture into the regulatory pathways that potentially underlie the behavioral traits of interest. We observed tight linkage of eQTL for multiple transcripts associated with the same trait as well as for transcripts associated with different (male and female) traits. This intriguing result suggests linked regulatory variation may contribute to coevolution of song and preference. Thus, there are multiple dimensions by which physical linkage may contribute to maintaining LD between signals and preferences. However, because physical linkage is incomplete, the striking coevolution of male and female traits is likely aided by sexual selection resulting from the shape of the pulse rate preference function in relation to the male signal distribution. We hypothesize that these mechanisms jointly facilitate trait‐preference coevolution and the maintenance of a strong prezygotic reproductive barrier despite gene flow.

INTEGRATED SONG SIGNALS

We showed strong, positive genetic covariance between two male song traits, pulse rate and carrier frequency, that are known to be strongly correlated phenotypically (Blankers et al. 2015b, 2017). The strong covariance observed here would allow for a correlated response to selection. Although we also report pQTL that are unique to only one trait, the overlapping QTL on LG 1 and LG 3 have relatively high effect sizes (>10% of difference between species and >5% of the backcross variance), suggesting that phenotypic effects of linkage may be substantial. This linkage may have resulted in indirect selection on carrier frequency due to strong selection on pulse rate. This process would result in the coevolutionary patterns observed for carrier frequency and pulse rate across closely related Gryllus species, despite broadly overlapping preference functions for carrier frequency (Blankers et al. 2015a; Hennig et al. 2016). Physical linkage between loci underlying the traits of an integrated sexual signal potentially facilitated signal divergence in multiple dimensions (pulse rate and carrier frequency) even though preference has diverged only in one dimension (pulse rate).

INTEGRATED FEATURES OF PULSE RATE PREFERENCE

Aspects of female preference are also tightly linked and are likely to cosegregate. It may seem trivial that pQTL scans for pulse rate peak preference and pulse rate preference function (i.e., linear discriminant scores) are concordant and that focus should be on the discordance instead of the similarity. However, the two measures incorporate different aspects of mate choice. The peak preference score is determined solely by the stimulus eliciting the strongest phonotactic response. This is what in theoretical literature of sexual selection and mate choice behavior is generally considered “preference” (Edward 2015). The linear discriminant function captures multiple aspects of the preference function shape (e.g., peak, width, and skew) by including responses to stimuli in heterospecific ranges and outside of the ranges of either species (Fig. S2, Table S2). Thereby, this measure is a composite representation of the interspecies differences in the preference function. The fact that the pQTL scans associated with these distinct measures of preference gave qualitatively similar results, differing only in the magnitude of correlation between genotypes and phenotypes and the presence of a small‐effect X‐linked QTL, suggests that the genetics of peak preference and preference to all tested stimuli are highly integrated. This provides rare empirical evidence for the idea that the genetic underpinnings of difference aspects of mate preference (e.g., peak preference and choosiness, responsiveness) cannot be straightforwardly separated (Kopp et al. 2018). However, we acknowledge peak preference contributes substantially to the first discriminant function indicating nonindependence of the two measures. We were unable to directly test the genetic architecture of preference strength, tolerance, or responsiveness separately because these properties have diverged only marginally between G. rubens and G. texensis.

SIGNAL‐PREFERENCE COEVOLUTION

Genetic covariance between signal and preference is expected if traits coevolve within populations and necessary for sexual selection to drive phenotypic divergence (Fisher 1930; Lande 1981; Kirkpatrick 1982; Kirkpatrick and Hall 2004). Empirically, it is unclear whether the dominant mechanism of genetic covariance is physical linkage (proximate loci or pleiotropy) or directional mate preference. Theoretically, both mechanisms would lead to LD between signal and preference and accentuate effects from directional selection (Andersson and Simmons 2006), but LD without physical linkage has been shown to be sensitive to gene flow between partially isolated populations (Servedio and Boughman 2017; Kopp et al. 2018). However, the extent to which physical linkage is required to maintain LD between signals and preferences is also sensitive to the shape of female preferences: for example, speciation proceeds more readily with open‐ended or relative preferences or with strongly divergent unimodal preferences (Kondrashov and Kondrashov 1999; Doebeli 2005)

We show that in G. rubens and G. texensis, for which detailed demographic analysis has demonstrated divergence in the face of (primary) gene flow, there is some physical linkage between song and preference loci, but also pQTL unique to each trait. Linkage was always between a significant pQTL and a suggestive pQTL or between two suggestive pQTL, which are of comparably weak phenotypic effect and associated with more statistical uncertainty. Compared to previous examples in crickets (Shaw and Lesnick 2009), flies (Marcillac et al. 2005), lepidopterans (Kronforst et al. 2006), and fish (Fukamachi et al. 2009), we observe a lesser degree of linkage, showing that divergence in male and female traits does not always require physical linkage (e.g., see Ting et al. 2001; Ritchie et al. 2005; Smadja and Butlin 2009). However, we observe equally strong or stronger levels of genetic covariance between song and preference. We suggest that this is partly explained by the shape of female preference. Although learned preferences could similarly generate LD between signals and preferences (Verzijden et al. 2005; Servedio et al. 2009), parental imprinting in Gryllus is unlikely because generations are nonoverlapping and crickets can typically only learn from siblings or unrelated individuals. In that case, LD would be weakened rather than strengthened (Servedio et al. 2009; Verzijden et al. 2012). If the shape of the preference function relative to the population distribution of the signal results in directional selection on the signal, sexual selection can generate strong covariance between traits and preferences without the need for tight physical linkage. In our system, pulse rate preferences are nonoverlapping, unimodal, and sharply tuned to the male song distribution (Blankers et al. 2015a,b). Comparisons across related Gryllus species (Blankers et al. 2015a; Hennig et al. 2016) suggest small differences in the preference result in strong selection on the signal. Depending on the ancestral distributions of the male trait, this may have been enough to drive divergence of signal and preference without strong physical linkage.

Comparing the location of the pQTL to the annotated transcripts on the dense linkage map, we were able to identify nine candidate genes that may underlie these traits. A promising candidate is nervana2, a sodium/potassium‐exchanging ATPase vital for hearing in D. melanogaster (Roy et al. 2013), which is located near a QTL for both carrier frequency and pulse rate on LG1. We also found several candidates associated with song production: Elongin B, a subunit of the elongin complex, is found near a QTL for pulse rate on LG 10 and osa is found at the QTL peak for carrier frequency on LG 3. In D. melanogaster, both osa and the elongin complex contribute to wing vein patterning (Terriente‐Félix and de Celis 2009; Rougeot et al. 2013). Crickets produce song by rubbing specialized structures on the forewings together and wing vein patterns control the resonant properties and thus the resulting sound (Bennet‐Clark and Ewing 1968; Nocke 1971; Bennet‐Clark 2003). Previous work has found that forewing shape varies between G. rubens and G. texensis and covaries with spectral and temporal song variation, although it is unclear how song rhythm would be affected by wing shape (Blankers et al. 2018). Our results further highlight two clock genes, period (about 5 cM away from the pulse rate and female preference QTLs on LG 12) and timeless (directly under the QTL on LG 1 for carrier frequency and pulse rate; however, tim is not a population genetics outlier). These clock genes regulate circadian rhythms, including the timing of mating behavior in fruit flies and crickets (Sakai and Ishida 2001; Fergus and Shaw 2013). Additionally, clock genes, specifically period, have been implicated in regulating biological processes acting on shorter time scales, including interpulse intervals in vibrational signals in D. melanogaster (Medina et al. 2015) and rhythmic fluctuations of pulse intervals of the courtship song (Kyriacou and Hall 1980). However, the role of circadian rhythm genes in the fly's song rhythm is heavily debated (Stern 2014; Kyriacou et al. 2017; Stern et al. 2017). This small set of candidate genes should be further investigated to determine their functional role in Gryllus crickets.

We acknowledge there are some caveats to the results discussed here. Statistical (i.e., the Beavis effect) (Beavis 1998) and experimental (selective breeding to optimize phenotypic space in backcross generations) considerations cause our results to be somewhat biased toward loci of large effect. This may either obscure additional linkage among (small effect) QTL or overestimate the total amount of linkage. Given the phenotypic distances in these closely related species, the sample sizes were not sufficient to detect smaller effect pQTL (<10%) of which there are likely plenty (e.g., Shaw et al. 2007; Blankers et al. 2018a). The fact that only carrier frequency and not pulse rate has X‐linked QTL is particularly puzzling, especially in the light of strong signatures of X‐linkage for pulse rate in reciprocal interspecific hybrids. This likely reflects the difficulties we had in reconstructing linkage on the X‐chromosome because only few markers segregated following expectations for XX‐XO mating systems. Additionally, limitations in power to detect pQTL (due to limited phenotypic divergence) make it difficult to distinguish a single pleiotropic locus from multiple loci in close genomic proximity. However, pQTL and eQTL results consistently point toward a mixture of linked and unlinked loci for the coevolving male and female traits, suggesting that it is unlikely that these caveats have falsely led us to reject completely linked or completely independent segregation of loci.

eQTL OVERLAP WITH pQTL

We found strong overlap between pQTL and eQTL, supporting a central role for regulatory variation in behavioral evolution and reproductive isolation (Wray 2007). In some cases, pQTL and eQTL peaks map to proximate or even identical locations. The genes at these overlapping pQTL/eQTL are good candidates for controlling the traits in question, although unfortunately we were not able to annotate any of these. In other cases, we detect more distantly located loci, as well as eQTL on LGs with no pQTL and vice versa. One reason pQTL and eQTL might be linked is because both detect a single regulatory variant: many trait‐associated SNPs in QTL and genome‐wide association studies are regulatory variants rather than protein coding variants (Nicolae et al. 2010) and this is particularly likely for behavioral and sexually dimorphic traits (Wray 2007; Williams and Carroll 2009). For example, small changes in the balance of excitation and inhibition within the neuronal recognition network can rapidly change the phenotype of female preference in crickets and katydids (Hennig et al. 2014). An alternative is that linked pQTL and eQTL represent tightly linked regulatory and coding variants. With the current data, we cannot distinguish between these alternative explanations. The eQTL that did not overlap with pQTL potentially represent false positives, in which case these transcripts are not related to the phenotype in question. However, an alternative hypothesis is that they are true QTL that were not picked up by our pQTL scan, because these loci have only small phenotypic effects but nevertheless sufficiently strong effects on trait‐associated gene expression variation to be picked up in the eQTL scan. More data from across a variety of taxa are needed to test this hypothesis.

PLEIOTROPIC GENE EXPRESSION AND SIGNAL‐PREFERENCE COEVOLUTION

The pleiotropic nature of gene expression is well‐known as many eQTL detected in transcriptome‐wide studies are concentrated in narrow genomic regions (Chesler et al. 2005; Gibson and Weir 2005; Hubner et al. 2005). We detect multiple eQTL for trait‐associated transcripts at identical or proximate locations, although most of the colocalization is observed only when the more permissive set of trait‐associated transcripts is considered (i.e., all heritable transcripts, including those with lower magnitudes of trait covariation). The most striking colocalization events occur on LG 3, where we detected pQTL and eQTL for transcripts associated with both male song traits and female song preference, and LG 5 where we mapped loci controlling expression of multiple pulse rate and pulse rate preference associated transcripts. We suggest that linkage of regulatory variants reflects an underappreciated genetic mechanism that can affect LD between signals and preferences. There is limited theory explaining the effects of regulatory variation on the efficacy of sexual selection in the face of gene flow. Existing theory generally indicates that regulatory variation can enhance the effectiveness of assortative mating (Ten Tusscher and Hogeweg 2009) and that linked cis regulatory loci can enhance the evolution of sex‐biased gene expression (Williams and Carroll 2009) and sexual dimorphism (Connallon and Clark 2010). Our findings provide important empirical insight into the potential for physical linkage, shared regulatory variation, and mate preferences to reciprocally shape LD during divergence with gene flow. We suggest that behavioral, quantitative genetic, and gene expression data be more broadly integrated to understand the effects from sexual selection on diversity across different biogeographic contexts of speciation.

Associate Editor: G. G. Rosenthal

Handling Editor: Mohamed. A. F. Noor

Supporting Information

Click here for additional data file.

Figure S1. Breeding scheme for the QTL.

Figure S2. Preference functions for pulse rate stimuli played back on the trackball system.

Figure S3. LOD score profiles of marker associations with all four traits for single QTL interval mapping.

Figure S4. Heatmaps for between species differential expression.

Figures S5‐S7. Heatmaps for LOD scores along the genome for each transcript associated with pulse rate, carrier frequency, and pulse rate preference, respectively.

Click here for additional data file.

Table S1. Typical stimulus array for G. rubens female to test pulse rate preference.

Table S2. Correlation coefficients for each of the eight test stimuli in the pulse rate test (Table S1) with their corresponding pulse rates.

Table S3. Genetic map for G. rubens × G. texensis.

Table S4. Potential outlier loci locations.

Table S5. Transcripts with expression levels correlated to pulse rate.

Table S6. Transcripts with expression levels correlated to carrier frequency.

Table S7. Transcripts with expression levels correlated to pulse rate preference.

Table S8. eQTL locations.

Table S9. Dense linkage map.

Click here for additional data file.