Molecular evidence of species- and subspecies-level distinctions in the rare Orchis patens s.l. and implications for conservation

Characterizing genetic diversity and structure of populations is essential for the effective conservation of threatened species. Orchis patens sensu lato is a narrowly distributed tetraploid species with a disjunct distribution (i.e., Northern Italy, North Africa and the Canary Islands), which is facing a severe decline. In this study, we evaluated levels of genetic diversity and population structuring using 12 new nuclear microsatellite markers. Our analyses of genetic differentiation based on multiple approaches (Structure analysis, PCA analysis, and F-statistics using the ploidy-independent Rho-index) showed that gene flow is low across the range of O. patens s.l., particularly in the Canary Islands. Clear differences in allele frequencies between Italy, Algeria and the Canary Islands underlie the genetic differentiation retrieved. Our study provides support for the recognition of O. canariensis as a sister species to O. patens and the separation of the Italian populations as a new subspecies of O. patens. Despite the high heterozygosity values found in all populations (ranging from 0.4 to 0.7), compatible with the tetraploid status of the species, small population sizes and reduced gene flow will be likely detrimental for the different populations in the long term, and we recommend immediate conservation actions to counteract further fragmentation and population decline.


Introduction
As we are facing the sixth age of extinction Barnosky et al. 2011), an increasing number of species are threatened by decreased population sizes, habitat degradation and habitat fragmentation. Endangered species are in most cases characterised by small, fragmented and declining populations (Brooks et al. 2002). Such limited groups of individuals are at high risk of extinction owing to decreased genetic diversity and the growing effect of biological stochasticity (Fahrig 2003;Henle et al. 2004). There is strong evidence, indeed, that genetic factors, such as low genetic diversity and consequently reduced reproductive fitness, are driving threatened taxa to extinction (Spielman et al. 2004;Vellend and Geber 2005;Garner et al. 2020). Preserving genetic diversity is therefore one of the essential requirements to ensure population survival Laikre et al. 2020), especially under anthropogenic and environmental pressures, including climate change (Hewitt 2004;Jump et al. 2009). Rare species with low genetic diversity are more vulnerable to environmental and biological stresses compared to species with higher genetic diversity (Leimu et al. 2006;Honnay and Jacquemyn 2007;Aguilar et al. 2008), and small, isolated and declining populations of rare species, under climate change, may also become maladapted (Borrell et al. 2020). Habitat fragmentation or the complete loss of the habitat are major reasons for depletion of genetic variation in plants (Young et al. 1996;Frankham et al. 2002). Generally, species with narrow distribution ranges or small and/or isolated populations are more affected by loss of their natural habitat, and they are more susceptible to extinction due to depletion of genetic diversity than those species that are widely distributed and have relatively large, non-fragmented populations (Hamrick and Godt 1996;Frankham et al. 2002;Pandey et al. 2015), although life-history traits may play an important role in buffering genetic diversity loss (Honnay and Jacquemyn 2007).
Spatial distribution of genetic variation within a species may offer useful insights into fundamental ecological and environmental processes and information important for identifying distinct genetic groups (if present) across its range (Allendorf and Luikart 2007;Mkare et al. 2017;Gargiulo et al. 2019a). The delineation of conservation units for an endangered species is critical not only for their long-term survival, but also for directing the prioritization of conservation efforts (Petit et al. 1998;Schwartz et al. 2007;Volkmann et al. 2014;Mkare et al. 2017;Médail and Baumel 2018). Characterizing genetic diversity and structure of populations is useful in informing assisted gene flow (Borrell et al. 2020), or to identify high priorities for effective conservation of threatened species (Fay 2018;Väli et al. 2019). This is even more urgent for orchids, as approximately half of the 1641 currently assessed in the IUCN Red List are in the categories Vulnerable, Endangered and Critically Endangered (Fay 2020;IUCN, 2020).
The genus Orchis, which includes about 20-25 species, is divided into two subgenera, Orchis and Masculae, with the first commonly called the "anthropomorphic group" (Kretzschmar et al. 2007). The topology of Bateman et al. (2003) shows the non-anthropomorphic Orchis species are well-supported as monophyletic. Cytological examinations of different Orchis taxa have revealed a main chromosome number of 2n = 42 (e.g. Mrkvicka 1992;Aedo and Herrero 2005), but exceptions with aberrant values have been found in section Robustocalcare (2n = 40; Hautzinger 1978), which includes the only tetraploid species in the genus, O. patens. However, more recent counts indicated tetraploid O. patens with 2n = 84 (Pridgeon et al. 2001;Bernardos et al. 2006).
Orchis patens sensu lato is a species with a disjunct geographic range, with O. patens subsp. patens occurring in Italy (Liguria region) and North Africa (Algeria and Tunisia), 1 3 and O. patens subsp. canariensis occurring in the Canary Islands (Orsenigo et al. 2016 Bateman et al. 2003;Bernardos et al. 2006), unlike WCSP (2020, in which O. canariensis is considered as the basionym of O. patens subsp. canariensis. Although not resolving relationships in O. patens s.l., ITS-based phylogenetic analyses showed that O. patens s.l. is the sister group of O. spitzelii (Bateman et al. 2003), a species predominantly found in Mediterranean regions, including northern Africa and western Asia, with outlying populations on the island of Gotland. A recent study on seed micromorphology of O. patens subsp. patens, however, identified some similarities with the O. mascula group (Calevo et al. 2017), a group of taxa distributed in Europe, north-western Africa and western Asia, thus opening up the hypothesis of an allotetraploid origin. Orchis patens is assessed in the Italian (Orsenigo et al. 2016) and European IUCN Red Lists (Rankou 2011) as Endangered (EN), and in the Mediterranean IUCN Red List (Calevo et al. 2018) as Vulnerable (VU), but all the assessments agree that, given the clear morphological differences between subspecies and the lack of studies investigating genetic diversity, clarifying taxonomic uncertainties at species and subspecies levels should be a priority for the conservation of O. patens s.l.
Moreover, although O. patens is reported as the only tetraploid species of the genus (Mrkvicka 1992;Pridgeon et al. 2001), it is not clear whether different cytotypes (ploidy levels) occur (Bernardos et al. 2006) and codominant markers, especially microsatellites, may provide indirect information through the allele counts observed (Besnard and Baali-Cherif 2009;Gompert and Mock 2017). In this study, we designed new nuclear microsatellite markers to elucidate genetic structure and diversity in O. patens s.l., which will aid in the development of appropriate conservation strategies. In particular, our aims were to i) ascertain the occurrence of gene flow among populations, with a focus on indirectly detecting differences in ploidy; ii) evaluate the infraspecific delimitation and taxonomic rank of O. patens subsp. patens and O. patens subsp. canariensis; iii) verify the potential use and informativity of the markers herein developed in other congeneric species.
We refer to Canarian populations as O. canariensis and to African and Italian populations as O. patens s.s. from here on.

