A maladaptive combination of traits contributes to the maintenance of a stable hybrid zone between two divergent species of Drosophila

Geographical areas where two species come into contact and hybridize serve as natural laboratories for assessing mechanisms that limit gene flow between species. The ranges of about half of all closely related Drosophila species overlap, and the genomes of several pairs reveal signatures of past introgression. However, only two contemporary hybrid zones have been characterized in the genus, and both are recently diverged sister species (D. simulans-D. sechellia, Ks = 0.05; D. yakuba-D. santomea, Ks = 0.048). Here we present evidence of a new hybrid zone, and the ecological mechanisms that maintain it, between two highly divergent Drosophila species (Ks = 0.11). On the island of Bioko in west Africa, D. teissieri occupies mostly forests, D. yakuba occupies mostly open agricultural areas, and recently, we discovered that hybrids between these species occur near the interface of these habitats. Genome sequencing revealed that all field-sampled hybrids are F1 progeny of D. yakuba females and D. teissieri males. We found no evidence for either advanced-generation hybrids or F1 hybrids produced by D. teissieri females and D.yakuba males. The lack of advanced-generation hybrids on Bioko is consistent with mark-recapture and laboratory experiments that we conducted, which indicate hybrids have a maladaptive combination of traits. Like D. yakuba, hybrids behaviorally prefer open habitat that is relatively warm and dry, but like D. teissieri, hybrids have low desiccation tolerance, which we predict leaves them physiologically ill-equipped to cope with their preferred habitat. These observations are consistent with recent findings of limited introgression in the D. yakuba clade and identify an ecological mechanism for limiting gene flow between D. yakuba and D. teissieri; namely, selection against hybrids that we have documented, in combination with hybrid male sterility, contributes to the maintenance of this narrow (~30m), stable hybrid zone centered on the forest-open habitat ecotone. Our results show how a deleterious combination of parental traits can result in unfit or maladapted hybrids.

