Differential introgression suggests candidate beneficial and barrier loci between two parapatric subspecies of Pearson's horseshoe bat Rhinolophus pearsoni.

Observations that rates of introgression between taxa can vary across loci are increasingly common. Here, we test for differential locus-wise introgression in 2 parapatric subspecies of Pearson's horseshoe bat (Rhinolophus pearsoni chinensis and R. p. pearsoni). To efficiently identify putative speciation genes and/or beneficial genes in our current system, we used a candidate gene approach by including loci from X chromosome that are suggested to be more likely involved in reproductive isolation in other organisms and loci underlying hearing that have been suggested to spread across the hybrid zone in another congeneric species. Phylogenetic and coalescent analyses were performed at 2 X-linked, 4 hearing genes, as well as 2 other autosomal loci individually. Likelihood ratio tests could not reject the model of zero gene flow at 2 X-linked and 2 autosomal genes. In contrast, gene flow was supported at 3 of 4 hearing genes. While this introgression could be adaptive, we cannot rule out stochastic processes. Our results highlight the utility of the candidate gene approach in searching for speciation genes and/or beneficial genes across the species boundary in natural populations.

Introgression describes the exchange of alleles via hybridization and backcrossing among distinct taxa (Anderson and Hubricht 1938; Anderson 1949), including species (Larson et al. 2014) and subspecies (Hird and Sullivan 2009; Sullivan et al. 2014), which can lead to incongruence among gene genealogies and produce shared and/or closely related alleles between distinct taxa. Similar genealogical patterns can also result from incomplete lineage sorting of ancestral polymorphism and it is very difficult to distinguish these 2 processes from each other (Wendel and Doyle 1998). Recently, coalescent isolation-with-migration (IM) models (Hey and Nielsen 2004) has been developed and proved to be successful in resolving this issue (e.g., Mao et al. 2010; Morgan et al. 2010; Choleva et al. 2014).

Growing numbers of studies have shown that rates of genetic introgression between taxa can vary across different loci (Dopman et al; Kronforst et al; Maroja et al; Ohshima and Yoshizawa 2010; Baldassarre et al; Taylor et al). Such differential introgression has commonly been attributed to the effects of natural selection or genetic drift, but may also reflect variable recombination rates due to structural features such as chromosome inversions (reviewed in Nachman and Payseur 2012) , or differences in the extent of linkage between genetic markers and genes under divergent selection (Payseur et al. 2004). As a consequence, it has been proposed that where introgressive hybridization occurs, some loci may be able to flow freely between taxa whereas others will be less able to do so, leading to a semipermeable genome (Wu 2001).

Loci that resist introgression are often considered to be putative “speciation genes” (Wu 2001; Dopman et al. 2005; Noor and Feder 2006), defined as genes “whose divergence made a significant contribution to the evolution of reproductive isolation between populations” (Nosil and Schluter 2011). In contrast to speciation genes, loci that introgress more readily may be selectively neutral or alternatively may confer fitness benefits in their new populations, although cases of so-called adaptive introgression are rare in wild animal populations (but see Fitzpatrick et al. 2009; Song et al. 2011; Pardo-Dieaz et al. 2012). Efforts to identify speciation genes and/or beneficial genes in natural populations of non-model organisms have frequently taken a differential introgression approach, whereby speciation genes show reduced levels of introgression, whereas beneficial genes show increased levels of introgression (Payseur 2010).