Plant material
We collected leaf samples from populations of both subspecies, representative of the distribution range of the species. Samples of O. patens s.s. were collected in north-western Italy (in figures referred to as OP samples, Fig. 1a and d), in four localities of the Liguria region (16 samples from Romaggi, 11 samples from Breccanecca, 12 samples from Portofino and six samples from Capreno), and in two localities of Algeria ( Fig. 1b and e): ten samples from Ait Zikki, in the Bejaia region (in figures referred to as OPAZ samples), and ten samples from Ain Talhi, in the Souk-Ahras region (in figures referred to as OPSA samples), carefully avoiding natural hybrids. Samples of O. canariensis ( Fig. 1c and f) were collected from Tenerife (ten samples, in figures referred to as TN), Gran Canaria (ten samples, in figures to referred as GC) and La Gomera (nine samples, in figures referred to as LG). Collection of plant material was performed after requesting permission from the competent authorities and respecting regional rules.

DNA extraction and Illumina sequencing
Total genomic DNA was extracted from silica-dried leaf tissues using the DNeasy Plant Mini Kit (QIAGEN) according to the manufacturer's instructions, after disruption with TissueLyser. The quantity and quality of DNA samples were assessed by spectrophotometry (ND-1000 Spectrophotometer NanoDropH; Thermo Scientific, Wilmington, Germany) and by electrophoresis on a 1% agarose gel.
A genome-skimming approach (Straub et al. 2012;Dodsworth 2015) was used for the detection of microsatellites (Viruel et al. 2018;Landoni et al. 2020). Genomic libraries were prepared using NEBNext® UltraTM II DNA Library Prep Kit for Illumina® (New England Biolabs, Ipswich, MA, United States) with AMPure XP magnetic beads for size selection (300-350 bp) and NEBNext® Multiplex Oligos for Illumina® (Dual Index Primer Sets I and II) as barcodes for simultaneous sequencing ). Library quality was evaluated using a QuantusTM fluorometer (Promega Corp.) and an Agilent 4200 TapeStation (Agilent Technologies, Santa Clara, CA, United States). Multiplexed  (Bolger et al. 2014) was used to remove low-quality reads and adapter sequences flagged by FastQC v0.11.7 (Andrews 2010).

