# Contrasting patterns of differentiation among three taxa of the rapidly diversifying orchid genus Ophrys sect. Insectifera (Orchidaceae) where their ranges overlap

Pascaline Salvado, Anaïs Gibert, Bertrand Schatz, Lucas Vandenabeele, Roselyne Buscail, David Vilasís, Philippe Feldmann, Joris Bertrand

► To cite this version:

Pascaline Salvado, Anaïs Gibert, Bertrand Schatz, Lucas Vandenabeele, Roselyne Buscail, et al. Contrasting patterns of differentiation among three taxa of the rapidly diversifying orchid genus Ophrys sect. Insectifera (Orchidaceae) where their ranges overlap. Botanical Journal of the Linnean Society, 2025, 208 (2), pp.142-151. 10.1093/botlinnean/boae062. hal-04893550

HAL Id: hal-04893550

https://hal.science/hal-04893550v1

Submitted on 17 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Original Article

Contrasting patterns of differentiation among three taxa of
the rapidly diversifying orchid genus *Ophrys* sect. *Insectifera*
(Orchidaceae) where their ranges overlap

Pascaline Salvado<fnref n="1" />,*, ID, Anaïs Gibert<fnref n="1" />, ID, Bertrand Schatz<fnref n="2" />, ID, Lucas Vandenabeele<fnref n="1" />, ID,
Roselyne Buscail<fnref n="3" />, ID, David Vilasís<fnref n="4" />, Philippe Feldmann<fnref n="5" />, ID and Joris A. M. Bertrand<fnref n="1" />, ID

ABSTRACT

In rapidly diversifying groups, taxa defined on the basis of typological criteria can be difficult to support with genetic data. The diversity observed in the insect-mimicking orchid genus *Ophrys* perfectly illustrates this situation; among 400 described species only 9–10 lineages are detectable by genetic markers such as nrITS. The three taxa described in the *Ophrys insectifera* group, *O. insectifera*, *O. subinsectifera*, and *O. aymoninii*, can be clearly distinguished by their flowers, which have evolved different phenotypes as a result of adaptation to specific pollinator insect species from three different families. However, genetic differentiation between these three taxa has never been really supported by population genetic data and their taxonomic status is still debated. Using population genomic approaches, we found a clustering consistent with the existence of three genetic entities where the geographic distributions of the three taxa overlap. Two of these clusters correspond to France populations of the widespread *O. insectifera* and the micro-endemic *O. aymoninii*. However, the last cluster grouped together all the Iberian individuals, suggesting that individuals phenotypically identified as either *O. insectifera* or *O. subinsectifera* are genetically weakly differentiated there. Populations of the two pairs of taxa thus may have experienced different patterns of inter-specific gene flow.

**Keywords:** biogeography; genotyping by sequencing; inter-specific gene flow; morphometrics; plants; population genomics; RAD-sequencing; speciation; systematics; taxonomy

INTRODUCTION

Adaptive radiations can lead to rapid accumulation of species
through adaptation to various selection pressures. In such cases,
reconstructing the evolutionary history based solely on gen-
etic data can be challenging. The advent of high-throughput
sequencing technologies has enabled to work at the whole-
genome scale, thus providing significant gain in the resolution
of genetic data, even in non-model organisms (e.g. Léveillé-
Bourret et al. 2020). For species with large and/or complex gen-
omes however, genome complexity reduction protocols remain
the most efficient way to simultaneously obtain genomic data on

a sufficient number of individuals and populations to investigate
how evolutionary lineages diverge from each other.

*Ophrys* orchids are renowned for their exceptional ability to attract male insect pollinators through their insect-mimicking flowers. The labellum, one of the petals, is especially effective in luring the insects to the point where they attempt to copulate with it, ultimately ensuring the plant's reproduction. The degree of specification required for the emergence and the maintenance of this plant-insect interaction is likely to be responsible for the high rate of diversification observed in *Ophrys* orchids [see Baguette *et al.* (2020) for a recent review]. The flower

<footnote n="1">¹Laboratoire Génome et Développement des Plantes (LGDP), UMR 5096, Université de Perpignan Via Domitia (UPVD) – Centre National de la Recherche Scientifique (CNRS) – Institut de Recherche pour le Développement (IRD), EMR 269 MANGO, Perpignan, France</footnote>

<footnote n="2">²CEFE, CNRS, Univ Montpellier, EPHE, IRD, Montpellier, France</footnote>

<footnote n="3">³Centre de Formation et de Recherche sur les Environnements Méditerranéens (CEFREM), Université de Perpignan Via Domitia (UPVD) – Centre National de la Recherche Scientifique (CNRS), Perpignan, France</footnote>

<footnote n="4">⁴Grup Oracidològic de Catalunya GOC HCHN, Barcelona, Spain</footnote>

<footnote n="5">⁵Association Abiome, Olivet, France</footnote>

<footnote>*Corresponding author. Pascaline Salvado and Joris Bertrand, 58 Avenue Paul Alduy, Bâtiment T, F-66860 Perpignan Cedex 09, France: Tél: +33.(0)4.68.66.17.74; E-mail: pascaline.salvado@univ-perp.fr; joris.bertrand@univ-perp.fr</footnote>
phenotypes of *Ophrys* orchids are subject to strong selection pressures exerted by male pollinator insects, which is one of the main drivers of this extreme adaptive radiation. The morphological diversity observable among *Ophrys* species has been well-documented and an increasing number of studies have reported subtle variations in flower odour compounds associated with pollinator species identity (e.g. Schlüter *et al.* 2011, Joffard, Buatois and Schatz 2016, Gervasi *et al.* 2017). However, genetic studies have often failed to characterize patterns and propose relevant scenarios to explain evolutionary divergence between closely related *Ophrys* species.

The inability of traditional genetic markers to delineate *Ophrys* lineages can be explained by divergence times that may be too recent to display allele sorting through genetic drift (i.e. incomplete lineage sorting) and/or reticulate evolution that may blur any emerging phylogenetic signal as well as population structure that may mislead the interpretation of the pattern inferred. In spite of the promising results from a pioneering study using a genotyping-by-sequencing (GBS) protocol to investigate patterns of differentiation between closely related species of the *Ophrys sphegodes* group (Sedeek *et al.* 2014), population genomic approaches have not been used in *Ophrys* since. Bateman *et al.* (2018) have successfully used a RADseq-like protocol to clarify the phylogenetic relationships between the main *Ophrys* lineages, see also Piñeiro Fernández *et al.* (2019) who used RNAseq data and Bertrand *et al.* (2021) who employed whole plastomes. However, none of these studies addressed patterns of species differentiation at lower taxonomic scale (but see Bateman *et al.* 2021). Thus, the actual taxonomic status of the so-called 'microspecies' of *Ophrys* is still regularly contested (e.g. Bateman and Rudall 2023) although the suitability of population genomic approaches to investigate evolutionary divergence between *Ophrys* lineages has rarely been properly tested. It was recently confirmed that a population-based sampling design coupled with an appropriate genome subsampling protocol (i.e. a GBS protocol called nGBS, which stands for normalized GBS) allowed to accurately characterize patterns of diversity and differentiation among little divergent populations and taxa of *Ophrys aveyronensis* (Gibert *et al.* 2024).