Hybridizing taxa offer promising systems in which to search for putative speciation genes and/or beneficial genes using the differential introgression approach (Payseur 2010; Larson et al. 2014). The Pearson′s horseshoe bat Rhinolophus pearsoni is such a system in which hybridization and introgression have been previously documented between its 2 subspecies (R. p. pearsoni and R. p. chinensis) (Mao et al. 2010). R. p. chinensis is restricted to eastern China, whereas R. p. pearsoni has a wider distribution across the southeast Asia including north of India, Nepal, Myanmar, west of China, Vietnam, Lao P.D.R, and Thailand (Csorba et al. 2003); currently these 2 subspecies are parapatric in central China (see Figure 1). Introgression of mtDNA had been detected from R. p. chinensis to R. p. pearsoni, whereas an analysis of 2 nuclear genes did not support the evidence for nuclear introgression between them (Mao et al. 2010). One explanation for this discordant introgression pattern between mitochondrial and nuclear markers is that those 2 nuclear markers examined previously were directly related to reproductive isolation or linked to loci relating to reproductive isolation. However, in Mao et al. (2010), no formal tests (e.g., likelihood ratio tests) were done to compare the fit of models with and without gene flow at those 2 loci. Thus, patterns of nuclear introgression between these 2 focal subspecies are still not clear.

Figure 1.

Map showing the sample sites of R. p. chinensis and R. p. pearsoni modified from Mao et al. (2010). Populations are presented as circles in which individuals are colored based on the taxon membership (R. p. chinensis: orange; R. p. pearsoni: blue).

To efficiently identify putative speciation genes and/or beneficial genes in our current system, we used a candidate gene approach by including loci that more likely contribute to reproductive isolation or undergo adaptive introgression on the basis of studies in other organisms. Studies from both laboratory crosses and wild populations have revealed that loci on the X chromosome contribute disproportionately to reproductive isolation (reviewed in Coyne and Orr 2004) and often exhibited reduced gene flow across the hybrid zone comparing with autosomal loci (Macholan et al. 2007; Geraldes et al. 2008), possibly due to Haldane’s rule (see Orr 1997; Payseur et al. 2004). Here, 2 X-linked markers were chosen as candidate loci possibly involved in reproductive isolation between our 2 focal subspecies.

For candidate beneficial genes, we chose hearing genes because hearing genes play important roles in echolocation, and echolocation in bats has been proved to evolve specifically for the detection of flying insects, which is essential for bats to survive and adapt to different environments (Jones and Holderied 2007). Currently, at least 4 hearing genes (e.g., FoxP2, Prestin, Tmc1, and Kcnq4) have been documented in echolocating taxa (Li et al. 2007, 2008; Davies et al. 2011; Liu et al. 2012). Three of them (FoxP2, Prestin, and Kcnq4) have been recently analyzed in another congeneric species Rhinolophus affinis and that study revealed gene flow at Prestin across the hybrid zone between 2 subspecies of R. affinis (Mao et al. 2014). Here, we tested for differential locus-wise introgression in our current system by performing phylogenetic and coalescent analyses on 2 autosomal genes (Chd1 and Sws1) taken from our previous study (Mao et al. 2010), 2 X-linked genes and 4 hearing genes. Among these loci, Sws1 was suggested to be pseudogenized in Rhinolophus (Zhao et al. 2009), thus, it can be used as a neutral control in the comparison with other loci. We predicted that hearing genes might exhibit higher levels of introgression because of their benefits for animals in adapting to new environments, whereas introgression rates of X-linked genes would be reduced, compared to other autosomal genes.

Materials and Methods

Ethics statement and sampling

All tissues used in this study were sampled from bats for our former project (Mao et al. 2010) (see details in Figure 1 and Table 1). The non-lethal procedure of sampling consisted of taking wing membrane biopsies from bats, and was approved by the National Animal Research Authority, East China Normal University (approval ID 20080209). Bats were initially assigned to R. p. pearsoni or R. p. chinensis on the basis of taxon-specific and non-overlapping call frequencies (Zhang et al. 2009; Mao et al. 2010). One congeneric species R. affinis was included as an outgroup.

Table 1.

GenBank accessions for all samples used in the molecular analysis. N means the location number as shown in Figure 1

N	Sample locations	Coordinates	Code	Prestin	Tmc1	FoxP2	Kcnq4	Usp9x	Pola1
	R. p. chinensis