Microsatellite design and analysis
Microsatellites were identified and primers designed with MSATCOMMANDER 1.0.8 (Faircloth 2008); duplicate reads and those containing more than one SSR array were manually checked and removed. Primers fulfilling the following criteria were selected: optimal size 18-22 bp, not directly flanking the SSR motif, lacking ambiguous bases, low selfand pair-product complementary parameters, amplicon expected size < 300 bp, and melting temperature (Tm) difference < 1.5 °C (Viruel et al. 2018).
Preliminary screening was carried out for 40 nuclear primer pairs on a 2% agarose gel. After testing and optimization of primers with a subset of samples, the final set of markers consisted of 12 microsatellite primer pairs (Table 1);other primer pairs were discarded due to issues with reproducibility and PCR product yield, possibly resulting from the tetraploid status of the species. Representative sequences for the microsatellite sequences were deposited in GenBank (accession numbers MT799099-MT799110).
PCR mix was prepared by adding 6 μL of (2 ×) DreamTaq PCR MasterMix (Ther-moFisher Scientific), 0.5 μL of 0.4% (w/v) bovine serum albumin (BSA), 0.5 μL of each primer (10 μM), 1 μL of template DNA (ca. 10 ng/μL), and H 2 O up to a final volume of 10 μL. Amplifications were performed with the following PCR conditions: 4 min of initial denaturation at 94 °C followed by 25 to 35 cycles (see Table 1 for primer details) of 30 secs of denaturation at 94 °C, 30 secs of annealing (see Table 1 for temperature details), 45 secs of elongation at 72 °C and an extra extension of 10 min at 72 °C, with the exception for primer pairs Opat_27 for which it was necessary to perform a touch-down PCR by decreasing annealing temperature of 0.2 °C at each cycle (see Table 1). One μL of PCR product was added to 10 μL HiDi™ formamide (Applied Biosystems) and 0.15 μL 500 ROX Size Standard (Applied Biosystems) and capillary electrophoresis was run on an ABI3730 DNA Analyzer (Applied Biosystems).
Fragment sizing was performed in GeneMapper v.5 (Applied Biosystems) and peaks were visually checked before assigning final allele sizes. We inferred the minimum chromosome number from the maximum number of peaks detected; this procedure has been shown to give results consistent to flow cytometry (Besnard and Baali-Cherif 2009;Stark et al. 2011;Nemorin et al. 2013;Teixeira et al. 2014;Donkpegan et al. 2015;Gompert and Mock 2017). We used the R package polysat v1.7 (Clark and Jasieniuk 2011; Clark and Schreier 2017) to explore segregation patterns in our microsatellite loci. Understanding whether loci are polysomic or are inherited as separate subgenomes can give some clues about the origin of polyploidy in O. patens and, most importantly, can inform about the correct treatment of the dataset for subsequent analyses (Dufresne et al. 2014;Hardy 2016;. To avoid the bias induced by population structure, we kept individuals from Italy, Algeria and the Canary Islands as three separate populations before performing the analyses in polysat; we used the function processData-setAllo with default parameters and excluding the locus Opat_38 (see results). To

Genetic diversity and structure
As we did not have prior information about the extent of gene flow in O. patens s.l., we conducted the analyses both by grouping individuals based on geographical regions (i.e., Italy, Algeria and the Canary Islands) and on sampling localities (i.e., nine sampling localities across the three regions). Standard population genetic parameters such as mean number of alleles (Num), effective number of alleles (Eff_num), observed heterozygosity (H O ), gene diversity (H S ) and inbreeding coefficient G IS based on the 12 nuclear microsatellites were estimated using the program GenoDive v.3 (Meirmans and Van Tienderen 2004;Meirmans 2020). In addition, we performed a hierarchical analysis of molecular variance (AMOVA; Meirmans 2012) based on the Rho-statistics (ρ-statistics), which is ploidy-independent (Ronfort et al. 1998). All analyses in GenoDive were performed by applying the correction for unknown dosage of the alleles. Test for deviations from Hardy-Weinberg proportions was not attempted since it was not possible to record allelic dosage information for all the loci Gargiulo et al. 2019b). We used the programme SPAGeDI v.1.5 (Hardy and Vekemans 2002) to produce a matrix of pairwise Rho values testing significance with 20,000 permutations; in addition, we estimated the selfing rate per locality using the standardized identity disequilibrium (Hardy 2016). Monomorphic loci were excluded from the computations.
We further analyzed genetic structure using different approaches. In polysat, we constructed a neighbor-joining tree based on a square dissimilarity matrix and we performed a principal component analysis (PCA) based on Bruvo distances (Bruvo et al. 2004). The model-based clustering procedure implemented in Structure v2.3.4 (Pritchard et al., 2000) was used to evaluate the number of different genetic clusters, K. We applied the Recessive Alleles model (Falush et al. 2007) to handle missing allelic dosage information. Preliminary analyses were conducted both with the Admixture model and the Non-Admixture model, testing 1 < K < 8. The final analysis was conducted with K ranging from 1 to 8, with a 5 × 10 5 burn-in period, 5 × 10 5 MCMC runs, and ten iterations per each K value, using the Non-Admixture model for the total dataset, and the Admixture model for subsets based on sampling locations and infraspecific delimitations. The most likely value of K was evaluated using ΔK (Evanno et al. 2005) and the estimators introduced by Puechmaille (2016) for uneven datasets (MedMedK, MedMeanK, MaxMedK and MaxMeanK) implemented in STRU CTU RESELECTOR (Li and Liu 2018). Structure results were summarized with CLUMPAK (Kopelman et al. 2015). Simulations studies have shown that Structure is particularly robust compared to other clustering methods when dosage information is incomplete, even in the presence of mixed ploidies (Stift et al. 2019).