In this study, we used a protocol similar to the one of Gibert *et al.* 2024 to test whether this approach also provides sufficient resolution to distinguish taxa within another *Ophrys* clade whose systematics and taxonomy are problematic: the *Ophrys insectifera* group. The latter consists of three described taxa: the Euro-Mediterranean widespread *Ophrys insectifera* L. and two endemic taxa: *Ophrys subinsectifera* C.E. Hermos. & Sabando, whose repartition spans the Southern flanks of the Pyrenees and *Ophrys aymoninii* (Breistr.) Buttler, which is (micro-)endemic to the Grand Causses region in France. The three taxa, currently recognized as species or subspecies, (*Ophrys insectifera* subsp. *subinsectifera* (C.E. Hermos. & Sabando) (O. Bolòs & Vigo) and *Ophrys insectifera* subsp. *aymoninii* Breistr.) are known to be pollinated by insect species belonging to three different families. Moreover, these taxa can be easily distinguished in the field on the basis of their phenotypes and have been shown to have different odour bouquets. So far, however, genetic studies have failed to delineate these taxa based on genetic data (Triponez *et al.* 2013, Gervasi *et al.* 2017). Previous phylogenetic studies only included one or two individuals per taxon often with only two of

the three taxa of the group, which did not allow to properly test their reciprocal monophyly (Soliva and Widmer 2003, Devey *et al.* 2008, Breitkopf *et al.* 2015, Bateman *et al.* 2018). With the future aim of investigating the ecological and evolutionary causes of species boundaries in this system, we tested whether population genomic data could allow to characterize species, especially where *O. insectifera* co-occur with one of the two endemics. Specifically, we tested whether population genomics allows us to distinguish the three taxa despite incomplete lineage sorting and possible inter-specific gene flow. If so, we planned to assess, for the first time, levels of genetic diversity within taxa and genetic differentiation between taxa based on genomic data.

# MATERIAL AND METHODS

## Population sampling