1	Qingyang	N30:20:511 E117:50:128	AH	JX502283	KC874587, 93	JX502243	KC874518, 20	JX502378 -79	JX502319, 20,32,33
2	Jingxian	N30:26:785 E118:24:783	AH	JX502282	KC874583, 84,86	JX502244	KC874512, 21	JX502374 -75	JX502322, 28
3	Huangshanjinjiao	N29:45:107 E118:23:171	AH	JX502284	KC874603	JX502245	KC874511	JX502392 -93	JX502329 -31
4	Huangshanxinming	N30:23:181 E118:14:116	AH	JX502285	KC874585, 88	JX502246	KC874514, 28	JX502397, 99	JX502334
5	Fuchunsanling	N29:22:112 E117:34:324	JX	–	KC874589- 90	JX502247	KC874519	JX502427 -28	JX502321, 36
6	Fuchunqingfeng	N29:22:262 E117:39:357	JX	JX502291	KC874598- 600	JX502248	KC874524, 29-30	JX502422 -23	JX502325 -26
7	Fuchunqinhui	N29:22:662 E117:32:335	JX	JX502294	KC874602	JX502249	KC874523	JX502424, 26	JX502327
8	Guwang cave	N27:42:664 E117:41:531	FJ	JX502292	KC874591- 92	JX502251	KC874513,22	JX502408, 12	JX502340 -41
9	Yanzijiao	N27:48:511 E117:42:505	FJ	JX502290	KC874594, 96	JX502253	KC874515, 16,33	JX502405 -07	JX502338
10	Taining	N26:42:236 E117:29:867	FJ	JX502286-89	KC874595	JX502252, 54	KC874531 -32	JX502403 -04	JX502337
12	Liancheng	N25:12:404 E117:15:066	FJ	–	KC874606	–	KC874517	JX502402	JX502339
	Total in chinensis			10	18	12	18	22	20
	R. p. pearsoni
14	Zhangjiajie	N29:21:410 E110:34:783	HN	JX502300-01	KC874610, 11,16-18	JX502255	KC874546- 48,51-52	JX502440 -42	JX502359
15	Longshan	N29:12:7865 E109:18:454	HN	JX502299	KC874615	–	KC874538- 39	JX502437	JX502352
16	Yongshun	N29:03:720 E109:38:358	HN	JX502303	KC874608- 09	–	KC874540- 41	JX502438	JX502357
17	Jishou	N28:18:208 E109:39:175	HN	JX502302	KC874612- 13	JX502257	KC874542- 45	JX502445, 46,48	JX502355 -56
18	Fenghuang	N27:59:580 E109:33:786	HN	JX502304	KC874614	JX502256	KC874549- 50	JX502439	JX502358
19	Tianquan	N30:10:671 E102:75:831	SC	JX502295	KC874635	–	KC874555- 56	JX502483	JX502343
20	Baoxing	N30:54:779 E102:65:219	SC	JX502296-97	KC874637- 41	JX502264 -65	KC874557- 58,62-63	JX502485- 87	JX502342, 44-45
21	Emei	N29:34:803 E103:24:708	SC	JX502298	KC874636	JX502266 -67	KC874559 -60	JX502484	JX502346
22	Zhenfeng	N25:27:807 E105:29:977	GZ	–	KC874622- 23	JX502263	KC874537	JX502449	JX502353
23	Anlong	N25:16:577 E105:31:931	GZ	JX502305	KC874620- 21	JX502260 -62	KC874534- 36	JX502450- 52	JX502348 -51
24	Xingyi	N25:04:374 E104:53:067	GZ	–	–	JX502259	–	JX502456	JX502354
25	Jingchengjiang	N24:70:924 E108:15:947	GX	–	–	JX502258	–	JX502457	–
26	Wuming	N23:43:161 E108:37:979	GX	–	KC874619	–	KC874553- 54	JX502458	JX502347
27	Meizi	N22:98:356 E103:68:761	YN	JX502310-11	–	JX502268	KC874562- 63	–	JX502362
31	Bac Kan,Vietnam	N22:30:329 E105:87:600	VN	JX502309	–	JX502279	–	JX502459	–
29	Lang Son,Vietnam	N21:40:881 E106:23:058	VN	JX502306-08	KC874624- 32	JX502269- 73	KC874564- 69	JX502460- 64	JX50236061,63,64
	Total in pearsoni			16	24	20	23	26	26
	Total in this study			26	42	32	41	48	46