Microsatellite design and analysis
Selected primer pairs, microsatellite motifs, and size ranges are reported in Table 1. The 12 newly designed nuclear microsatellite primer pairs amplified all samples and the resulting loci were variable with some exceptions. In particular, locus Opat_38 was not variable for the Italian O. patens but showed at least two alleles (Table 2) in the

3
Algerian populations, and it did not show any amplification product in the samples from the Canary Islands. Locus Opat_39, with multiple alleles in O. patens from Italy and Algeria (Table 2), was not variable for samples from the Canary Islands. The number of alleles as inferred by the maximum number of peaks in the electropherograms was four in all populations. Every genotyped individual had at least one locus displaying three peaks, except GC8 in the Canary Islands, which exhibited two peaks at every locus (Table 2). Therefore, we did not find strong evidence of populations with mixed ploidy. The analysis in polysat showed that polysomic inheritance prevailed at all loci, with one exception represented by locus Opat_36 in the Algerian samples (Fig. 2). Consequently, data were treated as for autopolyploid species in the subsequent analyses.
To verify the transferability of these microsatellite loci, we tested the primers in other Orchis taxa (Table 3)

Genetic diversity and structure
Genetic diversity parameters for each population and locality are summarized in Tables 2 and 4. The results indicate that genetic diversity indices are similar across all populations and localities, except Tenerife, which shows slightly lower heterozygosities and number of alleles (Table 4). The number of effective alleles (Eff_Num) per population varied between 3.2 (in Italy, averaged over 11 polymorphic loci) and 4.4 (in the Canary Islands, averaged over ten polymorphic loci). The lowest Ho was found in the Canary Islands populations (0.49, reflecting the low H O observed in Tenerife; Table 4) and the highest in the Italian populations (0.69); H S ranged between 0.64 in Italy and 0.72 in the Canary Islands (Table 2). Negative values of G IS were found in all populations, possibly reflecting the tetraploid status of the species; however, locus-specific high G IS were found for Opat_20 and Opat_29, possibly indicating null alleles or dropout (i.e., alleles that do not produce amplification by PCR). The highest selfing rate, as  computed in SPAGeDI was found in the Algerian locality OPAZ; absence of selfing was estimated in Italian and Canarian localities (Table 4). The hierarchical AMOVA (Table 5) showed most of the variation at the individual level, which is expected when using polymorphic (and polyploid) microsatellite loci. Moreover, differentiation among groups was the highest for the group formed respectively by Algerian/Italian populations and the Canary Islands populations (Rho_ct: 0.421) and the lowest when one of the groups included Italian and Canarian populations (Rho_ct: 0.162). In parallel, the matrix of pairwise Rho-values indicated strong differentiation among populations in Italy, Algeria and the Canary Islands and also among the individuals collected in different islands of the archipelago (Fig. 3). When performing AMOVA on the Italian subpopulations, for which we had more localities and individuals (i.e., 45), the results (Table 6) pointed to higher levels of variation within (89.8%) rather than between the subpopulations (10.2%).
The neighbor-joining tree based on the square matrix of genotype dissimilarities obtained with polysat indicated three clear clusters (Fig. 4). In particular, the cluster including all the samples from the Canary Islands (light blue) is divergent from the other two. The same pattern was observed in the principal component analysis (PCA) based on Bruvo distances (Fig. 5).
DeltaK values computed on the output of the Structure analysis indicated a strong signal for K = 2 (Supplementary information Fig. S1), whereas the other estimators (MedMedK, MedMeanK, MaxMedK and MaxMeanK) indicated three clusters (Supplementary information Fig. S2). The clusters underlined, once again, a strong genetic differentiation between O. patens s.s. and O. canariensis) (Fig. 6).
The analysis conducted on the Italian and Algerian populations identified four potential clusters (Supplementary information Fig. S3), whereas the analysis conducted on the samples from the Canary Islands identified three clusters (Supplementary information Fig. S6). In both cases, it was possible to observe an important division driven by regionality (Fig. 7) which, in the Canary Islands, are consistent with the three islands where samples were collected. For each of the Algerian and the Italian subsets, two clusters were identified (Supplementary information Figs. S4-S5) (Fig. 7).