Genetically distinct populations of closely related species often come into contact, mate, 27 and produce offspring in narrow zones of hybridization (Barton and Hewitt 1985;Hewitt 1988). 28 These geographical areas offer windows on the evolutionary processes that lead to the origin and 29 maintenance of species (Harrison 1990), and have been of interest to evolutionary biologists for 30 over a century (Harrison and Larson 2016). Incipient species could potentially fuse into a single 31 entity if fertile hybrids are produced and successfully backcross into one or both parental species, 32 but because selection has not tested the combinations of alleles found in divergent species, 33 hybrids are often less fit than either parental species (Dobzhansky and Dobzhansky 1937;Muller 34 1942;Turelli et al. 2001). For example, hybrids have been shown to have both developmental 35 (Mendelson 2003;Coyne and Orr 2004;Moyle et al. 2004;Maheshwari and Barbash 2011) and 36 behavioral defects (Coyne and Orr 1989;Funk et al. 2006;Turissini et al. 2017). These 37 reproductive isolating mechanisms should contribute to the maintenance of species in areas of 38 contact if they prevent the formation of hybrids, or if they lead to the production of sterile and/or 39 in Drosophila (Turissini and Matute in revision). In the most extreme cases, successful hybridization seems possible between species pairs of lancelets (120 my diverged) (Holland et al. 2015), ferns (60 my diverged) (Rothfels et al. 2015), and fungi (at least 3% nuclear divergence) 48 (Stukenbrock et al. 2012). Thus, while estimated divergence times remain error prone and 49 inexact, these data indicate that hybridization between highly diverged species is at least 50 possible. 51 About half of all closely related Drosophila species have overlapping geographical 52 ranges (Turelli et al. 2014), and several species pairs show evidence of past introgression 53 (Shoemaker et al. 1999;Jaenike et al. 2006;Kulathinal et al. 2009;Garrigan et al. 2012;Brand et 54 al. 2013;Lohse et al. 2015), but only two contemporary hybrid zones have been identified and 55 well described in the genus. Both hybrid zones involve closely related sister species in the D. in west Africa Llopart et al. 2005a). Along the altitudinal transect of Pico 61 de São Tomé, D. yakuba occurs at low elevations (below 1,450 m), D. santomea occurs at high 62 additional set of 100 males of each genotype and calculated Mahalanobis distances for each 115 individual (Mahalanobis 1936). The Mahalanobis distance between a focal individual (i) and the 116 average for a given genotype (i.e., centroids for D. teissieri, D. yakuba, and the two hybrid F 1 117 genotypes estimated using the training set) was calculated as: 118 where the super-index T denotes matrix transpose, S denotes the covariance matrix of a given 122 dataset, F i is the vector of phenotypic observations in a focal individual, i, and μ LB is the vector 123 of average phenotypic observations of the training set. To estimate the accuracy of this approach, 124 we calculated the proportion of our blind assignments that were correct. Mahalanobis distances 125 were later used in the field to distinguish genotypes, and the field assignments were later verified 126 using genomic sequencing (see below). 127 each sampled altitude. On day 15, the traps at 2,020 m became inaccessible for safety reasons. extracted from each individual using the Beckman-Coulter DNAdvance magnetic bead protocol using the Tagmentase protocol detailed in Picelli et al. (2014), and barcoded for multiplexing 163 (see supplement for barcode sequences). Libraries were then sequenced to low coverage with 164 Illumina 100 bp single end reads (Cornell Genomics Facility). In total, 20 parental D. teissieri, 165 17 parental D. yakuba, and 19 individuals of putative hybrid phenotype were sequenced to 166 sufficient coverage (i.e., > 10,000 reads). The mean number of markers per individual was 6.869, 167 and the mean coverage was 2.29X. 168 Genotypes of the individuals were determined by the Multiplexed Shotgun Genotyping 169 (MSG) pipeline described by Andolfatto et al. (2011), which uses a hidden Markov model 170 (HMM) to assign ancestry along a genome with low-coverage read data. Because this approach 171 utilizes linkage disequilibrium on a large physical scale, it is well-suited for assigning genome-172 wide ancestry in recently formed hybrids (within several generations of backcrossing). The D. 173 yakuba Flybase assembly (version 1.05) was repeat-masked (Smit et al. 2013(Smit et al. -2015 and used as 174 the first parental reference. The second parental input for MSG was made by mapping reads from 175 the outbred individual Bioko_cascade_2_2 (Turissini and Matute in revision) to the D. 176 yakuba flybase assembly and creating an updated FASTA file with genotype calls from GATK 177 (McKenna et al. 2010). The updated FASTA file uses only single nucleotide polymorphisms 178 (SNPs), and masks all inferred indels plus 5 bases both up and downstream. Some regions of the 179 genome are error-prone in terms of assigning ancestry, due in part to low sequence divergence, 180 poor reference genome assembly, or high polymorphism in the parental species. To reduce the 181 rate of miscalled ancestry, known intermediate-frequency SNPs in each parental population were masked in the corresponding reference genomes using sequence data from several wild-

Distributions of D. yakuba, D. teissieri, and their hybrids across the ecotone 240
The male individuals that we recaptured during our mark-recapture experiment enabled 241 us to estimate the rate of hybridization and to assess the distributions of each genotype across the 242 ecotone. Flies collected during this experiment included both dusted flies and any other D. 243 yakuba clade flies visiting our traps at the time of recapture. We combined these data with our 244 original capture data to estimate the prevalence of D. yakuba genotypes across the ecotone. We 245 first calculated the dilution factor, for each genotype, defined as the number of non-dusted flies 246 captured divided by the number of dusted flies recaptured (White et al. 1982). We then 247 multiplied the number of flies that we originally released by the dilution factor to estimate the 248 census size of each genotype at this single altitudinal sliver of habitat. The estimated census size 249 of hybrids divided by the sum of the three estimated census sizes (D. yakuba, D. teissieri, and 250 hybrids) provides an estimate for hybrid prevalence. To assess variation in the prevalence of D.
open habitat at five of the seven altitudes where flies were sampled (approximately 200, 650, 257 1,200, 1,650, and 2,020 m). Temperature was measured at 4 times of the day using a Digi Sense 258 field thermometer (Catalog number: 86460-05; Cole-Parmer Instrument Company, Vernon Hills, 259 IL). All measurements were made 10 cm away from the floor. In parallel, and at the same time, 260 we measured air humidity using a field hygrometer (Dwyer MST2-01 Moisture Meter; Michigan 261 City, IN). We recorded these data at five locations per altitude, for each of the two environments, 262 at five altitudes, resulting in a total of 50 measured sites. At each site, we recorded temperature 263 and humidity at 0600, 1200, and 1700 hrs. Each site was measured on three different days for a 264 total of 450 observations. 265 To assess variation in temperature and humidity, we fitted two linear models, one for 266 temperature and one for humidity. The two models followed the same form and assessed the 267 effect of altitude, type of habitat, and time of the day on each of the two environmental factors. 268 We also included all possible interactions. We compared temperature and humidity at open 269 versus forest habitats using the Tukey's Honest Significant Difference (HSD) post hoc pairwise 270 comparisons in the 'multcomp' library (Hothorn et al. 2008). 271 ♂ yak) 274 and F 1 (♀yak × ♂ tei) hybrids in the laboratory. Flies were randomly placed in a thermocline that 275 consisted of a plexiglass chamber (12 cm wide × 45 cm long × 1 cm high) with an aluminum 276 floor. The thermocline was placed in an 18ºC room and heat plates (120 VAC Thermo Scientific 277 Cimarec Hot Plate, Thermo Scientific Cimarec, # UX-04600-01, Waltham, MA, USA) were 278 used to generate a thermal gradient ranging from 18ºC to 30ºC, with a change in temperature of 279 approximately 2ºC every 6 cm. This range of temperatures encompasses the majority of the 280 temperatures flies might regularly experience on Bioko (see below). Flies were allowed to move 281 freely along this gradient over a period of one hour. At the end of each trial we isolated flies into 282 seven chambers, each 10.5 × 6 × 1 cm, by pushing a rod connected to six plexiglass partitions 283 across the width of the chamber. We then recorded the temperature within each partition using a 284 Digi-Sense thermometer equipped with a type-T thermocouple (Cole-Parmer Instrument Co., 285 Chicago, IL; catalog number: 86460-05) and counted the number of flies within each partition. 286 Males and females of each genotype were evaluated separately to avoid sexual attraction that 287 might influence results (N = 3 separate replicates per sex). 288 Because temperature preference was not normally distributed, and remained non-289 Gaussian after both log (Shapiro Wilk test, W = 0.914, P < 0.0001) and square root (Shapiro 290 Wilk test, W = 0.919, P < 0.0001) transformations, we used the 'ARTool' library to complete 291 analyses of variance on aligned rank transformed data (Kay and Wobbrock 2016). We 292 specifically assessed the effects of genotype, sex, and their interactions on temperature 293 preference. We then used lsmeans to conduct post hoc pairwise comparisons (Lenth 2016a). 294 Relative humidity (RH) preferences of D. teissieri, D. yakuba, and F 1 (♀tei × ♂ yak) and 297 F 1 (♀yak × ♂ tei) hybrids were evaluated by giving flies the choice of orienting themselves along saturated salt solutions: LiCl, NaCl, or KH 2 PO 4 . Each of these solutions generates a RH of 301 ~20%, ~70%, and ~85%, respectively, in the headspace above the rows. We filled the wells of 302 two adjacent end rows with LiCl, the three following rows with NaCl, and the three remaining 303 rows with KH 2 PO 4. This generates a gradient of RH ranging from ~20% to 85%. The top of each 304 plate was covered with 300 micron nylon netting (MegaView Science Co., Ltd. Taichung, 305 Taiwan) and covered with the culture plate lid on top of the mesh; this left ~1 cm for flies to 306 move freely around the plate. This experimental design is a modified approach from Enjin and 307 colleagues (2016). We lightly anesthetized approximately 50, 4 to 7-day old virgin, males or 308 females of a given genotype and placed them along the long axis of a plate. We ran eight plates 309 simultaneously: one plate for each of the parental and hybrid genotypes, with sexes evaluated 310 separately. Flies were allowed to orient themselves for the first hour after which pictures were 311 taken every 15 min, for an additional two hours. This procedure was repeated on four separate 312 days with the position of each plate and the orientation of the 'low' and 'high' RH end of the 313 plates relative to the room randomized each day. (This avoids confounding effects of non-314 uniform lighting and other conditions among days in the laboratory.) In total, we assayed ~200 315 individuals of each genotype and sex. To score preference, we counted the number of flies over 316 each well of the plates, for each of the eight images generated over the 2 hr assay, and summed 317 counts of flies oriented over wells containing the same super-saturated salt solution. All scoring 318 was double blinded.
Because there was no effect of sex on the number of flies choosing a given RH (LRT: χ 2 320 = 0.25; P = 0.62), sexes were analyzed together. We first tested whether humidity preference 321 varied across genotypes by modeling the mean number of individuals choosing a given RH 322 (Poisson) as a function of genotype, RH, and the interaction between genotype and RH. A 323 dummy 'plate identity' variable was included as a random effect to account for among-plate 324 variation such as the location of the plate in the room, slight differences in solution volume 325 within individual wells, or variation in temperature among days. We tested whether each fixed 326 effect influenced the number of flies choosing a given RH using LRTs that compared models 327 that included and excluded each term. We dropped each fixed effect independently, starting with 328 the interaction term. To explicitly test whether different genotypes displayed preferences for 329 different RH, we also modeled the mean number of individuals choosing a given humidity, for 330 each genotype separately, as a function of RH (fixed effect) and plate identity (random effect). 331 We tested for variation in humidity preference using LRTs comparing these genotype-specific 332 models to those lacking the fixed effect of RH. For genotypes that had significant variation in 333 humidity preference we carried out pairwise contrasts between each RH using the lsmeans and 334 pairs functions of the 'lsmeans' library (Lenth 2016b). 335 336 Desiccation tolerance in the laboratory 337 Desiccation resistance was measured by placing ten, 4-day old virgin, females or males in 338 30 mL empty vials (N = 11 vials per sex), which in turn were placed in a glass desiccator with 339 generated with the 'survplot' library. To assess any effects of body size on desiccation tolerance, 346 we measured the thorax length of 4-day old males and females (N = 50 of each sex from each 347 genotype) reared at 24ºC. We then fitted a linear model to evaluate the fixed effects of sex and 348 genotype, and the interaction between these factors, using the function lm function in the 'stats' 349 package. 350

Morphological characteristics differ among D. yakuba, D. teissieri, and F 1 hybrids 353
Before sampling the island of Bioko, we first confirmed whether D. yakuba, D. teissieri, 354 and F 1 hybrids could be reliably distinguished based on morphology in the laboratory. We found 355 that while D. teissieri males always have at least three strong chitinized anal spines ( significantly in Tukey's HSD post hoc pairwise comparisons (P < 0.037 for all statistically 363 significant comparisons), with the exception of reciprocal hybrids (P = 0.203) and the comparison between F 1 (♀yak × ♂ tei) hybrids and D. teissieri (P = 0.888). In contrast, while the 365 number of teeth in sex combs differed among all groups (F 3,1996 = 220.11, P < 0.0001), F 1 (♀yak 366 × ♂ tei) hybrids had fewer teeth in their sex combs (4.810 ± 1.328 SD) than did D. teissieri 367 (5.238 ± 0.826 SD), D. yakuba (6.224 ± 0.955 SD), and F 1 (♀tei × ♂ yak) hybrids (6.114 ± 0.952 368 SD). The latter two genotypes were statistically indistinguishable (P = 0.332), but all other post 369 hoc comparisons were significant (P < 0.0003 for all statistically significant comparisons). Thus, 370 sex chromosomes and/or the cytoplasmic factors likely influence the number of teeth in sex 371 combs. 372 We next assessed our ability to reliably identify pure species and hybrids using 373 Mahalanobis distances. Figure 1 shows the distributions for the number of anal spines ( Figure  374 1A), the number of teeth in sex combs ( Figure 1B), and the lengths of tibia ( Figure 1C) for the 375 four male genotypes. We were able to reliably identify pure species (D. yakuba: 100/100, D. 376 teissieri: 100/100) and F 1 hybrids (single class, 196/200). We attempted to gain further resolution 377 and assessed our ability to discriminate between F 1 (♀tei × ♂ yak) and F 1 (♀yak × ♂ tei) reciprocal 378 hybrids, but due to their similar morphology our assignments were incorrect 53% of the time. 379 Thus, for all later field experiments we treated hybrids as a single class, and relied on subsequent 380 genome sequencing to genotype hybrids. 381  To confirm our field genotype assignments based on Mahalanobis distances alone, we 406 sequenced a subset of the male individuals sampled from Bioko, including putative hybrids. All 407 19 sequenced hybrid males are heterozygous for D. yakuba and D. teissieri ancestry across their 408 autosomes and hemizygous for the D. yakuba X (Figure 2). This indicates that these hybrids are 409 the sons of D. yakuba females and D. teissieri males. There are many small regions for which 410 genotypes could not be confidently called (gray blocks in Figure 2). Taken together, these data 411 confirm our discovery of hybrids on Bioko based on Mahalanobis distances in the field, but these 412 data also suggest that F 1 (♀tei × This included 91 D. teissieri, 137 D. yakuba, and 43 hybrid (dusted) males that we released. We 462 estimate that with our measured recapture rates (11.31% for D. yakuba and 8.45% for D.
clade individuals in this area. These are clearly rough estimates, but they suggest that hybrids are 475 not common. 98% of the hybrids that we recaptured traveled less than 30 m from the ecotone, 476 suggesting that this hybrid zone is very narrow (Figure 2). Drosophila teissieri was mostly 477 absent from the open habitat and its prevalence increased with distance into the forest habitat  Table 1. Statistical analyses are presented in Table 2. 509

Hybrids have a maladaptive combination of parental phenotypes 512
The F 1 (♀yak × ♂ tei) hybrids we sampled on Bioko behaviorally prefer warm and dry 513 open areas 68% of the time in the field-these areas are also chosen by D. yakuba 92% of the 514 time. In contrast, D. teissieri chooses these areas only 3% of the time. We next sought to 515 determine whether these preferences exist in the laboratory, and if so, whether the physiological 516 tolerances of D. yakuba, D. teissieri, and hybrids to low humidity conditions corresponds to their 517 behavioral preferences in the field. 518 Drosophila yakuba and F 1 hybrids prefer warmer temperatures than D. teissieri. We 519 assessed the temperature preference of D. yakuba, D. teissieri, and their F 1 hybrids in the 520 laboratory. Temperature preferences for the eight lab-produced genotypes [(2 pure species + 2 521 reciprocal hybrids) × 2 sexes] are shown in Figure 5. The mean temperature preference of D. 522 teissieri (21.285ºC ± 2.905 SD, N = 1,200) was approximately 13% lower than D. yakuba 523 (24.51ºC ± 3.68 SD, N = 1,200), 12% lower than F 1 (♀tei × ♂ yak) (24.12ºC ± 3.72 SD, N = 524 1,200), and 10% lower than F 1 (♀yak × ♂ tei) hybrids (23.74ºC ± 3.82 SD, N = 1,200) (F 3,4792 = 525 209.927, P < 0.0001). Tukey's HSD post hoc comparisons revealed the temperature preference 526 of F 1 (♀tei × ♂ yak) hybrids differed from D. teissieri (P < 0.0001), but not D. yakuba (P = 0.05). 527 In contrast, F 1 (♀yak × ♂ tei) hybrids differed from both D. teissieri (P < 0.0001) and D. yakuba 528 (P < 0.0001). These patterns stem from a slight, but statistically significant, difference between 529 the two reciprocal hybrids (P = 0.03). Across all genotypes females preferred higher 530 temperatures (23.87ºC ± 3.70 SD) than did males (22.95ºC ± 3.78 SD) (F 1,4792 = 57.361, P < 531 0.0001). Finally, we found a modest, but statistically significant, interaction between genotype 532 and sex (F 3,4792 = 3.207, P < 0.022). These results indicate that F 1 hybrid genotypes have a preference for warm temperatures that is similar to D. yakuba in the laboratory, in addition to the field, suggesting that D. yakuba's preference for high temperatures is dominant or semi-dominant 535 to D. teissieri's preference for cool temperatures. χ 2 = 24.4; P < 0.0001; 553 Figure 6], but F 1 (♀yak × ♂ tei) hybrids show a weaker preference than F 1 (♀tei × ♂ yak) hybrids. 554 This is indicated by the only marginally significant pairwise contrasts of the number of F 1 (♀yak 555 × ♂ tei) individuals at different humidity (P = 0.06 for both 20% -70% RH and 20% -85% RH), 556 but highly significant contrasts for F 1 (♀tei × ♂ yak) hybrids (P = 0.0004 for both 20% -70% RH 557 and 20% -85% RH). These data indicate that F 1 hybrid genotypes have a preference for low 558 humidity that is similar to D. yakuba, suggesting that D. yakuba's preference for low humidity is 559 dominant or semi-dominant to D. teissieri's lack of preference. 560 561 FIGURE 6. Drosophila teissieri prefers higher relative humidity (%) than does D. yakuba, 562 we assessed whether D. yakuba, D. teissieri, and their F 1 hybrids differ in physiological tolerance 571 to osmotic stress. The body size of all genotypes was similar suggesting that any differences in 572 desiccation tolerance do not depend on size (F 3, 392 = 2.4335, P = 0.065). As expected, males 573 were more prone to suffer from desiccation than were females (~between 1% and 10% less 574 tolerant depending on the genotype; Cox hazards regression, sex effect: χ 2 = 8.39, df = 1, P = 575 0.004) (Table S2)-the smaller relative body size of males likely influences this difference (F 1, 576 392 = 295.25, P < 0.0001). The largest differences in desiccation tolerance were observed 577 between genotypes with D. yakuba having the highest desiccation tolerance (Cox hazards 578 regression, genotype effect: χ 2 = 143.55, df = 3, P < 0.0001) (Figure 7). The tolerance of D. 579 teissieri to osmotic stress (4.73 hr ± 1.69 SD) is 40% lower than the tolerance of D. yakuba (6.63 580 hr ± 1.54 SD), consistent with the distribution of these species in nature. Together, hybrids from 581 both reciprocal crosses have approximately 28% lower desiccation tolerance than D. yakuba, 582 despite behaviorally preferring such conditions. F 1 (♀yak × ♂ tei) hybrids that occur naturally on 583 Bioko had the lowest desiccation tolerance of all genotypes (4.50 hr ± 1.47 SD), however, these 584 hybrids did not differ statistically from F 1 (♀tei × ♂ yak) hybrids (4.76 hr ± 1.50 SD). This 585 suggests little influence of the cytoplasm or sex chromosomes on desiccation tolerance. Finally, 586 F 1 (♀yak × ♂ tei) hybrids have significantly lower desiccation tolerance than does D. teissieri, but 587 F 1 (♀tei × ♂ yak) hybrids do not differ statistically from D. teissieri. Means, standard deviations, 588 sample sizes, and statistical comparisons are reported in Table 4. Taken together, these results 589 indicate that D. teissieri alleles influencing desiccation tolerance are dominant or semi-dominant 590 to D. yakuba alleles. Thus, F 1 hybrids, and especially F 1 (♀yak × ♂ tei) hybrids, are predicted to FIGURE 7. Both male (A) and female (B) D. yakuba individuals have higher desiccation 593 tolerance than do D. teissieri, F 1 (♀tei × ♂ yak), and F 1 (♀yak × ♂ tei) genotypes. Across sexes 594 F 1 (♀yak × ♂ tei) have the lowest desiccation tolerance indicating that these hybrids are ill-595 equipped to cope with the dry environments they prefer in the laboratory and in nature. Summary 596 statistics are reported in Table S2 for males and females. Table 4 reports summary statistics and 597 the results of Tukey's HSD post hoc pairwise comparisons between genotypes. 598 conditions (Turissini et al. 2015), and notably, both F 1 (♀yak × ♂ tei) and F 1 (♀tei × There are other fitness costs that likely contribute to the maintenance of this narrow 627 hybrid zone. First, all hybrid males are sterile (Turissini et al. 2015). Second, teissieri hybrids have behavioral defects that compromise their foraging behavior in ways 629 predicted to reduce hybrid fitness (Turissini et al. 2017) revealed any obvious differences in the fertilization success or sterility/viability between F 1 643 hybrids produced by reciprocal crosses (Turissini et al. 2015;Cooper et al. 2017), as observed in 644 other systems (i.e., "Darwin's corollary", Turelli and Moyle 2007). Llopart et al. 2005b;Bachtrog et al. 2006). This discrepancy seems to be due to the number of 651 loci included in different analyses and on the methods used to assess amounts of introgression. 652 For example, the most recent published whole genome analysis of D. yakuba and D. santomea 653 mtDNA suggested pervasive introgression (Llopart et al. 2014), but this study relied on IMa2 654 Nielsen 2004, 2007), which is known to provide false positives for migration 655 (Cruickshank and Hahn 2014;Hey et al. 2015). 656 Only two other hybrid zones have been identified in the Drosophila, and both of these 657 examples are for recently diverged sister species in the D. melanogaster subgroup. The sister 658 species D. simulans and D. sechellia hybridize in the central Seychelles. F 1 hybrid males 659 between these two species have been found in ripe Morinda ripe fruits. Such evidence indicates 660 that the ability to breed on Morinda could be transferred from D. sechellia to D. simulans, and 661 preliminary analyses of interspecific introgression indicate that D. simulans and D. sechellia may 662 currently be exchanging genes (Solignac and Monnerot 1986;Kliman et al. 2000;Garrigan et al. the secondary forest where D. teissieri dwells is also recent, having grown after the downfall of the military dictatorship in Equatorial Guinea . Thus, while the exact dates are 670 uncertain, agricultural expansion suggests that secondary contact in the D. yakuba clade may be 671 recent. 672 The discovery of the D. yakuba-D. teissieri hybrid zone on Bioko (versus other areas of 673 range overlap) may not be surprising, as hybridization is thought to be common in this region of 674 Africa. For example, the reed frogs Hyperolius molleri and H. thomensis hybridize in the 675 midlands of Pico de São Tomé in an area that directly overlaps with the D. yakuba-D. santomea 676 hybrid zone (Bell et al. 2015). Limited data suggest that Caecilians (Schistometopum thomense 677 and S. ephele) may also hybridize on São Tomé (Stoelting et al. 2014 showed that many of the recombinant genotypes present in germinated seed-sets are not recruited 686 into the adult population. Selection acting on recombinant genotypes seems to contribute to RI 687 between P. alba and P. tremula (Christe et al. 2016). More subtle patterns of genetic variation 688 observed within advanced-generation hybrids, such as the overrepresentation of certain parental 689 alleles in linkage disequilibrium can be used to infer how selection acts to maintain species 690 boundaries (Schumer et al. 2014). Our results add to these studies by suggesting a phenotypic mechanism that likely contributes to a lack of advanced-generation hybrids in the Bioko D.