DNA sequencing

In this study, we amplified and sequenced introns from 2 X-chromosomal genes (Usp9x and Pola1) and 4 hearing genes (Prestin, Tmc1, FoxP2 and Kcnq4) in bats sampled for an earlier study (sample information and primer details are summarized in Table 2). Sequence data from 2 additional genes, Chd1 and Sws1, were taken from our previous study (Mao et al. 2010).

Table 2.

Primers information for nuclear markers used in this study

Name of markers	ID	Length (bp)	Primers (5′->3′)	References
The nucleosome remodeling factor gene	Chd1	556	F: GATAARTCAGARACAGACCTTAGA CG R: TTTGGCATTCACCTGYACTCC	Lim et al. (2008)
The short-wavelength-sensitive opsin gene	Sws1	645	F: CACAGGCTATGGTGCTGACTT R: GCCCGTGGGGATGGCTATTGA	Mao et al. (2010)
Prestin intron 4	Prestin	536	F: GAGGAGTAAATGCGACCAA R: ATCCCACTGTACCGCTTTG	Mao et al. (2014)
Transmembrane cochlear-expressed gene 1	Tmc1	515	F: AGACAACAAATTCAATTCTATCACA R: GTTAGCGAGAAACCTCAGGAATC	This study
FoxP2 intron 3	FoxP2	530	F: GCTTACCTCAAACCCCTACCA R: CCTGAAGTAAGCAAATGTCCG	Mao et al. (2014)
The voltage-gated potassium channel subfamily KQT member 4	Kcnq4	646	F: GCGTGGTCAAGGTGGAGA R: GCAGGCAGCGTGAATAGAA	Mao et al. (2014)
Ubiquitin specific protease 9 X	Usp9x	674	F:GGCAGACAGGTTGATGACTTGGA R: AGGTCTGCAACTTGCCAAAGGAA	Lim et al. (2008)
Polymerase (DNA directed) alpha 1	Pola1	549	F:GAAACTGGTAGAGCGGAGAA R: ACCTCCCTTCCTTTGTATG	Mao et al. (2014)

Polymerase chain reactions (PCR) were performed in 50 µl reaction mixtures containing 10–50 ng DNA, 0.25 mM of each primer, and 25 µl Premix Taq polymerase (TaKaRa). The thermal profiles for Usp9x, Pola1, FoxP2, Kcnq4, and Prestin have been described previously (Lim et al. 2008; Mao et al. 2014). For Tmc1, we used: 95 °C for 5 min; 34 cycles of 30 s at 94 °C, 40 s at 61 °C, 90 s at 72 °C; 72 °C for 10 min. PCRs were carried out on a PTC-220 thermal cycler (Bio-Rad). DNA sequencing was undertaken with either the forward primer for Tmc1, Usp9x, Pola1, and Kcnq4 or both forward and reverse primers for FoxP2 and Prestin. PCR products were analyzed on an ABI PRISM 3700 automated sequencer (Applied Biosystems). When multiple heterozygous sites were present in the sequences, haplotypes were resolved probabilistically using PHASE 2.1 (Stephens et al. 2001) in the package DnaSP v5 (Librado and Rozas 2009). Sequences were aligned using CLUSTAL_X 1.83 (Thompson et al. 1997) and edited manually. All sequences generated in this study have been deposited in GenBank (see detailed accessions in Table 1).

