{'fraction_non_alphanumeric': 0.03790650637456537, 'fraction_numerical': 0.024793764061541258, 'mean_word_length': 5.187289088863892, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Genic rather than genome‐wide differences between sexually deceptive Ophrys orchids with different pollinators

High pollinator specificity and the potential for simple genetic changes to affect pollinator attraction make sexually deceptive orchids an ideal system for the study of ecological speciation, in which change of flower odour is likely important. This study surveys reproductive barriers and differences in floral phenotypes in a group of four closely related, coflowering sympatric Ophrys species and uses a genotyping‐by‐sequencing (GBS) approach to obtain information on the proportion of the genome that is differentiated between species. Ophrys species were found to effectively lack postpollination barriers, but are strongly isolated by their different pollinators (floral isolation) and, to a smaller extent, by shifts in flowering time (temporal isolation). Although flower morphology and perhaps labellum coloration may contribute to floral isolation, reproductive barriers may largely be due to differences in flower odour chemistry. GBS revealed shared polymorphism throughout the Ophrys genome, with very little population structure between species. Genome scans for FST outliers identified few markers that are highly differentiated between species and repeatable in several populations. These genome scans also revealed highly differentiated polymorphisms in genes with putative involvement in floral odour production, including a previously identified candidate gene thought to be involved in the biosynthesis of pseudo‐pheromones by the orchid flowers. Taken together, these data suggest that ecological speciation associated with different pollinators in sexually deceptive orchids has a genic rather than a genomic basis, placing these species at an early phase of genomic divergence within the ‘speciation continuum’.


Introduction
Speciation is characterized by the formation of barriers to gene flow between populations that will ultimately result in genetic divergence (Coyne & Orr 2004;Lowry et al. 2008). Whereas speciation research historically focused upon the geographic mode of speciation (e.g. sympatric vs. allopatric), focus has recently shifted on the processes of divergence (Rundle & Nosil 2005;Wolf et al. 2010;Smadja & Butlin 2011;Feder et al. 2012;Nosil 2012). Barriers to gene flow can be mediated by many genes of small effect that are distributed throughout the genome, or by few genes with large effect on reproductive isolation (RI); this latter situation is conceptualized in the genic view of speciation (Wu 2001;Wu & Ting 2004;Lexer & Widmer 2008). In genic speciation, only few loci under divergent selection are responsible for species differences, whereas the rest of the genome may be 'porous' to gene flow (i.e. subject to allele exchange between diverging species) (Wu 2001;Wu & Ting 2004;Lexer & Widmer 2008). This conceptually overlaps with ecological speciation, in which divergent selection by the (biotic or abiotic) environment drives species divergence, a process which may in principle proceed in the presence of gene flow (Nosil 2012;Rundle & Nosil 2005; Smadja & Butlin 2011; but see Cruickshank & Hahn 2014). It has been proposed that the overall accumulation of genomic divergence between incipient species best be seen as a 'speciation continuum', which proceeds from a genic stage with direct selection on few loci, via 'divergence hitchhiking' (Via & West 2008) and 'genome hitchhiking' to completed speciation and postspeciation genome divergence (Feder et al. 2012;Nosil 2012). However, the question of the relative importance of these postulated phases of speciation is debated, and it is therefore important to assess, in as many cases as possible, the phase of genomic divergence that has been reached in any given case of ecological speciation (Feder et al. 2012;Nosil 2012).
Early stages of speciation may allow the identification of the initial barriers that caused species divergence, along with their genomic architecture (Rundle & Nosil 2005). Likewise, incipient speciation may shed light on the order in which different reproductive barriers evolve (Rundle & Nosil 2005;Lexer & Widmer 2008). It has been postulated that the evolution of premating barriers, such as pollinator-mediated RI, often precedes the emergence of postmating barriers (Coyne & Orr 2004;Lowry et al. 2008). Moreover, such early evolving barriers probably involve alleles of large effect, the action of which can later be modified by minor-effect alleles (epistasis) as species divergence continues (Coyne & Orr 2004;Widmer et al. 2009). Putative cases of pollinator-driven speciation (Johnson 2006;Schiestl 2012), in which divergence of plant species is due to strong pollinator-mediated RI (floral isolation; Schiestl & Schl€ uter 2009), may therefore be prime candidates for genic ecological speciation processes, especially where the genetic control of pollinator attraction is simple (Schl€ uter et al. 2011b;. Sexually deceptive plants achieve high pollinator specificity by mimicry of pollinator females (Schiestl 2005) and provide good study systems for ecological speciation and the evolution of RI (Peakall et al. 2010;Schl€ uter et al. 2011a;Xu et al. 2012b;. In the best-studied cases of sexual deception by plants, namely the orchid genera Ophrys in the Mediterranean and Chiloglottis in Australia, speciation is often associated with pollinator shifts (Paulus & Gack 1990;Peakall et al. 2010;Ayasse et al. 2011;Xu et al. 2012b). In particular, Ophrys bears all the hallmarks of ecological speciation (sensu Nosil 2012), namely (i) divergent selection and (ii) RI due to pollinators, and (iii) a genetic mechanism linking these two (Schl€ uter & Schiestl 2008;Xu et al. 2012b). Both Ophrys and Chiloglottis attract pollinators by chemical mimicry of insect sex pheromones (Schiestl 2005;Schl€ uter & Schiestl 2008;Ayasse et al. 2011). In Ophrys, a blend of cuticular hydrocarbons, especially alkanes and alkenes, is responsible for pollinator attraction (Schl€ uter & Schiestl 2008;Ayasse et al. 2011), and simple genetic changes affecting hydrocarbon biosynthesis can have a drastic effect on pollination; simple mutations may therefore lead to pollinator shifts, floral isolation and speciation (Schl€ uter et al. 2011b;Xu et al. 2012a,b).
This study investigates four closely related Ophrys species (O. exaltata, O. garganica, O. incubacea and O. sphegodes), among which the published literature indicates a lack of pollinator sharing that is expected to result in floral isolation (Paulus & Gack 1990;Gaskett 2011). Previous work suggests low levels of genetic differentiation, but odour differentiation and strong floral isolation (Soliva & Widmer 2003;Mant et al. 2005;Xu et al. 2011;Breitkopf et al. 2013), and has revealed candidate genes for pollinator attraction (Schl€ uter et al. 2011b;Xu et al. 2012a;Sedeek et al. 2013). However, as those studies examined different subsets of species, a complete picture of (genetic and phenotypic) species differences and reproductive barriers in this species group is still lacking. Moreover, previous population genetic studies (Soliva & Widmer 2003;Mant et al. 2005;Xu et al. 2011;Breitkopf et al. 2013) were limited by a small number of genetic markers and hence do not provide detailed information on the genomics of species divergence. This study therefore aims to provide a comprehensive picture of species and speciation in Ophrys by addressing the following questions: (i) are the study species reproductively isolated and if so, due to which reproductive barriers; (ii) to what extent are the species phenotypically differentiated and distinct; (iii) what proportion of the genome is associated with species differences (i.e. which stage of divergence in the speciation continuum has been reached?); and (iv) can genes associated with species divergence be identified?

Study species and plant material
This study examines sympatric populations of closely related Ophrys species (Fig. S1 (Paulus & Gack 1990;Xu et al. 2011). Initial species identification in the field was based upon floral morphology, and photographs and odour samples were taken for later re-assessment wherever possible. As a previous study (Xu et al. 2011) had collected a considerable amount of data (particularly on RI) on Exa, Gar and Sph, the sampling of this study prioritized data collection that enabled incorporation of Inc into the data set. All study species co-occur in the Gargano area of southern Italy and coflower in spring (March-May); plants of all species occur in close proximity to each other and well within the expected distance of pollinator-mediated pollen transfer (see Peakall & Schiestl 2004;Xu et al. 2011). Field experiments were performed at populations near Marina di Lesina (MDL; previously used by Xu et al. 2011) and Mattinata (MTT; N41.7349°, E16.1055°; chosen because O. incubacea appeared commoner there). To achieve a more comprehensive genetic sampling, additional populations in the area were sampled, namely San Nicandro Garganico (SNG; N41.8063°, E15.5075°), Cagnano Varano (CGV; N41.8476°, E15.6955°) and Capoiale (CAP; see Xu et al. 2011). The four study species occur in mosaic sympatry (sensu Mallet et al. 2009) at these sites (map in Fig. S2, Supporting information), each being sympatric for at least three of the four study species; in particular, Inc was not observed at populations CAP and SNG, and Exa was absent from CGV and MTT.

Reproductive isolation
Different components of premating and postmating RI between Ophrys species were quantified as detailed in Appendix S3 (sample sizes: Table S1, Supporting information). Briefly, for premating RI, flower phenology and pollinator-mediated floral isolation (and pollination success) were assessed; floral isolation was measured by direct tracking of pollen transfer (see Xu et al. 2011). Assuming absence of other prezygotic barriers, the frequency of interspecific pollen in the pollen pool will depend exclusively upon the relative frequency of flowers of different species. Therefore, following Martin & Willis (2007), pairwise temporal RI indices for the three studied species were calculated based on the frequencies of flowers of the different species during their flowering seasons. For this calculation, we used the spreadsheet provided in Table S3 of Lowry et al. (2008). To evaluate potential postmating barriers, ploidy was measured and crossing experiments conducted to determine seed capsule formation (fruit set) and the proportion of developed (potentially viable) embryos in seeds (as in Scopece et al. 2007;Xu et al. 2011). As orchid female gametophyte development is triggered by pollination (Zhang & O'Neill 1993), fruit set is considered to be an estimate of prezygotic postmating RI. These data were used to (i) estimate temporal and floral RI, and pre-and postzygotic postmating RI and (ii) to calculate total pre-/postmating RI and (iii) the contribution of different barriers to overall RI, as described previously (Scopece et al. 2007).

Phenotypic analysis of floral traits
Floral traits potentially contributing to pollinator-mediated RI, particularly flower size, morphology, colour, odour and speculum shape, were evaluated as described in detail in Appendix S3 (sample sizes in Table S1, Supporting information). In brief, flower size, colour and speculum shape were measured from photographs, and spectral reflectance data were collected and mapped into the honeybee colour-space (as in Chittka & Kevan 2005); three-dimensional (3D) flower morphology was determined by micro-computed tomography (lCT) (see Staedler et al. 2013; Table S2, Supporting information); floral odour was analysed by gas chromatography as described previously (Xu et al. 2011).

Population genomic data generation and analysis
Population genomic data for the study species were obtained using a previously published genotyping-bysequencing (GBS) protocol (Elshire et al. 2011) with modifications as detailed in Appendix S3 (Supporting information). Paired-end Illumina HiSeq 2000 sequence data (two lanes; 32 + 96 samples; Tables S1, S3 and S4, Supporting information) from 127 plant individuals and a replicate were analysed as described in Supporting information. Briefly, after demultiplexing, single nucleotide polymorphisms (SNPs) were identified using STACKS (Catchen et al. 2011), population structure analysed by STRUCTURE (Falush et al. 2003) and principal coordinate analysis (PCoA), and genome scans for F ST outlier loci conducted using FDIST2 and BAYESCAN (Beaumont & Nichols 1996;Foll & Gaggiotti 2008). The latter analysis was conducted both as a global analysis and between species pairs, identifying outlier loci that were repeatedly detected in either two or three sympatric populations. Linkage disequilibrium (LD) was compared between the entire data set and the outlier loci.

Flowering phenology and temporal isolation
Flowering time was recorded in 2012 using a total of 1061 flower observations in the field, the lower number of observations for Inc (Tables S1 and S5, Supporting information) reflecting the lower abundance of this species at the monitored populations. Phenology data from four patches were combined and showed a similar trend as (qualitative) data from 2011. Exa and Sph flowering peaked a week before Gar and Inc, the number of flowers declining thereafter, with the exception of Exa which displayed a 2-week flowering peak. Exa and Sph flowering time overlapped by 81%, Gar and Inc by 94%, although ≥46% of flowers of all species were open simultaneously (Table S6, Supporting information). Corresponding mean RI phenology ranged from 0.025 to 0.5 (Table 1), with an average of 0.26. Temporal RI was weakest between Gar/Inc and Exa/Sph and strongest between Sph/Gar and Sph/Inc.

Floral isolation and pollination success
Average pollination success was markedly different for the species between the years (Fig. S3, Supporting information); pollination success for Sph was higher than for Gar or Inc in 2011, whereas the opposite was true in 2012. In total, 2500 stained flowers were used in pollen tracking plot experiments. Of 108 pollinated flowers (4.32%), the majority received only unstained massulae (from newly opened flowers on the experimental plants or from the surrounding naturally occurring plants). In total, 46 flowers received stained massulae (Table S7, Supporting information), 45 of which indicated withinspecies transfers (15 Sph, 12 Gar, 18 Inc), and one interspecies pollination event (Table S7, Supporting information). Specifically, one Gar flower had received massulae from Inc. The mean floral isolation index (RI floral ) was 0.96 between Gar and Inc and one for all other tested species combinations (Table 1). Likewise, complete floral isolation (RI floral = 1) was previously observed for all species combinations involving Exa, Gar and Sph (Xu et al. 2011).

Ploidy level
All Ophrys species were clearly diploid, consistent with previous reports (D'Emerico et al. 2005;Xu et al. 2011); only a single Sph individual (population CAP) appeared to be triploid (Fig. S4A, Supporting information). A small but significant difference in nuclear DNA A/T content between Inc and Gar was observed, consistent with karyological data (D'Emerico et al. 2005).

Postmating isolation
After hand pollination, 90% of intraspecies (Inc 9 Inc) and 64-94% of interspecies crosses developed capsules (Table S8, Fig. S4B, Supporting information). Similarly, 68% of intra-and 63-76% of interspecies seeds contained embryos (therefore considered viable; Table  S8, Fig. S4C, Supporting information). Although mean fruit set for Inc 9 Gar and Inc 9 Exa appeared comparatively low, this finding was not significant. Neither was there a significant difference in the proportion of viable seeds in any cross. Fruit set and viable seed proportion between intra-and interspecies crosses were used to estimate the strength of postmating prezygotic and postzygotic barriers as described previously (Scopece et al. 2007); negative values were assumed to indicate zero. Postmating isolation was largely weak or absent, although the RI prezygotic values for the aforementioned Inc 9 Gar and Inc 9 Exa crosses were also comparatively high (Table 1). Weak postmating isolation is in line with previous data for other species combinations (Xu et al. 2011), which were re-analysed and included in Table 1 for comparison.

Overall RI
In nature, species may be isolated by several barriers, which work in a sequential manner; each barrier therefore also affects the action of subsequent barriers. Data from this study and Xu et al. (2011) were used to estimate the absolute strength and the contribution of each barrier to total RI (Table 1). The early acting barriers (phenology and floral isolation) were strong and responsible for the majority of total isolation. Combined premating RI was 1 (i.e. complete) for all species comparisons except for Gar and Inc (RI premating = 0.96). Combined postmating RI was generally weak (max. 0.19 for Exa/Gar); however, in all cases, postmating barriers effectively did not contribute to overall RI.

Flower morphology and size
Comparison of flower size revealed significant interspecies differences (Fig. S5, Supporting information): Exa and Gar differed from Inc in total flower width, width of the stigmatic cavity and the position of the widest point of the labellum. Gar and Inc had a wider labellum than Sph and Exa; Exa had the longest flowers; and Sph had a smaller labellum area than Gar and Inc. The largest overall differences were in labellum 'slenderness' (length/width ratio; Fig. S5I, Supporting information), which was highest for Exa, medium for Sph and Inc, and lowest for Gar.  A more detailed view of floral morphology was obtained by lCT scanning of 56 flowers (Figs 1 and S6 to S8; movies: Files S1 to S4, Supporting information). Figure S6 (Supporting information) presents average and 'extreme' reconstructed 3D flowers for each species; the geometric landmarks used for analysis are depicted in Fig. S7A (Supporting information). Canonical variate analysis (CVA) based on these landmark data showed the four species to be clearly morphologically distinct (Fig. 1). Within the reconstructed morphometric shape space of the pooled four-species data set, the first two principal components (PCs) of shape variation both represent relative changes of position of the viscidia with respect to the labellum. Along PC1 (36% variance explained), the viscidia move more or less in parallel to the plane of the labellum, whereas along PC2 (30% variance explained), the angle between column and labellum varies, thereby bringing the viscidia closer to, or farther from the labellum (Fig. S7B, Supporting information). These two relatively simple deformations allow a large flexibility in the positioning of the viscidia with respect to the labellum. A striking detail that differs between Inc and the other species are Ophrys incubacea's longer trichomes on the side of the labellum (Fig. S8, Supporting information; cf. Cortis et al. 2009).

Flower colour
Differences in colour, as perceived by a human observer, were quantified from colour-adjusted flower photographs. In the field, flower labella of Gar and Inc appear darker than those of Sph, and this was reflected in a significant difference in labellum R[ed], G[reen] and B [lue] values (Fig. S9, Supporting information). Likewise, the whitish perigon (see Fig. S1, Supporting information) of Exa, unsurprisingly, constituted a significant difference; other differences were minor. Spectral reflectances of Sph and Inc flower parts were collected and mapped into a colour hexagon representing the bee visual space (Fig.  S10, Supporting information). This analysis suggests that while the dark labella of both species are uncoloured (achromatic) to their pollinators, their specula provide a UV/UV-blue colour signal; perigon parts appeared green.

Speculum shape
The colour difference between speculum and the rest of the labellum implies that speculum shape should be detectable by pollinators. Speculum shape is highly variable (Fig. S11, Supporting information), but elliptic Fourier-descriptor-based PC analysis (EFD-PCA) revealed no clear shape differences between species (Fig. S11B, Supporting information). However, species differed significantly in the proportion of labellum area taken up by the speculum (Fig. S11F, Supporting information; highest: Exa, medium: Gar, lowest: Inc). Speculum pattern complexity (number reflective patches/ nonreflective holes) differed between species, with the largest difference between Sph and Inc (Fig. S11G, H, Supporting information).

Floral odour
Ophrys floral odour has been well-studied because of its importance for pollinator attraction. Floral odours of Exa, Gar and Sph in the study area had previously been investigated (Mant et al. 2005;Xu et al. 2011); our data are in agreement with previous findings. Here, we analysed 24 compounds, comprising six n-alkanes with different chain lengths (C 21 -C 31 ), 16 alkenes with different chain lengths (C 21 -C 31 ) and double-bond positions [(Z) 7-, (Z)9-, and (Z)12-alkenes], and two esters (2-nonyl palmitate and 2-nonyl oleate) (Fig. S12, Supporting information). Amongst these, esters were only detectable in Gar and, at low levels, in Inc (Fig. S12F, Supporting information). Sph had the lowest overall alkene level. Gar, Sph and Inc had higher levels of 9-and 12alkenes than Exa, and conversely, Exa had the highest proportion of 7-alkenes, consistent with previous findings (Schl€ uter et al. 2011b;Xu et al. 2012a), Inc also containing a high proportion of 7-alkenes as compared to Gar/Sph (Fig. 2A). The latter two species had a similar alkene double-bond composition and primarily differed in carbon chain length, Gar featuring the longest hydrocarbons on average and Exa the shortest (Fig. 2B). Linear discriminant analysis (LDA) based on all identified compounds clearly separated the four species into four corresponding clusters (Fig. 2C).

Population genomic data
Genotyping-by-sequencing data were collected for 127 biological individuals, one of which was repeated to assess reproducibility. In total, 141 million assignable paired-end reads were obtained, and 67 752 (raw: 2 538 827) unique sequence tags with 158 871 (raw: 707 059) SNPs identified. Unfortunately, the number of reads per sample was highly uneven (especially in the 96-plex run), and ranged from 43 345 to 14 738 130, with SNPs/sample ranging from 550 to 78 729 (Table  S4, Supporting information). Nevertheless, genotyping repeatability was high for loci sampled twice among technical replicates (99.2%; 63 conflicting genotypes out of 3037 SNPs at coverage ≥10). Overall, 2 823 364 genotype calls were observed at ≥10 read coverage (157-68 306 calls/individual; mean 22 058). Uneven sequence coverage implied a potentially large amount of missing data in STRUCTURE input files (88-119 individuals at 585-4807 loci), the effect of which was investigated by varying the stringency of data filtering. However, best K estimates were highly inconsistent among analyses (K = 2-15) and showed no apparent patterns; the analyses generally suggested a lack of strong genetic structure in the study species. Often, when clusters corresponding to the study species were identified (e.g. Fig. S13B, Supporting information), they only accounted for a small proportion of the genome, whereas the majority of markers were typically shared among species. Exa and Gar were most frequently separated from other species; Inc was frequently split into two groups.
Principal coordinate analysis revealed a weak but discernible clustering according to species identity (Fig. 3), the first three axes however only explaining 19.6% of variation. This analysis qualitatively confirmed STRUC-TURE analysis, including the suggestion of genetic structure within Inc and the presence of a few Sph/Inc individuals that were not clearly assignable to a species. Although improbable, we note that sample misidentification cannot be completely excluded. Like STRUCTURE, mean F ST -based analysis also suggested Inc and Sph to be the most similar, and Exa the most dissimilar study species (Fig. S13, Supporting information), although no in-group relationship had any bootstrap support. Overall, although genetic structure was broadly consistent with species groups, species were very similar and only weakly differentiated.
Outlier analysis (FDIST2) on 95 079 SNP markers typically revealed 9-13 tags that both contained global F ST outliers and repeated outliers in pairwise species comparisons in at least three sympatric populations (21-197 replicated in ≥2 populations) (Figs 4A and S13F, Supporting information). The stringent global BAYESCAN analysis of the same 95 079 markers revealed 19 (0.02%) F ST outlier SNPs (in 17 tags) at FDR <0.05 (Fig. 4B). A significantly elevated extent of LD was found for all sets of outlier loci (all P < 0.0005; Fig. 4C), and LD in pairwise species outliers increased with divergence (Fig. S14, Supporting information). As expected for GBS markers sampled from a large genome, BLAST searches of tags identified by threefold population replication and/or BAYESCAN analysis (Appendices S1 and S2, Supporting information) did not identify confirmed (i.e. not 'predicted') or function-   (Fig. 4D). Both outlier SNPs found in VPS45, like the outlier SNP in CER1, were synonymous. However, the SNP in CER1 was placed at a predicted exon/intron junction (Fig. 4D).

Reproductive barriers
The emergence of RI between populations is critical for speciation. All Ophrys study species were found to be strongly reproductively isolated from each other, as was expected given their different pollinators. Our survey of reproductive barriers demonstrates the importance of premating RI (Table 1), that is temporal and floral RI. Flower phenology differences between early (Exa/Sph) and late-flowering species (Gar/Inc) accounted for considerable temporal RI (Table 1), for instance between the most closely related species pair (Sph/Inc). As efficient gamete transfer requires the presence of pollinators, one may expect pollinatorimposed selection to synchronize flower phenology with pollinator flying times. Interestingly, pollination success was higher in 2011 for the early flowering spe-  cies (Sph), whereas in 2012, the late-flowering species (Gar, Inc) had higher success (Fig. S3, Supporting information), suggesting a difference in the abundances of the respective pollinators. This is consistent with the notion of a spatiotemporal pollinator mosaic (Johnson 2006;Breitkopf et al. 2013).
In absolute terms, the strongest barrier was floral isolation, which was complete (RI floral = 1) in all measured cases (including those published by Xu et al. 2011) except for Gar/Inc, where a single interspecies pollen transfer was observed. Despite this, floral isolation was also very high (RI floral = 0.96) in this case. Floral trait-based RI is mediated by pollinator behaviour, which in Ophrys has been shown to be largely determined by floral odour that acts as a pseudo-pheromone Schl€ uter & Schiestl 2008). The finding of strong floral isolation is entirely consistent with previous data on sexually deceptive orchids (Schl€ uter et al. 2007Peakall et al. 2010;Xu et al. 2011;), but we can now conclude that strong floral isolation is consistent over several seasons (Sph/Gar: 4 years) and locations (cf. Xu et al. 2011). Strong, constant floral isolation supports the role of pollinators in species divergence, and notably, such strong floral isolation may also be found in more generalized pollination systems (Schiestl & Schl€ uter 2009 and references therein). Although floral isolation is clearly the most effective barrier, the finding of considerable temporal RI-which acts before floral RI-suggests that flower phenology may play an important role in speciation in sexually deceptive orchids that merits closer attention (cf. Xu et al. 2012b).
There were only small genome size differences between Gar and Inc, and postmating barriers were effectively absent between species (Table 1, Fig. S4, Supporting information and Xu et al. 2011). However, even with stronger postmating RI, the large premating RI indices would certainly have rendered the contribution of postmating to total RI negligible when species are sympatric. Nonetheless, this also indicates that rare interspecies pollination events are expected to result in gene flow between the species. Overall, our data suggest that early acting premating barriers, but not later-acting barriers, have evolved between the four study species, similar to findings in the distantly related Australian Chiloglottis ) and other plant groups (Rieseberg & Willis 2007;Lowry et al. 2008;Widmer et al. 2009). As less closely related Ophrys species have been shown to build up postmating barriers (Cortis et al. 2009), this supports the idea that premating barriers may generally evolve earlier (Coyne & Orr 1989;Widmer et al. 2009).

Floral trait differences
Floral isolation, the strongest barrier between the study species, is mediated by floral traits. We investigated a suite of floral traits that may potentially be involved in floral isolation for differences between species. Given the apparent variability of Ophrys flowers in the field (Fig. S1, Supporting information), previous investigations of our study species (e.g. Mant et al. 2005;Xu et al. 2011) largely ignored variation in floral shape and morphology. Study species showed only minor floral size differences; labellum length may contribute to floral isolation because of the need to match pollinator body size (Paulus & Gack 1990;Paulus 2006;Vereecken 2009); however, no significant labellum length differences were found between species. Nonetheless, 3D morphometric analysis revealed a clear separation of species in overall flower shape and suggests that an important component of shape variation is the positioning of the viscidia with respect to the labellum, thereby potentially affecting placement of pollinia on the pollinators (Fig.  S7B, Supporting information). While in Ophrys incubacea the lateral protrusions of the labellum are conspicuously developed and trichomes are longer than in other species (Fig. S8, Supporting information), the functional significance of these features remains unclear. Although trichome direction has been implied in mechanical floral isolation between Ophrys groups pollinated by the same insects by pollinaria attached to either the head or the abdomen (e.g. Agren et al. 1984), this may constitute only a relatively weak barrier (Cortis et al. 2009). Likewise, although potentially preventing interspecies pollination, slight differences in pollinarium placement would be expected to contribute only modestly to overall floral isolation.
Colour cues may aid specific pollinator attraction and thereby contribute to floral isolation. One may expect flowers (i) to share features of the mimicked female insect and (ii) to be distinguishable from other species by their pollinators. Among the three Andrena-pollinated study species (Gar, Inc, Sph), labellum coloration of Inc and Gar was markedly blacker than of Sph. This intriguingly mirrors their pollinators' black (Inc, Gar) and brownish (Sph) body colour, suggesting that this trait may be under divergent selection by pollinators. As labella of both Inc and Sph appear achromatic to pollinators (Fig.  S10E, Supporting information), we hypothesize that pollinator-mediated selection on labellum coloration acts on the level of brightness rather than colour (cf. Renoult et al. 2013). A coloured perigon has been shown to increase pollinator visitation in Ophrys heldreichii (e.g. Spaethe et al. 2007;Streinzer et al. 2010). However, this is less likely in the species studied here, given that no effect of perigon coloration on attractiveness to pollinators was found in (the more closely related) Ophrys arachnitiformis ). In contrast to overall labellum coloration, the speculum likely constitutes a visual colour signal (Fig. S10E, Supporting information), implying that its highly variable shape is detectable by pollinators against the uncoloured labellum background. No clear interspecies differences in speculum shape, however, were found between species, suggesting that speculum shape does not aid the between-species discrimination of flowers by pollinators. However, large between-individual/within-species shape variation may well help pollinators to learn to avoid individual 'false females' they have already visited, thereby promoting outcrossing . Although species differences in speculum complexity and relative area were detected, their functional significance remains unclear. Taken together, labellum coloration may well reflect selection by pollinators, but overall evidence for the contribution of colour signals to floral isolation is limited.
Floral odour underlies the chemical mimicry of pollinators' sex pheromones in sexually deceptive orchids and thus is a key contributor to floral isolation Schiestl & Schl€ uter 2009;St€ okl et al. 2009;Ayasse et al. 2011;Xu et al. 2011Xu et al. , 2012b. Clear differences in floral odour were observed, the study species differing in carbon chain length and alkene double-bond position (Fig. 2). This is consistent with previous data documenting floral odour differentiation between Exa, Gar and Sph (Mant et al. 2005;Xu et al. 2011). Similarly, O. incubacea produces a distinct odour blend, which interestingly includes 7-alkenes (only present in small amounts in Gar and Sph), albeit at lower levels than Exa. Such alkene double-bond differences have previously been linked to changes at desaturase-encoding genes (Schl€ uter et al. 2011b;Xu et al. 2012a). Besides hydrocarbons, two esters were included in our odour analysis because of their high levels in Gar. These two esters were also present at low levels in Inc but undetectable in Sph and Exa. Unfortunately, the sex pheromone of Ophrys garganica's pollinator, Andrena pilipes, has not yet been determined and the function of these esters requires experimental testing.
Overall, several traits may contribute to floral isolation, but phenotypic differences suggest a role of floral odour, morphology and perhaps even coloration. Among these, based upon previous data from our study group (e.g. Schiestl et al. 2000;Mant et al. 2005;Xu et al. 2012a), floral odour is certainly the trait with the largest contribution and with pronounced species differences. This is consistent with findings in Chiloglottis, where, however, no clear morphological differences were found . Hence, one may hypothesize that changes in floral odour precede changes in other floral traits during species divergence. Likewise, one may therefore expect genetic changes affecting floral isolation (odour) and/or temporal isolation (flowering time) to initiate species divergence. This prediction is not restricted to the study species but may more broadly be applied to other plants which diverge via the evolution of strong premating barriers.

Genomics of species divergence
A genome-wide set of SNP markers revealed all four study species to be only weakly differentiated, with the majority of allelic variation shared among all species. Weak genetic interspecific differentiation mirrors earlier population genetic reports that utilized a small number of microsatellite or AFLP markers (Soliva & Widmer 2003;Mant et al. 2005;Xu et al. 2011;Breitkopf et al. 2013). Although initially interpreted as high levels of gene flow among distinct species in secondary contact (Soliva & Widmer 2003), those genetic results may also be explained by shared ancestral polymorphism and recent species divergence (Mant et al. 2005;Schl€ uter et al. 2011a,b;Xu et al. 2011;. More evidence now points in the direction of recent and perhaps ongoing speciation, specifically: (i) phylogenetic dating analysis suggests that the species group has indeed only recently diversified in the Pleistocene (Inda et al. 2012;Breitkopf et al. in press), (ii) the measurement of consistently strong RI in the field makes rampant gene flow unlikely, (iii) ancestral polymorphism has been documented in Ophrys (Schl€ uter et al. 2011a) and (iv) large effective population sizes, together with plant longevity (including clonal vegetative reproduction), make the maintenance of low nuclear differentiation plausible (Mant et al. 2005;. Despite the large amount of shared polymorphism, the study species could largely be distinguished genetically, the separation of Inc/Sph being least clear. The latter is inferred to be the most closely related species pair, with the caveat that there was essentially no statistical support for the inference of species relationships. Only a very small proportion of the genome (<0.05%) was repeatedly identified as more strongly differentiated between pairs of species than expected and is therefore interpreted as being associated with species divergence. These outlier loci may be expected to be linked to targets of divergent selection, for example by pollinators. A caveat to this analysis is that the GBS data only cover a small fraction of the~10 Gbp (Leitch et al. 2009) genome (on average, one analysed tag every <1 kbp); however, it should be noted that the GBS method biases against sampling of methylated, highly repetitive parts of the genome (Elshire et al. 2011). A second caveat is that moderately divergent parts of the genome might not be detected as outliers and that we may therefore underestimate the true extent of genomic differentiation that is linked to species divergence. It is, however, surprising that differentiation between species is so low, given that they differ in a number of phenotypic traits and consequently would be expected to differ at a number of genetic loci involved in specifying these traits. Moreover, considerable genetic structure has previously been reported for other closely related Ophrys species based on fewer markers (Schl€ uter et al. 2011a), making it unlikely that significant genome-wide differentiation has gone undetected in our study. It is currently unclear whether several genes involved in phenotypic differences are closely linked in the genome, for instance due to variation in chromosome structure or divergence hitchhiking. Alternatively, it is also conceivable that indeed only few genetic polymorphisms differ between species if change in several traits (many of which concern features of epidermal cells) is affected by variation in, for example, a key transcriptional regulator controlling phenotypic differences at the level of gene expression. Nonetheless, the small proportion of F ST outliers suggests that the Ophrys species studied here have only reached an early phase of the speciation continuum, that is, either direct selection on genic targets or at best the stage of divergence hitchhiking. Moreover, elevated LD among F ST outlier loci indicates that the few outlier loci identified probably stem from even fewer divergent genomic regions, among-outlier LD increasing with species divergence. The small number of outlier regions therefore certainly suggests genic rather than genome-wide patterns of divergence between Ophrys species and notably fits the model of species divergence proposed for other sexually deceptive orchids .
Genic targets of pollinator-mediated selection would be expected to include genes involved in specifying floral odour-based RI. Such odour changes may only require few genetic changes, as exemplified by the gene SAD2 (not covered by GBS data) involved in alkene differences between Exa and Sph (Schl€ uter et al. 2011b;Xu et al. 2012a). The product of this gene catalyses specific double-bond insertion into precursors of alkenes, and SAD2 alleles with different enzymatic activities are differentially expressed between the two species; SAD2 thereby contributes to a major difference in hydrocarbon composition between Exa and Sph/Gar that demonstrably affects pollinator attraction and thus RI (Schl€ uter et al. 2011b;Xu et al. 2012a). In principle, other genes affecting pseudo-pheromone production (see Sedeek et al. 2013) may also be expected to respond to divergent selection by pollinators. Here, two annotated genes, VPS45 and CER1, were identified among the GBS tags containing F ST outliers. Arabidopsis VPS45 is an SM protein family member involved in vesicle trafficking and regulation of a SNARE complex localized on the trans-Golgi network (TGN) that is required for membrane fusion (Zouhar et al. 2009). Interestingly, secretion of alkanes to the plant surface by epidermal cells has recently been found to require Golgi-and TGN-mediated vesicle trafficking (McFarlane et al. 2014), thereby raising the possibility that VPS45 is involved in producing the adaptive, pollinator-attractive odour phenotype. Although the outlier SNPs in VPS45 were synonymous, evidence for selection on, and functional relevance of, synonymous sites (e.g. for miRNA binding) is mounting (Gu et al. 2012;Lawrie et al. 2013). It may, however, be more likely that these SNPs are linked to other polymorphisms of functional importance. The finding of an F ST outlier in CER1 is especially striking, considering that this gene has been considered an a priori candidate for affecting pollinator attraction (Sedeek et al. 2013). Arabidopsis CER1 is an aldehyde decarbonylase acting in the last step of hydrocarbon biosynthesis (Bernard et al. 2012); the outlier SNP at the Ophrys CER1 exon/intron junction may well have functional consequences by affecting pre-mRNA splicing, which may in turn result in pollination-relevant phenotypic change. It must be borne in mind that hypotheses on gene function based upon homology to Arabidopsis are inherently uncertain and, even if correct, that the outlier SNPs identified here may simply be linked to other, more important, polymorphisms. Although intriguing, the involvements of VPS45 and CER1 in pollinator-mediated species differences therefore require testing. Nevertheless, the fact that one (or perhaps even both) of the two outlier-containing genic regions has putative functions linked to odour production suggests that the outliers identified are not spurious and that these genic regions might represent direct targets of selection, as would be expected in the early genic phase of species divergence.

Conclusion
This study examines pollinator-driven ecological speciation between sexually deceptive orchids with a specialized pollination system. Species are isolated by premating barriers, namely temporal isolation and, especially, strong floral isolation. Floral odour and, to a lesser extent, flower geometry and perhaps labellum coloration probably underlie floral isolation. Population genomic data revealed the majority of polymorphism to be shared across species. Only few repeatable F ST outliers were found in a genome scan, two genes with potential involvement in odour production being among the outliers. Overall, these data suggest that these species are defined mostly by genic rather than genome-wide differences, suggesting that Ophrys orchids have only reached an early stage in the speciation continuum. We hypothesize that this also applies to other sexually deceptive orchids, where odour genes may be among the first to diverge. More generally, genic speciation may well be common in systems where subtle trait changes strongly affect prezygotic RI.

Supporting information
Additional supporting information may be found in the online version of this article.
Appendix S1 FASTA file of tag sequences containing the best outlier SNPs in the global BAYESCAN F ST outlier analysis.
Appendix S2 FASTA file of tag sequences containing the best outlier SNPs from pairwise FDIST2 F ST outlier analyses, replicated in three populations.