Discussion
Preserving genetic variation is imperative to ensure the existence of any species, especially under anthropogenic and environmental stresses, including climate change (Hewitt 2004;Jump et al. 2009). For this reason, rare species with a narrow distribution or a fragmented habitat, and therefore potentially low genetic diversity, are a priority for studies aimed at investigating genetic diversity to establish appropriate conservation strategies. Orchidaceae    species ). Clear species delimitations are important for orchids because many of the currently recognized species are of conservation concern and threats continue to increase due to habitat degradation (Cozzolino et al. 2003;Fay 2018). In this study we used 12 newly developed nDNA markers to determine the level of genetic differentiation among populations of the O. patens polyploid group, with a view to resolving taxonomic uncertainties for conservation purposes (Rankou 2011;Swarts et al. 2014;Orsenigo et al. 2016;Calevo et al. 2018;Médail and Baumel 2018).

Genetic structure and taxonomic delimitation
Our analyses of genetic differentiation based on different approaches revealed the same clear clustering of samples (Figs. 3,4,5,6), suggesting that the microsatellite data presented here carry a clear signal of population subdivision. In particular, samples from the Canary Islands are included in a genetic cluster divergent from the Italian and Algerian populations and are, in turn, characterized by genetic subdivision depending on the island of provenance (Figs. 3, 7; Table 5). The failure to amplify a product using Opat_38 in the Canary Island populations further implies genetic divergence between O. canariensis and O. patens s.s. This locus was successfully amplified in the samples from Italy and Algeria. Even though both the Italian and Canarian plants were recently found to share the same mycobiont showing a high level of phylogenetic conservativism in their symbiotic interactions (Calevo et al. 2020), our results support the separation of O. canariensis as a distinct (probably sister) species, as initially described from Tenerife by Lindley (Orchis canariensis, Lindley 1835) based on morphological characters and as recognised by some authors (e.g., Aceto et al. 1999;Bateman et al. 2003) based on phylogenetic analyses. Analysis of the genetic structure based on microsatellites has proved to be informative in revealing cryptic infraspecific structure and subspecies delimitation, for example in the genera Greenwayodendron (Lissambou et al. 2019) and Lophira (Ewédjè et al. 2020). The genetic structure of O. canariensis, as characterized by the number of subpopulations and the frequencies of different genetic variants (alleles) in each subpopulation (Chakraborty 1993), might be indicative of a high degree of genetic isolation among the subpopulations on the different islands. The strong clustering of Canarian populations (that mirrors the island of provenance) has been observed in other endemic plants of the islands such as Lotus sessilifolius (Yang et al. 2018) and Bethencourtia spp. (Rodríguez-Rodríguez et al. 2018), and the genetic differentiation among island populations may reflect the very beginning of speciation as proposed for the endemic Kleinia neriifolia by Sun and Vargas-Mendoza (2017). This, as suggested by the latter authors, might be expected as the ocean acts as a natural barrier among the islands.
At the same time, there is also some evidence of differentiation between the Italian and Algerian populations of O. patens s.s. (Figs. 3, 4, 5, 6, 7), which provides a genetic basis for the treatment of the Italian populations as a subspecies of Orchis patens. As mentioned for O. canariensis, microsatellite amplifications revealed a subpopulation-specific pattern (Tables 2 and 6) even among samples from Italy and Algeria, which may suggest that mutations in the priming sites are consistent with intraspecific divergence. This may be promoted by the tetraploid status of O. patens, in which the two chromosome sets undergo a range of somatic mutations and diverge if gene flow is reduced (Moody et al. 1993;. Morphologically, the Italian populations are characterized by slightly spotted leaves (Fig. 1d), whereas spots are absent in African individuals (Fig. 1e) and were not mentioned in the first description of the species from Algeria (Desfontaines 1799).
Further morphological studies should be focused on detecting morphometric differences among the subspecies based on their entire phenotype, in order to strengthen their taxonomic circumscription and to collect more data about functional differences towards their conservation.

Influence of polyploidy on current genetic diversity patterns
Although populations of O. patens have experienced natural and anthropogenic fragmentation, and fruit set is low and affected by hybridization with the sympatric species O. provincialis, the genetic structure found here reflects historical gene flow, especially between Italy and Algeria. Orchis patens may have been relatively common in the past within its narrow range and therefore the genetic effects of recent fragmentation and reduction in population number may yet to be detected. In polyploid species, loss of genetic diversity by drift is slowed down because of the availability of multiple allele copies compared to diploid species (Moody et al. 1993;. In addition, the effect of gene flow is amplified as more gene copies are involved. Consequently, the effects of fragmentation and genetic erosion may be not evident until more generations have elapsed. The levels of heterozygosity we found in O. patens are compatible with that of other polyploid species (Gargiulo et al. 2019b;Hagl et al. 2020), and in particular tetraploid orchid populations, such as Gymnadenia conopsea (H O = 0.77; Gustafsson and Thorén 2001) and the Dactylorhiza majalis complex (H O = 0.77 averaged across tetraploids; Balao et al. 2016). It is also possible that the current distribution of O. patens reflects long-distance dispersal from one of the three areas towards the other two, although our results do not allow us to speculate further on the most likely explanation.
Regarding the origin of O. patens s.l., the group has been considered sister of the O. spitzelii group based on an analysis of the ITS region (Bateman et al. 2003). These similarities have led to the hypothesis that O. patens could have been derived from the latter species as an autotetraploid, but a recent study based on seed micromorphology as a taxonomic tool also highlighted similarities with the O. mascula group (Calevo et al. 2017), suggesting the possibility of an allopolyploid origin. Our microsatellite analysis did not recover a segregation pattern typical of allopolyploid species (allopolyploidy is generally associated with disomic inheritance, at least for some loci). However, current segregation patterns are not always informative of past events, as the inheritance may shift from disomic to polysomic and vice-versa (Stift et al. 2008). However, more analyses are needed to better understand the origin of O. patens and the potential taxa involved in the allopolyploidization (Ramsey and Schemske 1998;Soltis et al. 2010;Barker et al. 2015).

Conclusions and implications for conservation
Combining molecular evidence from our newly developed nuclear microsatellite markers with the geographical distribution of the studied populations, we propose the recognition of O. patens subsp. canariensis as a separate (possibly sister) species, i.e., O. canariensis, and the treatment of the Italian populations of O. patens as a distinct subspecies. Moreover, the primer set for the 12 microsatellite loci with transferability potential represents a useful tool for interspecific research in the genus Orchis. Our study of species delimitation within the O. patens s.l. complex can offer guidelines for further morphological and genetic investigations and conservation studies. Therefore, we recommend the following actions: • Orchis canariensis should be assessed for the IUCN red list; • the Mediterranean assessment of O. patens (Calevo et al. 2018) should be revised taking into account the taxonomic delimitation suggested here (which will probably lead to a higher category of threat), the evidence of restricted gene flow, narrow range and low number of individuals (Orsenigo et al. 2016); • given that the availability of the main symbiotic fungus of O. patens s.l. at global scale seems not to be a limiting factor for their distribution (Calevo et al. 2020), conservation actions aimed at preventing further habitat fragmentation should be taken. In particular, local authorities and agencies in Algeria, Tunisia (where the natural population is already potentially extinct), Italy and the Canary Islands should institute measures to protect remaining individuals from damage and illicit collection. • Further non-destructive sampling, even if difficult due to the rarity and low density of the species, may be required to better investigate the level of inbreeding and its contribution to the reproduction of the populations. Once reproduction mechanisms are clarified, strategies which facilitate pollination should be taken into account, also considering the low fruit-set of the species.