We sampled a total of 239 individuals at 12 localities in May 2022 and 2023, where the widespread *O. insectifera* locally co-occurs with one of the endemic taxa: *O. subinsectifera* in Catalonia, Spain at six sites (1, Les Guilleries; 2, Puig Grifó; 3, Tarradell; 4, Castanyola; 5, L'Estanyol; and 6, Seva) and *O. aymoninii* in the Grands Causses Region, France at four sites (8, Le Buffre; 9, Col de Montmirat; 10, Nivoliers; and 11, Fretma) as well as at two additional localities in between these two regions where only *O. insectifera* is present (7, Versols-et-Lapeyre and 12, Rodome), as represented in Figure 1 (see also Table 1 and Table S1 in the Supporting Information). At each locality, we collected ~1–2 cm² of leaf tissues that were stored in 90% ethanol at 4°C until DNA extraction. *O. subinsectifera* and *O. insectifera* are not protected in the study area but *O. aymoninii* has a Near Threatened status on the French National Red List. We asked for and obtained special permits from the 'Parc National des Cévennes' to collect plants and drive within the national park (2022-0115, 2022-0117, 2023-0093, and 2023-0094). The research project also complies with the Nagoya protocol (TREL2206915S/562), which provides a framework for the fair and equitable sharing of benefits arising from the use of genetic resources. We photographed each individual and measured values over 10 phenotypic traits, as described in detail in Gibert *et al.* (2022). We recorded five 'whole-plant' traits: plant size, stem diameter, number of reproductive structures (flowers and buds), distance from the ground to the first flower, and distance between the first and the second flower. We also measured five floral traits: length and width of the labellum and the median sepal as well as the length of the right petal. We calculated the ratio between the length and width of the labellum to get an overall idea of its proportions. All the individuals were categorized as *Ophrys insectifera*, *O. aymoninii*, *O. subinsectifera* or intermediates (*O. aymoninii*/*O. insectifera* or *O. subinsectifera*/*O. insectifera*) based on established visual criteria (i.e. following Delforge 2021). We performed linear mixed models analysis on each trait to test for difference between the five categories. The fixed effect included the categories, and the random effect included the study site. We presented 95% confidence intervals, and used their lack of overlap with zero or lack of overlap between categories to define what effects were statistically significant, in addition to *P* values. This is a conservative approach that provides information on the statistical effect (range, direction, strength and reliability) that is not provided by *P* values (Ho *et al.* 2019). Models were implemented in R via the
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Map displaying the 12 study sites: four in France for *O. aymoninii*, *O. insectifera*, and their intermediates (8, Le Buffre; 9, Col de Montmirat; 10, Nivoliers; and 11, Fretma); two ‘control’ sites in France with only *O. insectifera* (7, Versols-et-Lapeyre and 12, Rodome); and six sites in Spain for *O. subinsectifera*, *O. insectifera*, and their intermediates (1, Les Guilleries; 2, Puig Grifó; 3, Tarradell; 4, Castanyola; 5, L’Estanyol; and 6, Seva). Pie charts represent the species composition (on the basis of identification on the basis of phenotypes) of each site, yellow corresponding to *O. aymoninii*, orange to *O. insectifera*, dark blue to *O. subinsectifera*, red for the intermediates *aym/ins*, and light blue for the intermediates *sub/ins*.</figcaption></figure>

<table><caption>**Table 1.** Geographic coordinates (in °) and elevation (in m above sea level) of the sampling localities, sample size (*N*), observed (*H*<sub>O</sub>), and expected (*H*<sub>E</sub>) levels of heterozygosity, deviation from panmixia (*G*<sub>IS</sub>), total number of polymorphic sites (*A*), and number of private alleles (*Ap*) per locality.</caption><thead><tr><th>Locality</th><th>Latitude</th><th>Longitude</th><th>Elevation</th><th>N</th><th>H<sub>O</sub></th><th>H<sub>E</sub></th><th>G<sub>IS</sub></th><th>A</th><th>Ap</th></tr></thead><tbody><tr><td>1, Les Guilleries</td><td>41.883</td><td>2.356</td><td>602</td><td>30</td><td>0.228</td><td>0.254</td><td>0.103</td><td>11158</td><td>2</td></tr><tr><td>2, Puig Grifó</td><td>41.857</td><td>2.320</td><td>786</td><td>20</td><td>0.226</td><td>0.249</td><td>0.096</td><td>11034</td><td>0</td></tr><tr><td>3, Tarradell</td><td>41.876</td><td>2.300</td><td>681</td><td>7</td><td>0.229</td><td>0.224</td><td>-0.023</td><td>10750</td><td>0</td></tr><tr><td>4, Castanyola</td><td>41.855</td><td>2.276</td><td>604</td><td>15</td><td>0.233</td><td>0.247</td><td>0.058</td><td>10980</td><td>0</td></tr><tr><td>5, Seva</td><td>41.824</td><td>2.277</td><td>742</td><td>21</td><td>0.232</td><td>0.245</td><td>0.052</td><td>10756</td><td>0</td></tr><tr><td>6, L'Estanyol</td><td>41.831</td><td>2.284</td><td>715</td><td>24</td><td>0.236</td><td>0.244</td><td>0.036</td><td>10890</td><td>1</td></tr><tr><td>7, Versols-et-Lapeyre</td><td>43.897</td><td>2.934</td><td>412</td><td>17</td><td>0.232</td><td>0.267</td><td>0.131</td><td>11066</td><td>0</td></tr><tr><td>8, Le Buffre</td><td>44.292</td><td>3.411</td><td>922</td><td>40</td><td>0.233</td><td>0.272</td><td>0.143</td><td>11182</td><td>0</td></tr><tr><td>9, Col de Montmirat</td><td>44.419</td><td>3.551</td><td>1068</td><td>24</td><td>0.226</td><td>0.253</td><td>0.105</td><td>10846</td><td>0</td></tr><tr><td>10, Nivoliers</td><td>44.255</td><td>3.453</td><td>1017</td><td>14</td><td>0.225</td><td>0.266</td><td>0.155</td><td>10679</td><td>0</td></tr><tr><td>11, Fretma</td><td>44.242</td><td>3.389</td><td>1117</td><td>15</td><td>0.235</td><td>0.232</td><td>-0.013</td><td>10868</td><td>4</td></tr><tr><td>12, Rodome</td><td>42.800</td><td>2.049</td><td>959</td><td>12</td><td>0.241</td><td>0.267</td><td>0.095</td><td>10825</td><td>8</td></tr><tr><td>Mean</td><td></td><td></td><td></td><td></td><td><strong>0.231</strong></td><td><strong>0.252</strong></td><td><strong>0.078</strong></td><td></td><td></td></tr></tbody></table>

'lme4' package (Bates et al. 2015). Assumptions of models (e.g. heteroscedasticity, over and under dispersion, significance of outliers) were checked both visually and by using the DHARMa R-package (Hartig 2018).

## Molecular procedures

Genomic DNA extraction and genotyping were subcon- tracted to LGC Genomics GmbH (Berlin, Germany) fol- lowing a protocol called nGBS for normalized GBS. This

double digest restriction site-associated DNA seq or RAD- seq-like protocol relies on two restriction enzymes (PstI and ApeKI in our case) to reduce *Ophrys* genome complexity and it includes a normalization step that aims at avoiding the overrepresentation of repetitive regions. Most *Ophrys* spe- cies are diploid: 2n = 36 (Turco et al. 2023). If preliminary estimates of genome size we obtained from flow cytometry techniques corroborate this ploidy level for *O. insectifera* and *O. aymoninii* (see also Gervasi et al. 2017), the values we
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Distribution (violin plots) of six floral and five whole-plant traits for each species: *O. aymoninii*, *O. insectifera*, *O. subinsectifera*, and individuals with intermediate phenotypes: *aym/ins* and *sub/ins*. The different sites are displayed in different colours. All the values are in mm except for the plant size, distance to first flower, and distance between first and second flower, which are in cm. Asterisks indicate that the values are significantly different from the generalist *O. insectifera*.</figcaption></figure>

measured for *O. aymoninii* (1C value = 6.56 pg, *n* = 1) and *O. insectifera* (1C values = 6.51 pg and 6.83 pg, *n* = 2) were consistently higher than the ones documented in the literature for other diploid *Ophrys* species (Gibert et al. 2024, Russo et al. 2024, see also Bou Dagher-Kharrat et al. 2013). This suggests that genome size could be ~20% larger in the *insectifera* clade compared to representatives of the *O. sphegodes* group. So far, no measures have been done for *O. subinsectifera*. The resulting 239 individually barcoded libraries were sequenced in paired-end mode (2 × 150 bp) on an Illumina NovaSeq 6000 with an expectation of a minimum number of 1.5 million read pairs per sample. We used Stacks v.2.62 (Catchen et al. 2011, 2013) to build loci from Illumina reads. After demultiplexing and cleaning the samples with *process_radtags*, we aligned them to the reference genome of *Ophrys sphegodes*, which is the most closely related reference genome available (Russo et al. 2024) using BWA-mem v.0.7.17 (Li and Durbin 2009). Then, the *ref_map.pl* script was executed for building RAD-loci and SNP calling. Finally, we ran the populations script included in Stacks, from which only loci that were genotyped in at least 80% of the individuals per population (-r 80) and 10 populations (out of 12) were kept (-p 10), while loci exhibiting estimates of observed heterozygosity greater than 70% (--max-obs-het 0.7) were filtered out to reduce the risk of including remaining paralogs. We also discarded sites whose minor allele frequency was <5% (--min-maf 0.05). The genotype matrix of the remaining 11 089 SNPs was then exported in several formats such as Structure, VCF, genepop, and radpainter.

## Population genomics analyses

To get an overview of the overall genetic diversity and differentiation among individuals and populations, we first performed a principal component analysis (PCA) based on a matrix of 239 individuals (as rows) and 11 089 SNPs (as columns) coded as a genlight object with the R-package 'adegenet' (Jombart and Ahmed 2011). We then used Genodive v.3.04 (Meirmans 2020) to compute expected and observed heterozygosity ($H_E$ and $H_O$ respectively). For all the analyses realized in Genodive, we kept only the most abundant species at each site. The number of polymorphic sites (A) and private alleles ($A_p$) were retrieved from the populations outputs given by Stacks. The deviation from panmixia was evaluated by computing $G_{IS}$ and the statistical significance of the obtained values was estimated based on 10 000 permutations, in Genodive. Overall and pairwise genetic differentiation was assessed based on $G_{ST'}$ (Nei 1987) as implemented in Genodive in a similar manner. Genodive was also used to carry out a hierarchical Analysis of Molecular Variance (AMOVA) to assess proportions of genetic variance (i) among series (here, Series A between *O. aymoninii* and *O. insectifera* and Series B between *O. subinsectifera* and *O. insectifera*) ($F_{CT}$) and (ii) among populations within series ($F_{SC}$). The statistical significance of the obtained values was estimated based on 10 000 permutations. To further investigate population structure and characterize putative migration/admixture event, we used ADMIXTURE (Alexande, Novembre and Lange 2009). This method computes 'ancestry coefficients' that represent the proportions of an individual genome that originate from multiple ancestral gene pools. The number of genetic clusters was varied from *K* = 1 to
12, analyses were run with 10 replicates and a cross-validation procedure was performed to determine the optimal value of K. To assess genetic co-ancestry among individuals within and across the different sites, we used haplotype-based population inference approach implemented in *fineRADstructure* v.0.3.1 (Malinsky et al. 2018) using the same filtering criteria described for the SNP dataset. Analyses were performed using default settings and results were visualized with dedicated R-scripts provided with the program.

# RESULTS

## Species distribution

In both France and Spain, we sampled sites with only the endemics (*O. aymoninii* or *O. subinsectifera*), others with only the widespread *O. insectifera*, and also sites with both one of the endemics and *O. insectifera* in close proximity (Fig. 1). In Spain, sites 3, Tarradell; 4, Castanyola; and 5, L'Estanyol showed the presence of *O. insectifera*, *O. subinsectifera*, and a non-negligible number of individuals with intermediate phenotypes (15 out of 117, 13%), which could not be clearly assigned to one species or the other. Conversely, sites 1, Les Guilleries; 2, Puig Grifó; and 6, Seva exclusively presented individuals with phenotypes corresponding to *O. subinsectifera*. In France, our survey revealed three sites (8, Le Buffre, 9, Col de Montmirat, and 11, Fretma) where *O. aymoninii* and *O. insectifera* co-occurred. However, we found only a few individuals with intermediate phenotypes (3 out of 123, 2.4%) at sites 8, Le Buffre and 9, Col de Montmirat (we did not observe *O. insectifera* at site 10, Nivoliers). At sites 7, Versols-et-Lapeyre and 12, Rodome, we observed only *O. insectifera*, which makes sense as these two sites are outside the known range of both endemic species. Therefore, these sites were considered as unadmixed with either one of the endemics and served as reference sites for *O. insectifera* in our study.

## Patterns of phenotypic variation

For morphometric analyses, we considered *n* = 50 *O. insectifera* individuals from sites 3, Tarradell; 5, L'Estanyol; 7, Versols-et-Lapeyre; 8, Le Buffre; 9, Col de Montmirat; 11, Fretma; and 12, Rodome, *n* = 95 *O. subinsectifera* individuals from sites 1, Les Guilleries; 2, Puig Grifó; 4, Castanyola; 5, L'Estanyol; and 6, Seva, *n* = 75 *O. aymoninii* individuals from sites 8, Le Buffre; 9, Col de Montmirat; 10, Nivoliers; and 11, Fretma, *n* = 14 intermediate individuals between *O. insectifera* and *O. subinsectifera* from sites 3, Tarradell; 4, Castanyola; and 5, L'Estanyol and *n* = 3 intermediate individuals between *O. insectifera* and *O. aymoninii* from sites 8, Le Buffre and 9, Col de Montmirat.

Our results confirm that the flowers of the three taxa show differences in morphometrics, with the endemics (*O. aymoninii* and *O. subinsectifera*) being the smallest overall. The labellum length of the two endemic species (*O. aymoninii* and *O. subinsectifera*) is significantly shorter (9.4 ± 1.2 and 10.6 ± 1.3 mm, respectively, mean ± SD, *P* < 0.01, Fig. 2A, Supporting Information Fig. S2A and Table S2A) compared to *O. insectifera* (13 ± 1.7 mm, Fig. 2A, Supporting Information Fig. S2A and Table S2A). The labellum ratio differs between the three species; *O. aymoninii* showed a ratio 42% smaller than *O. insectifera* (respectively, ratio = 1.4 vs. 2.4) and *O. subinsectifera* showed a ratio 15% smaller than *O. insectifera*.

(2, Fig. 2C, Supporting Information Fig. S2C and Table S2C). Petal length is significantly shorter for *O. subinsectifera* (2.7 ± 0.9 mm, *P* < 0.01) and *O. aymoninii* has the shortest sepal length (5.0 ± 0.8 mm, *P* < 0.01). *O. subinsectifera* has a narrower sepal width compared to *O. insectifera* (2.3 ± 0.6 mm, *P* < 0.01, Fig. 2E, Supporting Information Fig. S2E and Table S2E). Looking at whole-plant traits, we observe that both endemics tend to have smaller plant size compared to *O. insectifera* (19.4 ± 6.4 cm, *p* = 0.01 for *O. subinsectifera* and 19.7 ± 6.3 cm, *P* < 0.01 for *O. aymoninii* vs. 25.4 ± 5.9 cm for *O. insectifera*, Fig. 2G, Supporting Information Fig. S2G and Table S2G). In addition, endemics have more compact inflorescences, whereas those of *O. insectifera* appear looser, with greater distances between the ground and the first flower (16.8 ± 4.9 cm, *P* < 0.01, Fig. 2J, Supporting Information Fig. S2J and Table S2J) and between the first and second flower (2.9 ± 0.8 cm, NS, Fig. 2K, Supporting Information Fig. S2K and Table S2K). Intermediate individuals, termed ‘sub/ins’ and ‘aym/ins’, show intermediate values for labellum length, labellum ratio, petal length, and sepal width (Fig. 2C, E, F, Supporting Information Fig. S2A, C, E, F and Table S2A, C, E, F). In some cases, they also tend to have larger plant sizes and stem diameters compared to the three other categories, as well as greater distances between the first and second flower (*P* < 0.05, Fig. 2G, H, K, Supporting Information Fig. S2G, H, K and Table S2G, H, K).

## Population genomic dataset

We obtained a total of 1 051 266 532 of read pairs across all the individuals (496 296–7 535 054 of raw reads per individual, mean 4 380 277, SD = 1 044 261) and 357 155 reads (< 1%) were removed because of low quality. We genotyped a total of 9 766 159 loci (composed of 2 065 743 449 sites) including 794 414 variant sites (SNPs) with Stacks. Average read depth ranged from 7× to 102.7× (mean 63.6×, SD = 16.1X) based on the combination of parameters we used. After filtering the data with populations, we finally kept 13 772 loci from which we retained 11 089 variant sites (SNPs).

## Patterns of genetic diversity and differentiation

All sites were found to have similar levels of genetic diversity (A and $A_p$). The inbreeding index ($G_{IS}$) showed that the deviation from panmixia was higher in the French sites, with the exception of site 11, Fretma, compared to Spain (Table 1). The negative $G_{IS}$ values for site 3, Tarradell are considered unreliable as they were calculated based on the basis of only four individuals of *O. insectifera*. The AMOVA highlights a higher genetic differentiation between *O. insectifera* from the control sites (7, Versols-et-Lapeyre and 12, Rodome) and *O. aymoninii* from sites 8, Le Buffre; 9, Col de Montmirat; 10, Nivoliers; and 11, Fretma with a $F_{CT}$ of 0.072 compared to *O. insectifera* from the control sites and *O. subinsectifera* from sites 1, Les Guilleries; 2, Puig Grifó; and 6, Seva, which displays a $F_{CT}$ of 0.038 (Table S3 in the Supporting Information).

PCA, with PC1 and PC2 explaining 62.49% and 20.11% of the total genetic variance respectively (Fig. 3A), shows that all individuals morphologically assigned to *O. aymoninii* cluster together (although individuals from site 11, Fretma were slightly separated from those from the other sites 8, Le Buffre; 9, Col
<figure><img src="image_4.png" /><figcaption>**Figure 3.** A, PCA displaying the two first axes (PC1 and PC2) representing 62.49% and 20.11% of the total genetic variance. PCA was computed based on 239 individuals from 12 sites genotyped at 11 089 SNPs, and colours depict sampling localities. B, PCA displaying the two first axes (PC1 and PC2) representing 62.28% and 20.56% of the total genetic variance. PCA was computed based on 196 individuals by removing three sites (i.e. 3, Tarradell; 4, Castanyola; and 5, Seva) and 8611 SNPs. *O. insectifera* control sites are represented by squares, solid ellipses distinguish the three clusters and dashed ellipses delineate pure *O. subinsectifera* from *O. insectifera* and intermediates.</figcaption></figure>

de Montmirat; and 10, Nivoliers). In addition, all individuals assigned to *O. insectifera* from France (at sites 8, Le Buffre; 9, Col de Montmirat; 10, Nivoliers; 11, Fretma) cluster together with the reference sites of 7, Versols-et-Lapeyre and 12, Rodome. However, the individuals morphologically assigned to *O. insectifera* from Spain do not cluster with *O. insectifera* from France. Instead, they were found to cluster with other individuals from Spanish sites. Indeed, in Spain, the pattern is not as clear as in France since we can observe two subclusters. The first one groups the sites 1, Les Guilleries; 2, Puig Grifó; 6, Seva, four individuals from site 5, L'Estanyol (i.e. 22-Div-025, 22-Div-026, 22-Div-027, 22-Div-044) and 10 individuals from site 4, Castanyola (i.e. 22-Div-068, 22-Div-070, 22-Div-071, 22-Div-073, 22-Div-075, 22-Div-076, 22-Div-077, 22-Div-078, 22-Div-079, 22-Div-080) corresponding to *O. subinsectifera*. The second subcluster groups the remaining individuals from sites 4, Castanyola and 5, L'Estanyol, and all the individuals from site 3, Tarradell, corresponding to intermediate individuals between *O. insectifera* and *O. subinsectifera*. Removing sites containing both *O. insectifera* and *O. subinsectifera* and/or intermediate individuals (i.e. 3, Tarradell; 4, Castanyola; and 5, L'Estanyol) and leaving sites containing only individuals with *O. subinsectifera* phenotypes, results in an even clearer clustering of the species into three entities (with PC1 and PC2 explaining 62.28% and 20.56% of the total genetic variance respectively), this time corresponding to three species *O. subinsectifera*, *O. insectifera*, and *O. aymoninii* (Fig. 3B).

Similarly, the clustered co-ancestry matrix obtained from fineRADstructure (Fig. 4) shows three main clusters, grouping together *O. insectifera* and *O. subinsectifera* and suggesting less similarity with *O. aymoninii*. A first cluster included all *O. aymoninii* individuals together. A second one grouped all *O. insectifera* from France together, regardless their site of origin.

The third cluster was found to be divided into two subclusters, one consisting of intermediate individuals from the sites 3, Tarradell; 4, Castanyola; and 5, L'Estanyol, and the other one consisting of *O. subinsectifera* from sites 1, Les Guilleries; 2, Puig Grifó; and 6, Seva. However, four individuals from 5, L'Estanyol (i.e. 22-Div-025, 22-Div-026, 22-Div-027, 22-Div-044) and one individual from 4, Castanyola (22-Div-073) are clustered with *O. subinsectifera* from sites 1, Les Guilleries; 2, Puig Grifó; and 6, Seva. Spanish individuals exhibit higher rates of relatedness between them among sites compared to French individuals. Indeed, sites 3, Tarradell; 4, Castanyola; and 5, L'Estanyol, which contain intermediate individuals, display particularly high rate of relatedness. Conversely, in France, site 11, Fretma shows relatively high relatedness among individuals.

Furthermore, the clustering analysis carried out with ADMIXTURE suggests an optimal number of genetic clusters of *K* = 8 based on the cross-entropy criterion (Supporting Information, Fig. S3). We plotted the Q-matrices (that contain individual admixture coefficients) corresponding to the best replicate runs for *K* = 3 (to investigate whether the program could distinguish the three taxa), *K* = 4 and *K* = 8. Notably, the *K* = 3 graph distinguishes well *O. aymoninii*, *O. insectifera* from France, and *O. subinsectifera* from sites 1, Les Guilleries; 2, Puig Grifó; and 6, Seva in light blue as well as four individuals from the site 5, L'Estanyol (22-Div-25, 22-Div-26, 22-Div-27, 22-Div-44). However, the individuals from site 3, Tarradell and the remains individuals from sites 4, Castanyola (22-Div-73, 22-Div-77, 22-Div-78, 22-Div79, 22-Div-80) and 5, L'Estanyol displayed intermediate admixture proportions (between 20% and 80%) between the two last clusters described. The pattern is overall similar at *K* = 4 and *K* = 8 even if they highlight particular sites like the site 11, Fretma, which seems to be apart from
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Co-ancestry matrix between each pair of individuals inferred using fineRADstructure. Each pixel depicts the magnitude of the individual co-ancestry coefficient between two individuals. Low co-ancestry coefficients values are depicted by yellow colours, whereas high values are indicated by darker colours. Different sampling localities are represented by the different colours on the side.</figcaption></figure>

*O. aymoninii* of the other sites or group of individuals that may be considered as relatives.

## DISCUSSION

The population genomic approach used in this study, for the first time on the *Ophrys insectifera* group, confidently revealed differentiation patterns that were not discernible through morphology.

### Patterns of differentiation among the representatives of the *Ophrys insectifera* clade in Southwestern Europe

In terms of morphology, our results are consistent with previous findings showing that the flowers of both endemics (*O. aymoninii* and *O. subinsectifera*) are smaller than those of the widespread *O. insectifera* (*Triponez et al. 2013*, *Gervasi et al. 2017*, *Delforge 2021*). This pattern is likely to have evolved as a result of adaptation to smaller pollinator insect species: *Andrena combinata* (for *O. aymoninii*) and *Sterictiphora gastrica* (for *O. subinsectifera*) compared to the larger *Argogorytes mystaceus* (or *A. fargei*) for *O. insectifera*. In addition to dimensions, the flowers of the three taxa also show different proportions. The flowers of *O. subinsectifera*

look like small flowers of *O. insectifera* with very short petals whereas the flowers of *O. aymoninii* have labella that are relatively wider than those of *O. insectifera*, but with similar petal length. Both endemic taxa have yellow margins on their labellum and *O. aymoninii* tend to have yellowish or greenish petals, which contrast with the generally brown petals observed in the other two species. At non-floral traits, *O. insectifera* individuals tend to be taller and have more flowers but there is non-negligible overlap between species, as well as a significant amount of intra-specific variation, presumably influenced by local environmental factors, as already suggested by (Gervasi et al. 2017).

We clearly identified three genetic clusters (Figs 3, 4 and Supporting Information Fig. S3) whereas previous studies (*Triponez et al. 2013*, *Gervasi et al. 2017*) generally failed to discriminate the three species on the basis of genetic data. We can *a priori* rule out sampling bias to explain such discrepancies as these two previous studies were carried out in the same geographic context. Instead, this suggests that the increased resolution provided by population genomic approaches (nGBS) may allow us to better highlight subtle genetic differentiation between such emerging lineages. Two of these clusters correspond to *O. aymoninii* and French individuals of *O. insectifera*.
The two taxa are genetically distinct with no evidence of recent hybridization even at sites (e.g. 8, Le Buffre; 9, Col de Montmirat; and 11, Fretma) where individuals of both species co-occur and grow at < 1 m apart. Even the three individuals that we identified as phenotypically intermediate between *O. aymoninii* and *O. insectifera* (i.e. 22-Div-146; 22-Div-152, and 22-Div-195) could be genetically assigned to *O. insectifera* for 22-Div-146 and 22-Div-195, and to *O. aymoninii* for 22-Div-152. After a reexamination of the pictures, individuals 22-Div-146 and 22-Div-152 indeed present atypical floral dimensions but which are compatible with the natural intra-specific variation observable. Also, they display coloration pattern and labellum shapes that are consistent with *O. insectifera* for 22-Div-146 (strongly trilobed labellum and absence of yellow margins) and *O. aymoninii* for 22-Div-152 (less marked trilobed labellum and presence of yellow margins). Although genetically assigned to *O. insectifera*, the individual 22-Div-195 looks clearly phenotypically intermediate: its flower's shape is typical of *O. aymoninii* but its coloration resembles the one of *O. insectifera*. One possible explanation is that this individual could be the result of ancient hybridization and subsequent backcrossing events (with *O. insectifera*). The genomic signature of the other parental species (*O. aymoninii*) may thus have been erased except at putatively adaptive genomic regions that determine floral phenotypes that may have resisted gene flow. The species barrier between *O. insectifera* and its counterpart from the Grands Causses (*O. aymoninii*) thus appears to be significant. This situation contrasts with the pattern observed on the other side of the Pyrenees. On the Spanish side, we observed a clear clustering of individuals corresponding to *O. subinsectifera*, composed of the sites 1, Les Guilleries; 2, Puig Grifó; 6, Seva; and some individuals from site 4, Tarradell and 5, L'Estanyol (i.e. 22-Div-25, 22-Div-26, 22-Div-27, 22-Div-44, 22-Div-73, which are consistent with all the genetic analyses). All other Spanish individuals form a subcluster that appears as intermediate between *O. subinsectifera* and the French population of *O. insectifera*. This latter subcluster is composed of all the remaining individuals from sites 4, Tarradell and 5, L'Estanyol, and all the individuals from site 3, Tarradell. This suggests that the barriers to gene flow between *O. insectifera* and *O. subinsectifera* are not as strongly established as between *O. insectifera* and *O. aymoninii*. The situation becomes clearer when the individuals of those later three sites are excluded (i.e. if only sites 1, Les Guilleries; 2, Puig Grifó; and 6, Seva are retained) in a way that we observe a genetic cluster formed of all *O. subinsectifera*, which is distinct from the ones formed by *O. aymoninii* and the *O. insectifera* individuals from France.

## Pervasive gene flow suspected in the Iberian Peninsula

In contrast to the observations made for *O. insectifera* and *O. aymoninii* in France, the genetic differentiation pattern suggests that at least some sites in Spain may appear as a genetic hybrid swarm, despite individuals being clearly morphologically as-signable to either *O. insectifera* or *O. subinsectifera* in most cases. The intermediate proportions of admixture shown by these in-dividuals support this hypothesis. Notably, these admixture proportions do not allow predicting the phenotype category, whether it is *O. insectifera*, *O. subinsectifera*, or an intermediate between the two. Nevertheless, this result calls into question

the existence of genetically 'pure' *O. insectifera* on the Iberian Peninsula. In this respect, the inclusion of an individual identified as *Ophrys insectifera* from Bercedo (Castilla y Léon), a region of Spain where *O. subinsectifera* should not be present, and genotyped using the same protocol (Gibert et al. 2024) also presents an intermediate status between French (at least with those from sites 7, Versols-et-Lapeyre and 12, Rodome) and Spanish individuals (Supporting Information, Fig. S4), suggesting that patterns of inter-specific gene flow may have varied from east to west of the Pyrenees, but that the existence of genetically unadmixed *O. insectifera* in Northern Spain is questionable. The possibility of inter-specific gene flow may also be supported by $G_{IS}$ values that appear to be closer to panmixia expectations in Spain.

Patterns of presumed hybridization in Spain do not appear to be arranged according to geography nor ecology. All the sampling sites studied in this study are <10 km apart from each other and do not appear to show strong environmental differentiation. In *Ophrys*, reproductive barriers are generally thought to be pre-zygotic and maintained by specific insect pollinator species. Despite this relatively high degree of plant-insect relationship, it is not necessarily one-to-one (Joffard et al. 2019), and in some cases, secondary pollinator species may be responsible for inter-specific gene flow. To date, the sawfly *Sterictiphora gastrica* has been the only pollinator species documented for *O. subinsectifera*. Notably, on 11 May 2023, we observed and captured a small solitary bee with several pollinia attached to its head, later formally identified as *Andrena afzeliella* (by D. Genoud), pseudo-copulating on an individual of *O. subinsectifera* at site 1, Les Guilleries. Whether or not this species may be responsible for inter-specific gene flow between *Ophrys* species remains to be proven, but this observation documents a new pollinator species for the endemic *O. subinsectifera* and supports the hypothesis that secondary pollinators may play a role in the evolution of *Ophrys* (see also Claessens and Kleynen 2016 mentioning *Chalicodoma sicula* as additional pollinator).

## On the systematics and taxonomy of the *O. insectifera* clade

Our results confirm that population genomic approaches such as the one used in this study (i.e. nGBS) are fruitful for discriminating *Ophrys* lineages that have recently begun to diverge, and show for the first time that this is also the case for the *Ophrys insectifera* group. Even in the sampling localities where both species occur together, we could unambiguously distinguish both endemics from the widespread *Ophrys insectifera* genetically (as we could on the basis of morphology) with the notable exception of a couple Iberian individuals being intermediate between *O. insectifera* and *O. subinsectifera*. This situation confirms that pairs of closely related *Ophrys* taxa may show variable patterns of genetic differentiation between regions, as suggested by previous work (e.g. Soliva and Widmer 2003). Indeed, neutral genetic differentiation between allopatric populations of a given *Ophrys* taxon is sometimes higher than genetic differentiation between sympatric different *Ophrys* species. One hypothesis to explain this pattern is that even a small amount of pervasive inter-specific gene flow may be locally sufficient to erase the already weak level of differentiation between incipient species.
Furthermore, although the three taxa can be to some extent separated on the basis of population genomic data, we still could not reliably infer the systematic relationships between them. The widespread and generalist *O. insectifera* could intuitively be considered as occupying a basal position in the phylogenetic species tree of the group, but our results, as well as phylogenetic reconstruction rather suggest a basal position for *O. aymoninii* (*Supporting Information, Fig. S5*). This could support an unexpected status of relictual ancestral lineage for the micro-endemic *Ophrys aymoninii*. However, secondary contact and inter-specific gene flow after the initial divergence between *O. subinsectifera* and *O. insectifera* may also make them cluster together. Aside, gene flow between *O. aymoninii* and another species outside the *insectifera* clade could have given this pattern. This needs to be tested, but in the field, we observed putative hybrids between *O. aymoninii* and a local representative of another *Ophrys* lineage, the *O. sphegodes* group because of which introgression may also alter the position of *O. aymoninii* in the species tree.

## CONCLUSION

This study confirms previous research corroborating that the endemics (*O. aymoninii* and *O. subinsectifera*) can be easily distinguished from the widespread *O. insectifera* on the basis of size, proportion, and colour of their flowers. Similar patterns are observed between both endemics suggesting parallel adaptation to smaller but distinct pollinators, with notably shorter labellum length and the presence yellow margins. At the genome level, this population genomic approach allowed us to differentiate for the first time both endemics from *O. insectifera* and highlight a putative hybrid swarm between *O. subinsectifera* and *O. insectifera*. This pattern, observed only in Spain could be due to inter-specific gene flow linked to a secondary pollinator. Although further evidence is required to confirm this, this phenomenon could bias the inference of the systematic relationships between these three species by erasing the already weak phylogenetic signal.

## SUPPLEMENTARY DATA

Supplementary data is available at *Botanical Journal of the Linnean Society* online.

## ACKNOWLEDGEMENTS

We thank P. Aymerich, R. Souche, F. Hopkins, R. Rudelle, S. Coisne, L. Barré, D. Genoud, the Parc National des Cévennes, and Takh, as well as the Bertrand and Zuliani families for the support and information provided.

## FUNDING

This work was primarily supported by an Agence Nationale pour la Recherche Jeune Chercheur Jeune Chercheuse (ANR JCJC) grant to J.B., grant number ANR-21-CE02-0022-01 and is set within the framework of the ‘Laboratoires d’Excellences (LABEX)’ TULIP [ANR-10-LABX-41].

## CONFLICT OF INTEREST

None declared.

## DATA AVAILABILITY

Sequencing data have been submitted to the European Nucleotide Archive (ENA; https://www.ebi.ac.uk/ena browser/home) under study with primary accession number PRJEB75177 (and secondary accession number ERA2969608277).

## REFERENCES

Alexander DH, Novembre J, Lange K. Fast model-based estimation of ancestry in unrelated individuals. *Genome Research* 2009;19:1655–64. https://doi.org/10.1101/gr.094052.109

Baguette M, Bertrand JAM, Stevens VM *et al*. Why are there so many bee-orchid species? Adaptive radiation by intra-specific competition for mnesic pollinators. *Biological Reviews of the Cambridge Philosophical Society* 2020;95:1630–63. https://doi.org/10.1111/brv.12633

Bateman RM, Guy JJ, Rudall PJ *et al*. Evolutionary and functional potential of ploidy increase within individual plants: somatic ploidy mapping of the complex labellum of sexually deceptive bee orchids. *Annals of Botany* 2018;112:133–50. https://doi.org/10.1093/aob/mcy048

Bateman RM, Rudall PJ. Morphological continua make poor species: genus-wide morphometric survey of the European bee orchids (*Ophrys* L.). *Biology* 2023;12:136. https://doi.org/10.3390/biology12010136

Bateman RM, Rudall PJ, Murphy ARM, *et al*. Whole plastomes are not enough: phylogenomic and morphometric exploration at multiple demographic levels of the bee orchid clade *Ophrys* sect. *Sphegodes*. *Journal of Experimental Botany* 2021;72:654–81. https://doi.org/10.1093/jxb/eraa467

Bates D, Mächler M, Bolker B *et al*. Fitting linear mixed-effects models using lme4. *Journal of Statistical Software* 2015;67:1–48. https://doi.org/10.18637/jss.v067.i01

Bertrand JAM, Gibert A, Llauro C *et al*. Whole plastid genome-based phylogenomics supports an inner placement of the *O. insectifera* group rather than a basal position in the rapidly diversifying *Ophrys* genus (Orchidaceae). *Botany Letters* 2021;168:452–7. https://doi.org/10.1080/23818107.2021.1893216

Bou Dagher-Kharrat M, Abdel-Samad N, Douaihy B *et al*. Nuclear DNA C-values for biodiversity screening: case of the Lebanese flora. *Plant Biosystems - An International Journal Dealing with all Aspects of Plant Biology* 2013;147:1228–37. https://doi.org/10.1080/11263504.2013.861530

Breitkopf H, Onstein RE, Cafasso D *et al*. Multiple shifts to different pollinators fuelled rapid diversification in sexually deceptive *Ophrys* orchids. *The New Phytologist* 2015;207:377–89. https://doi.org/10.1111/nph.13219

Catchen JM, Amores A, Hohenlohe P *et al*. Stacks: building and genotyping loci de novo from short-read sequences. *G3* 2011;1:171–82. https://doi.org/10.1534/g3.111.000240

Catchen J, Hohenlohe PA, Bassham S *et al*. Stacks: an analysis tool set for population genomics. *Molecular Ecology* 2013;22:3124–40. https://doi.org/10.1111/mec.12354

Claessens J, Kleynen J. Orchidées d'Europe: Fleur et pollinisation, Mêze, France: *Biotope Editions* 2016.

Delforge P. Orchidées de France, de Suisse et du Benelux. 3e éd. *Les Guides Delachaux* Delachaux et Niestlé, Paris, p. 352, 2021

Devey DS, Bateman RM, Fay MF *et al*. Friends or relatives? phylogenetics and species delimitation in the controversial European Orchid Genus *Ophrys*. *Annals of Botany* 2008;101:385–402. https://doi.org/10.1093/aob/mcm299

Gervasi DDL, Selosse M-A, Sauve M *et al*. Floral scent and species divergence in a pair of sexually deceptive orchids. *Ecology and Evolution* 2017;7:6023–34. https://doi.org/10.1002/ece3.3147

Gibert A, Buscail R, Baguette M *et al*. Holocene climate change promoted allopatric divergence and disjunct geographic distribution in a bee orchid species. *Journal of Biogeography* 2024. https://doi.org/10.1111/jbi.14998

Gibert A, Louty F, Buscail R *et al*. Extracting quantitative information from images taken in the wild: a case study of two vicariants of the
*Ophrys aveyronensis* species complex. *Diversity* 2022; **14**:400. https://doi.org/10.3390/d14050400

Hartig F. DHARMA: Residual Diagnostics for Hierarchical (multi-level/ mixed) Regression Models. R Package version 0.4.6 2018.

Ho J, Tumkaya T, Aryal S *et al*. Moving beyond *P* values: data analysis with estimation graphics. *Nature Methods* 2019; **16**:565–6. https://doi.org/10.1038/s41592-019-0470-3

Joffard N, Buatois B, Schatz B. Integrative taxonomy of the fly orchid group: insights from chemical ecology. *Die Naturwissenschaften* 2016; **103**:77. https://doi.org/10.1007/s00114-016-1403-y

Joffard N, Massol F, Grenié M *et al*. Effect of pollination strategy, phylogeny and distribution on pollination niches of Euro-Mediterranean orchids. *Journal of Ecology* 2019; **107**:478–90.

Jombart T, Ahmed I. adegenet 1.3-1: new tools for the analysis of genome-wide SNP data. *Bioinformatics* 2011; **27**:3070–1. https://doi.org/10.1093/bioinformatics/btr521

Léveillé-Bourret E, Chen B-H, Garon-Labrecque M-E *et al*. RAD sequencing resolves the phylogeny, taxonomy and biogeography of Trichophoreae despite a recent rapid radiation (Cyperaceae). *Molecular Phylogenetics and Evolution* 2020; **145**:106727.

Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. *Bioinformatics* 2009; **25**:1754–60. https://doi.org/10.1093/bioinformatics/btp324

Malinsky M, Trucchi E, Lawson DJ *et al*. RADpainter and fineRADstructure: population Inference from RADseq Data. *Takezaki N* (ed.). *Molecular Biology and Evolution* 2018; **35**:1284–90. https://doi.org/10.1093/molbev/msy023

Meirmans PG. genodive version 3.0: Easy-to-use software for the analysis of genetic data of diploids and polyploids. *Molecular Ecology Resources* 2020; **20**:1126–31. https://doi.org/10.1111/1755-0998.13145

Nei M. *Molecular Evolutionary Genetics*. New York, NY: Columbia University Press, 1987.

Piñeiro Fernández L, Byers KJRP, Cai J *et al*. A Phylogenomic analysis of the floral transcriptomes of sexually deceptive and rewarding European Orchids, *Ophrys* and *Gymnadenia*. *Frontiers in Plant Science* 2019; **10**:1553. https://doi.org/10.3389/fpls.2019.01553

Russo A, Alessandrini M, El Baidouri M *et al*. Genome of the early spider-orchid *Ophrys sphegodes* provides insights into sexual deception and pollinator adaptation. *Nature Communications* 2024; **15**:1–12.

Schlüter PM, Xu S, Gagliardini V *et al*. Stearoyl-acyl carrier protein desaturases are associated with floral isolation in sexually deceptive orchids. *Proceedings of the National Academy of Sciences of the United States of America* 2011; **108**:5696–701. https://doi.org/10.1073/pnas.1013313108

Sedeek KEM, Scopepe G, Staedler YM *et al*. Genic rather than genome-wide differences between sexually deceptive *Ophrys* orchids with different pollinators. *Molecular Ecology* 2014; **23**:6192–205. https://doi.org/10.1111/mec.12992

Soliva M, Widmer A. Gene flow across species boundaries in sympartic, sexually deceptive *Ophrys* (Orchidaceae) species. *Evolution* 2003; **57**:2252–61. https://doi.org/10.1111/j.0014-3820.2003.tb00237.x

Triponez Y, Arrigo N, Pellissier L *et al*. Morphological, ecological and genetic aspects associated with endemism in the Fly Orchid group. *Molecular Ecology* 2013; **22**:1431–46. https://doi.org/10.1111/mec.12169

Turco A, Albano A, Medagli P *et al*. Comparative cytogenic of the 36-chromosomes genera of Orchidinae subtribe (Orchidaceae) in the Mediterranean region: a summary and new data. *Plants* 2023; **12**:2798. https://doi.org/10.3390/plants1215279810.3390/plants12152798