Gene networks

If differential introgression occurs, phylogenetic relationships between taxa may show differences depending on the type of marker used. At the intraspecific level, gene genealogies are often multifurcated and traditional tree-based phylogenetic methods may be difficult to represent true genealogies (Posada and Crandall 2001). We, therefore, performed network-based phylogenetic reconstructions for each nuclear marker by constructing statistical parsimony networks in the package TCS version 1.21 (Clement et al. 2000).

Gene flow

Shared or closely related haplotypes between R. p. pearsoni (excluding Sichuan) and R. p. chinensis were observed at several nuclear genes (see section ‘Results’), which could have resulted from either introgression or incomplete lineage sorting. To distinguish these 2 processes we ran IM models in the program IMa2 (Hey and Nielsen 2007; Hey 2010). We repeated the IM analysis for each of the 8 loci (Chd1, Sws1, Usp9x, Pola1, Prestin, FoxP2, Kcnq4, and Tmc1) individually. Data for Chd1 and Sws1 were taken from our previous study (Mao et al. 2010). Before performing the IM analysis, for each locus we used DnaSP to test for recombination using the 4-gamete test (Hudson and Kaplan 1985). For loci showing recombination, only those segments without recombination were used in the IM analysis. It was worth pointing out that nonrecombined regions of each marker still showed informative variation between R. p. chinensis and R. p. pearsoni (data not shown). DnaSP was also used to assess neutrality based on the Hudson–Kreitman–Aguade test (HKA, Hudson et al. 1987) and Tajima’s D test (Tajima 1989) whose values were not significant (see Supplementary Table S1). For this reason, and because recent simulations (Strasburg and Rieseberg 2010) have highlighted the robustness of IM models to selection, all of the focal genes were used in the IM models (also see Bull et al. 2006; Pardo-Dieaz et al. 2012). Inheritance scalars were set at 0.75 for 2 X-linked markers (Usp9x and Pola1) and 1 for autosomal markers. For all loci, the Hasegawa-Kishino-Yano (HKY) model was applied. Several preliminary runs were performed to establish upper bounds on prior distributions. To check for the convergence of the Markov chain, the IM analysis was run at least twice using different random seeds. Each run included 200 000 genealogies at every 100 steps after a burn-in of 106 steps including 20 Metropolis-coupled chains with a geometric heating scheme: -hfg -hn20 -ha0.96 -hb0.9. A total of 200 000 genealogies were used to perform likelihood ratio tests of the nested models for migration rates (Hey 2010).

Results

Haplotypes from the 2 X-linked genes (Usp9x and Pola1) were resolved into 3 subnetworks, corresponding to R. p. chinensis, R. p. pearsoni, and a divergent group of R. p. pearsoni from Sichuan (Figure 2A,B). However, 3 of the 4 hearing genes displayed contrasting results to this, with at least 1 haplotype of Prestin, FoxP2, and Tmc1 shared between R. p. pearsoni and R. p. chinensis (Figure 2C–E). It was notable that the shared FoxP2 haplotype between R. p. pearsoni and R. p. chinensis was from populations of their contact zone, Hunan and Fujian. For the fourth hearing gene, Kcnq4, we found a 63-bp deletion in R. p. chinensis compared to R. p. pearsoni (Figure 2F), indicating strong divergence between these 2 taxa at this locus. Like other nuclear genes, networks based on these 4 hearing genes showed that R. p. pearsoni haplotypes from Sichuan were strongly divergent from those from elsewhere. Consequently, individuals of R. p. pearsoni from Sichuan were excluded from estimates of migration rate in the IM analysis.

Figure 2.

Statistical parsimony networks for each nuclear marker used in this study. Haplotypes representing lineages of R. p. chinensis and R. p. pearsoni are shaded orange and blue, respectively. Each circle represents a single haplotype and the area of circle size is scaled by haplotype frequency. The filled black circles represent missing or unsampled haplotypes. Haplytopes were coded as population identities (AH, JX, FJ, SC, HN, GX, GZ, YN, VN) as shown in Figure 1. The arrow in Kcnq4 network denotes a 63-bp deletion (1 mutational step) between R. p. chinensis and R. p. pearsoni.

Two independent IM analysis gave similar posterior probability with the effective sample sizes of > 200 for the migration rate parameter, indicating convergence on the true stationary distribution. To test whether introgression contributed to the observation of shared or closely related haplotypes between R. p. pearsoni excluding Sichuan and R. p. chinensis at several nuclear genes, we compared the fit of models with and without gene flow for all 8 loci individually. Based on likelihood ratio tests, the model with zero gene flow was rejected at 3 of 4 hearing genes (Prestin, FoxP2, and Tmc1 but not rejected at the rest of five loci (see details in Table 3).

Table 3.

Tests of nested models for migration rates between R. p. chinensis to R. p. pearsoni based on the full dataset

Genes	Model	df	2LLR*	P
Chd1	mcp = mpc	1	0.466	0.495
	mpc = 0	1	0.544	0.461
	mcp = 0	1	0.001	1.000
	mcp = mpc = 0	2	2.366	0.306
Sws1	mcp = mpc	1	0.747	0.388
	mpc = 0	1	0.760	0.389
	mcp = 0	1	0.001	1.000
	mcp = mpc = 0	2	0.760	0.689
Prestin	mcp = mpc	1	3.585	0.058
	mpc = 0	1	0.001	1.000
	mcp = 0	1	7.462	0.006
	mcp = mpc = 0	2	12.00	0.003
Tmc1	mcp = mpc	1	3.366	0.067
	mpc = 0	1	0.001	1.000
	mcp = 0	1	2.685	0.101
	mcp = mpc = 0	2	11.56	0.003
FoxP2	mcp = mpc	1	2.80	0.094
	mpc = 0	1	0.001	1.000
	mcp = 0	1	4.197	0.040
	mcp = mpc = 0	2	7.264	0.026
Kcnq4	mcp = mpc	1	0.001	1.000
	mpc = 0	1	0.001	1.000
	mcp = 0	1	0.001	1.000
	mcp = mpc = 0	2	0.001	1.000
Usp9x	mcp = mpc	1	1.762	0.184
	mpc = 0	1	0.001	1.000
	mcp = 0	1	1.708	0.191
	mcp = mpc = 0	2	2.503	0.286
Pola1	mcp = mpc	1	0.523	0.470
	mpc = 0	1	0.001	1.000
	mcp = 0	1	0.401	0.527
	mcp = mpc = 0	2	0.472	0.790

mcp means migration rates from R. p. chinensis to R. p. pearsoni; mpc means migration rates from R. p. pearsoni to R. p. chinensis. * Model improvement was assessed using the log likelihood ratio (LLR) statistics calculated in IMa2, with P values estimated from a chi-squared distribution of 2LLR with the degree of freedom (df). Significant values were shown in bold.

Discussion

Patterns of differential introgression have been frequently used to search for putative speciation genes involved in reproductive isolation and/or beneficial genes which can spread across the species boundaries (see Payseur 2010). In this study, results from 2 X-linked markers (Pola1 and Usp9x) suggested no introgression between R. p. pearsoni and R. p. chinensis supporting previous findings from the other 2 nuclear genes Chd1 and Sws1 (Mao et al. 2010). Furthermore, IMa2 analysis based on likelihood ratio tests could not reject the model of zero gene flow at these 4 genes individually between these 2 subspecies, perhaps indicating these genes are involved in reproductive isolation either directly or via linkage to other genes.

In contrast to the above patterns, 3 of 4 hearing genes (Prestin, Tmc1, and FoxP2) exhibited shared and/or closely related haplotypes between R. p. pearsoni and R. p. chinensis. While this result could in theory be explained by either incomplete lineage sorting or introgression (Funk and Omland 2003; Ballard and Whitlock 2004), the results of the IMa2 analyses supported the latter scenario, with the rejection of the model of zero gene flow at these 3 hearing genes when analyzed individually based on likelihood ratio tests. This result was consistent with our previous finding that Prestin appeared to show gene flow across the hybrid zone between 2 subspecies of the congeneric species R. affinis (Mao et al. 2014). More horseshoe bat taxa need to be studied to test the generality of this pattern.

Several scenarios can be considered to explain the pattern of increased rates of introgression observed in 3 of 4 hearing genes examined. First, it is possible that these 3 hearing genes in fact provide an adaptive advantage in a heterospecific background (Arnold 2006; Pardo-Dieaz et al. 2012; Hedrick 2013). Indeed in mice (Mus), genes that function in olfaction were shown to be subject to adaptive introgression across a hybrid zone (Teeter et al. 2008). Our neutrality tests failed to support evidence of selection acting on genes examined here; nonetheless, it is known that strong adaptation can occur in the absence of detectable signatures of selection (e.g., Mc1r gene in mice, Domingues et al. 2012) and therefore we cannot rule this out completely. If introgression of these hearing genes was beneficial, these genes might not be involved in echolocation call frequency. Otherwise, hybrids would be particularly selected against due to quite different call frequency between their parental taxa (R. p. pearsoni and R. p. chinensis, see Mao et al. 2010). Alternatively, these hearing genes examined may be linked to loci that can cross the species boundaries due to positive selection. Ultimately, functional analysis on additional candidate hearing coding gene sequences from individuals of the 2 focal taxa would be needed to test more thoroughly for adaptive introgression associated with what is likely to be a complex phenotypic trait.

Third, the observed transfer of alleles across taxon boundaries may have arisen via stochastic processes (i.e., genetic drift), and it is often difficult to distinguish the roles of these 2 processes in introgression events (but see Payseur et al. 2004; Teeter et al. 2008; Fitzpatrick et al. 2009). This is especially likely to be the case if these hearing genes under study do not directly impact on echolocation call frequency per se, but rather function in other aspects of this complex trait.

Although not a focus of our study, the observed strong levels of both mitochondrial and nuclear differentiation between R. p. pearsoni individuals from Sichuan versus those from adjacent populations strongly point to the presence of a cryptic taxon. In addition, published differences in diploid chromosome number and chromosomal rearrangements between R. p. pearsoni from Sichuan (2N = 46, Wu et al. 2009) and ones from other regions (e.g., 2N = 44 in Guizhou, Mao et al. 2007) also support either different taxa or distinct chromosomal races. Such chromosomal rearrangements are well known to reduce gene flow and thus increase genetic differentiation, for example, by suppressing recombination (Ortiz-Barrientos et al. 2002; Navarro and Barton 2003).

In conclusion, parapatric taxa that undergo genetic exchange offer good opportunities to identify candidate loci that cross taxonomic barriers versus those that resist gene flow and thus might be related to reproductive isolation. By examining patterns of differential introgression among candidate loci, we revealed evidence of increased introgression from R. p. chinensis to R. p. pearsoni at 3 of 4 hearing genes and reduced introgression at 2 X-linked and 2 autosomal loci. However, we were unable to explicitly relate gene flow across species barriers to phenotypic differences in the relevant individuals. Although this study is one of the first to test for introgression of sensory genes among different taxa, our statistical power to find effects was limited by our low coverage of the genome. To address this issue, as well as known heterogeneity in genomic divergence (reviewed in Nosil et al. 2009), high-throughput sequencing approaches (e.g., whole-genome resequencing) offer promise for more thoroughly assessing genetic differentiation and introgression in these and other taxa (Twyford and Ennos 2012; Martin et al. 2013).

Supplementary Material

Supplementary material can be found at http://www.cz.oxfordjournals.org/

Click here for additional data file.