Spatial genetic structure in the Madeiran endemic Dactylorhiza foliosa (Orchidaceae)

Oceanic islands have isolated biota, which typically include many endemic species. However, island endemics are vulnerable due to small population sizes, and they are often threatened by habitat destruction or by introduced pests and predators. Adequate conservation planning requires good information on genetic variability and population structure, also when seemingly viable species are considered. Here, I analysed the genetic structure in the terrestrial orchid Dactylorhiza foliosa, which is endemic to Madeira. This species is a characteristic component of evergreen laurel forests occupying the northern slopes of the island. Levels of diversity in both the plastid genome and in the nuclear genome were comparable to levels of diversity found in congeners growing in continental regions. Within populations, plants separated by distances up to 256 m shared plastid haplotypes significantly more often than plants at random, but when nuclear markers were considered, only plants growing closer than eight metres were significantly more closely related. Analysis of plastid marker variation revealed that gene dispersal by seeds is not sufficiently large to counterbalance the accumulation of mutations that build up divergence between the most distant populations. However, differentiation in the nuclear genome was considerably smaller, suggesting that gene dispersal by pollen is much more efficient than gene dispersal by seeds in D. foliosa. The overall pollen to seed dispersal ratio, mp/ms, was 7.30. Considering genetic parameters, conditions for long-term persistence of D. foliosa on Madeira seem to be good.


Introduction
Plants and animals inhabiting isolated islands have long fascinated travellers and students in natural history. Islands hold disproportionally high numbers of endemic species compared to continental regions, and the numbers of endemics per unit surface area may be extremely high (Myers et al. 2000;Kiera et al. 2009;Veron et al. 2019). Studies of island biota inspired Wallace and Darwin to formulate the theory of evolution by natural selection (Cox et al. 1976), and islands provide excellent opportunities to study various aspects of speciation.
Species numbers on oceanic islands can be elevated as result of niche differentiation. Factors promoting such differentiation include climate variability associated with altitudinal gradients, expressed as e.g. differences in precipitation, temperature or predomint wind directions. A high species diversity can also depend on high speciation rates, which can be elevated as an effect of small population sizes (Cox et al. 1976;Irl et al. 2017). Organismal groups on archipelagos of oceanic islands may additionally have expanded the number of species by the process of adaptive radiation, i e., by diversification into multiple niches when colonizing new islands, often multiplied by back-colonization of islands where the lineage is already represented by ancestral forms.
Notwithstanding their importance, island biota occupy small areas and must be considered as vulnerable in a conservation context. Plants and animals inhabiting isolated islands typically (i) comprise restricted numbers of individuals, and (ii) contain reduced levels of genetic diversity as compared to related species living in mainland regions (Hamrick and Godt 1990;Frankham 1997). A small population size leads 1 3 30 Page 2 of 12 to an increased risk of extinction due to habitat destruction or stochastic fluctuations in numbers of individuals (Lande 1999), whereas a restricted genetic diversity may impede adaptive responses to environmental change in the future (Frankel and Soulé 1981;Frankham 1997). Many island species are narrow endemics and their island habitats have often been subject to habitat destruction to a much higher degree than corresponding mainland habitats. Furthermore, island species are often threatened by predators, pests or competitors unintentionally or intentionally introduced by man. For all these reasons, many island species are threatened and of concern in conservation (Frankham 1997).
Whereas data on population sizes can be sensored by ordinary methods of observation, data on genetic diversity are difficult to obtain and are lacking for most island organisms. However, genetic data are necessary for estimating levels of population subdivision, local migration patterns, gene flow and other parameters useful for conservation planning. Genetic variation is a prerequisite for biological organisms to respond to future changes in the environment, including effects of climate warming. For island organisms, genetic diversity is especially important, since these organisms are stuck where they are and do not have the option to migrate along climatic gradients as the global climate is successively changing.
In the present study, I analysed the genetic structure in the orchid Dactylorhiza foliosa (Soland. ex Lowe) Soó, which is restricted to the island of Madeira, Portugal. The island is situated in the Atlantic Ocean about 845 km SW of the Iberian Peninsula. On the island, D. foliosa is widespread in the evergreen laurel forest, "Laurisilva" covering the northern slopes at altitudes between ca 500 and 1100 m a. s. l. (Delforge 2001;Baumann et al. 2006;Capelo et al. 2005). Here, it grows most frequently in semi-open habitats, and is common on rocky slopes, on roadsides, along paths, and along the water-channels, "levadas", which run horizontally along the mountain slopes and which lead water to cultivated areas on the drier south-facing slopes of the island. The orchid also has scattered occurrences in the slightly drier habitats above and below the laurel forest zone at altitudes between ca. 300 to 1500 m a. s. l. (Jardim and Francisco 2000;Delforge 2001;Baumann et al. 2006;Rankou 2011).
The plant is a herbaceous perennial between 20 and 80 cm tall. The leaves are more or less linear, shining, usually unspotted and aggregated along the lower portion of the stem. The inflorescence contains about 10-50 pinkish flowers with a lip of ca 8-16 × 9-18 mm covered with spots and short streaks (Jardim and Francisco 2000;Delforge 2001;Baumann et al. 2006). The plant usually forms discrete individuals with one or a few flowering stems from an underground tuber, but sometimes whole clumps of plants are formed as the result of successive division of the tuber.
Dactylorhiza foliosa is a diploid member of the D. maculata s.l. species complex (Delforge 2001). Molecular studies show that the closest relative to D. foliosa is the western European subpopulation of D. maculata subsp. maculata (Hedrén et al. 2001;Bateman et al. 2003;Devos et al. 2005;Brandrud et al. 2020). However, D. maculata subsp. maculata is invariably tetraploid, and the distinct morphology of D. foliosa suggests that it has been separated from the mainland populations of D. maculata for a long period of time. Other members of the D. maculata complex are known from the mountain regions of NW Africa (Delforge 2001), but these taxa are poorly studied genetically. The complex is lacking from other Macaronesian islands.
Available data on D. foliosa are inconclusive for estimating its evolutionary potential and adaptability. On the one hand, the orchid is restricted to the ca 150 km 2 of laurel forests on Madeira, and it is not known from anywhere else, which suggests that the genetic diversity in the species might be restricted. On the other hand, the orchid is widespread within the laurel forest habitat, including secondary vegetation (Baumann et al. 2006), and the population size is large (Rankou 2011), which suggests that the species could store a considerable amount of genetic diversity. Moreover, laurel forests covers much of the northern slopes of Madeira without any major gaps and the geographic subpopulations are probably well connected (Fig. 1), which should allow for genetic connectivity and good opportunities for the maintenance of the genetic diversity that exists within the species.
Similar to other species of Dactylorhiza, D. foliosa is characterised by efficient mechanisms for seed and pollen dispersal. Like most orchids, D. foliosa has minute, partly air-filled seeds that can be transported long distances by wind. The flowers contain no reward but act by deceit. The syrphid fly Eristalix tenax and the bumble-bee Bombus maderensis have been observed as the most frequent pollinators (Fernandes et al. 2005).
A better understanding of genetic diversity in island endemics is interesting in itself and may also be important for implementation of conservation strategies. In the present study, I have chosen to analyse variation in both plastid markers and nuclear markers in D. foliosa. By combining data from the resulting data sets, it is possible to compare the effects of gene dispersal by seeds and gene dispersal by pollen. By sampling populations from different localities on Madeira and by mapping physical distance between individuals within sampling localities, gene dispersal and genetic diversity structure could be described for both the species level and the population level. Finally, I compare the results obtained for D. foliosa with similar data obtained from previously analysed species of Dactylorhiza, and ask whether the genetic structure in D. foliosa would be of concern in a conservation context.

Material and sampling
Samples of D. foliosa were taken from eight sites, separated by between 0.7 and 26 km (Online resource 1). Plants were sampled along levadas, paths or roadsides. Three to forty plants were sampled from each site, depending on the density and extent of the local population. To avoid repeated sampling of the same genet, samples were taken from plants separated by at least one meter. The average distance between adjacent samples within sites ranged between 12.8 and 151.2 m (Online resource 1). The minimum distance between samples within sites were 1 or 2 m, whereas the maximum distance ranged between 71 and 1051 m. Some populations were much larger than the sampling areas and might be connected without any major gaps in distribution to adjacent populations, for instance the populations at Portela and Ribeiro Frio.
A few flowers or a small amount of leaf tissue was collected for DNA extraction, leaving the rest of the plants for continued growth. The sampled material was dried in silica gel (Chase and Hills 1991) and stored at room temperature until extraction. DNA was extracted from one dried flower with bract or from ca 1 cm 2 dried leaf tissue according to the 2 × CTAB procedure (Doyle and Doyle 1990). From one individual per population an additional flower was sampled and deposited in the Lund botanical museum (LD) as voucher material.

Molecular analyses
To enable discrimination between plastid genomes, ten repeat regions described in Hedrén et al. (2008) were initially screened for variation; 1, 4, 5, 6, 6B, 8, 9, 10b, 11b, 17, 18, and 19B (Online resource 2). After initial screening, the four regions 6B, 10b, 11b, and 19B were selected for analysis of the entire material. Locus 6B is a complex indel region, 10b is a polyA-polyTA-polyT repeat region and 11b and 19B are mononucleotide repeat regions (Online resource 2). Each combination of fragment sizes at the four repeat regions was recognized as a separate haplotype. Finally, a haplotype distance matrix for downstream analyses was constructed by counting the number of repeat regions that differed in fragment size between any two haplotypes.
The ITS region was studied as in Pillon et al. (2007), in which ITS alleles were identified by their combined fragment lengths at two size-variable regions (Online resource 2). ITS has previously been shown to vary in the D. maculata s.l. complex, in some cases even within taxa (Devos et al. 2005;Pillon et al. 2007;Ståhlberg and Hedrén 2010).
Variation in the nuclear genome was studied by analysis of four microsatellite loci (Online resource 2). Two of these, ms3 and ms11 gave rise to fragments separated by steps of three base-pairs, whereas ms2 and D2501 gave rise to fragments separated by steps of two base-pairs. The first two loci were described by Nordström and Hedrén (2007), and the latter two were added in .
Size variable fragments in the plastid genome were amplified by an initial round of denaturing at 94 °C for 2 min, followed by 40 cycles of 94 °C for 1 min, 50-58 °C (Online resource 2) for 1 min, 72 °C for 1 min 30 s, ended by a final extension of 72 °C for10 min. Nuclear microsatellite loci were amplified by an initial denaturing at 94 °C for 3 min, followed by 36 cycles of 94 °C for 30 s, annealing at a temperature 54 °C for 30 s, extension at 72 °C for 45 s, and Queimados concluded by a final extension period at 72 °C for 10 min. The ITS alleles were amplified by means of the same program, but with an annealing temperature of 52 °C. The PCR product from each reaction was mixed with 8.5 μl formamide containing appropriate size markers to enable exact relative size determination of the amplified fragments. The fragments were then heated at 94 °C for ca 20 min before loading on a horizontal acrylamide gel and separated by size on an ALFexpress II automated sequencer (Amersham Pharmacia Biotech). Fragment sizes were determined by using ALFwin™ fragment analyser 1.03.01 software (Amersham Pharmacia Biotech).

Data analysis
For each population sample, the following statistics were calculated on basis of nuclear microsatellite data (Table 1): A, average number of alleles per locus; A e , effective number of alleles per locus; AR, allelic richness adjusted for a population size of 5 samples (10 gene copies); H e , expected diversity corrected for sample size; and F IS , inbreeding coefficients calculated according to Weir and Cockerham (1984).
For plastid haplotype data, the following parameters were calculated for each population sample: A, number of different haplotypes; A e , effective number of haplotypes; AR, haplotype richness adjusted for a population size of 5 samples (5 haplotype copies); and H e , expected diversity corrected for sample size (Table 1). All statistics were obtained in SPAGeDi 1.4 (Hardy and Vekemans 2002).
The general differentiation between population samples was described by F ST s calculated according to an ANOVA approach as described Weir and Cockerham (1984). Separate analyses were performed for the plastid and nuclear data sets. To test for an isolation by distance structure (Wright 1943), the differentiation between pairs of populations were calculated and compared to geographic distances by means of Mantel tests (Table 2). Differentiation was given as pairwise F ST /(1-F ST ) and geographic distances as the logarithm of geographic distance (Rousset 1997). I also tested for an association between pairwise F ST derived from nuclear data with pairwise G ST from the haplotype data. Pairwise genetic distances were calculated in SPAGeDi 1.4 (Hardy and Vekemans 2002). Associations between distance matrices were tested by10000 rounds of permutation in NTSYS-pc 2.2 (Rohlf 2005).
To analyse D. foliosa for patterns of genetic differentiation, I constructed a maximum parsimony network for the plastid haplotypes using the computer program Network 10.2.0.0 (Fluxus technology 2007), and I investigated the microsatellite data for potential clustering of nuclear genotypes using the computer program Structure 2.3.4 (Pritchard et al. 2000). For the plastid loci, I assumed a stepwise mutation model and weighted all mutations equally, irrespective of the amount of variability shown at the different marker sites. For the nuclear data, I used the mixed ancestry model, assumed independent allele frequencies, and made no a priori assumptions of population belonging. I tested the data for separation into K = 1 to 10 clusters and repeated the analysis 20 times for each value of K. I used a burnin period of 10,000 steps and made 10,000 iterations in each run. The optimum number of clusters was estimated by the Evanno method (Evanno et al. 2005) as implemented in the computer program Structure Harvester (Earl and vonHoldt 2012), and the optimal combined solution for the multiple runs of each value for K was estimated by the computer program CLUMPP (Jakobsson and Rosenberg 2007).
To test for a phylogeographic structure, i.e. whether different haplotypes or microsatellite alleles are significantly more divergent between populations than within populations (Hardy and Vekemans 2002), I compared differentiation measures that describe the degree of differentiation between markers with values obtained by permutation of the original data. For the plastid haplotype data the differentiation measure N ST was calculated (Pons and Petit 1996). Here, genetic divergence was estimated as the number of loci with different fragment sizes between any pair of haplotypes. A corresponding G ST , along with 95% confidence intervals, was obtained by 10,000 rounds of permutation of genetic distances in the haplotype distance matrix. For the nuclear microsatellites. I estimated the statistic R ST , which is based on the numbers of mutational steps between microsatellite alleles (Slatkin 1995;Michalakis and Excoffier 1996). A corresponding F ST , along with 95% confidence intervals, was estimated by 10,000 rounds of permutation of allele sizes among alleles within loci (see also Hardy and Vekemans 2002). Calculations were performed in SPAGeDi 1.4 (Hardy and Vekemans 2002).
The F ST and the F IS estimates from nuclear microsatellites and the F ST from plastid haplotypes were used to calculate the pollen to seed-flow ratio, mp/ms according to the formula given by Ennos (1994) and Petit et al. (2005): where F ST(B) and F IS(B) are derived from biparentally inherited microsatellite markers and F ST(M) is derived from maternal plastid markers.
The spatial genetic structure, SGS, within populations was estimated by use of the kinship coefficient F ij of Loiselle et al. (1995). I calculated mean kinship between individuals within distance intervals of 4, 8, 16, 32, 64, 128 and 256 m, which subdivided the material into categories containing appropriate numbers of pairwise comparisons, and analysed whether kinship decreased with distance. To test whether the kinship values were significant, I compared the observed values with the mean values obtained by 10,000 rounds of permutation of individual locations among all individuals. These results were also used to calculate the Sp statistic of Vekemans and Hardy (2004): bF is the linear regression slope of F ij (n) on the natural logarithm of the distance intervals, and F ij (1) is the mean kinship coefficient for the first distance interval.
Separate analyses were performed for plastid haplotype and nuclear microsatellite data. All within-population statistics were calculated in SPAGeDi 1.4 (Hardy and Vekemans 2002).

Results
A total of 33 plastid haplotypes was found in D. foliosa (Online resource 3). The small population of Encumeada West-Central comprising only three samples, was found with a single haplotype (Table 1, Online resource 4). In the remaining populations, between four (Encumeada West) and 18 (Encumeada North) haplotypes were encountered within populations. Correcting for sample size by rarefaction with k = 10 (excluding Encumeada West-Central and Rabaçal), between 3.82 (Encumeada West) and 7.53 (Encumeada North) haplotypes per populations was estimated (Table 1). Expected haplotype diversity (H e ) ranged from 0.0 (Encumeada West) and 0.69 (Ribeiro Frio) to 0.93 (Encumeada N).
In ITS all plants were fixed for the ITS allele designated as ITS-I in Pillon et al. (2007).
A total of 29 different alleles were encountered at the four examined nuclear microsatellite loci (Online resource 5). Locus ms13 was fixed for an allele of 90 bp in all populations. At the remaining loci, three alleles were found at ms 11 with between one and two alleles per population, six alleles at ms 14 with between two and four alleles per population, and 19 alleles at locus D2501 with between three and fifteen alleles per population (Online resource 5). Adjusting for sample size by rarefaction and comparing all populations with the population with the lowest number of individuals (Encumeada VC; n = 3, six gene copies), the number of alleles per population ranged from 1.00 to 1.53 at ms11, 1.25-2.97 at ms14, and 3.00 to 5.10 at D2501 (Online resource 5). The expected gene diversity (H e ) ranged from 0.35 (Encumeada West) and 0.391 (Pico das Pedras) to 0.575 (Rabaçal) ( Table 1, Online resource 5).
The global inbreeding coefficient over populations and loci was moderate in D. foliosa, with F IS = 0.092 (Table 3). Inbreeding coefficients varied over populations, ranging from values close to zero in the population samples from Ribeiro Frio and Portela to F IS = 0.8 and F IS = 0.523, respectively, in the small samples from Encumeada West-Central and Rabaçal (Table 1). The average unweighted F IS over populations samples was 0.222 (Table 1). Disregarding the smallest population samples, the average inbreeding coefficient in the remaining six samples was 0.075.
The overall F ST for plastid haplotypes was 0.137 (Table 3). A clearly lower amount of differentiation between populations was observed in the nuclear genome and the overall F ST estimated on basis of nuclear microsatellites was 0.018.
Comparison of pairwise F ST values derived from plastid haplotype data and geographic distances between sampling sites by Mantel tests (Table 2) did not reveal presence of any significant isolation by distance structure, IBD (Z = 0.716, p (random Z ≥ obs Z) = 0.204). For the nuclear genome, comparison of pairwise F ST values derived from microsatellite allele frequency data and geographic distances between sampling sites did also not reveal any significant isolation by distance structure (Z = -0.173, p (random Z ≥ obs Z) = 0.758). A somewhat stronger correlation was found between pairwise F ST s derived from plastid haplotypes and pairwise F ST s derived from nuclear microsatellites, but the association was still not significant (Z = 0.952, p (random Z ≥ obs Z) = 0.102).
The relationships of the encountered plastid haplotypes were investigated in a parsimony network analysis. Eight most parsimonious trees of 40 steps each were found. The resulting network is given as Fig. 2. To some extent, the more common haplotypes were often shared between neighbouring sites, or else neighbouring sites often shared related haplotypes with each other.
To test D. foliosa for presence of any substructure on basis of nuclear genomes, nuclear microsatellite genotypes were analysed for aggregation into one to ten (K) groups by means of Bayesian cluster analysis as implemented in the computer program Structure. However, whatever the number of K, each sample invariably obtained about equal probability of belonging to any of the groups tested for (Online resources 6). Applying the Evanno method to obtain the optimal number of subgroups, K = 4 had a slightly higher probability than K = 2 (Online resources 7), but none of these solutions were visibly any better than any of the other solutions.
Calculation of N ST and comparison with the corresponding G ST obtained by permutation of numbers of differences between pairs of haplotypes in the haplotype distance matrix (Table 3), revealed the presence of a significant general phylogeographic structure in the plastid genome of D. foliosa (0.400 vs 0.229; p < 0.000 *** ).
In contrast, comparison of the R ST , which considers the number of mutational steps between microsatellite alleles at locus, and the corresponding F ST , obtained by permutation of allele sizes among alleles within each locus, revealed no significant phylogeographic structure in the nuclear genome (R ST = 0.021, F ST = 0.014; p = 0.308 ns ).
Based on the global F ST values derived from plastid haplotype data (F ST (M) = 0.137), the global F ST s derived from nuclear microsatellite data (F ST (B) = 0.018), and the inbreeding coefficient derived from nuclear microsatellite data (F IS(B) = 0.092; Table 3), the pollen to seed flow ratio in D. foliosa was calculated to mp/ms = 7.30.  (Pons and Petit 1996) 0.400 Permuted value (G ST ) (permutation of rows and columns of distance matrix between haplotypes) 0.229 p = 0.000 *** 0.021 Global R ST (Michalakis and Excoffier 1996) Permuted value (F ST ) (permutation of allele sizes among alleles within loci) 0.014 p = 0.308 ns Spatial autocorrelation calculations within geographic distances up to 256 m (Fig. 3) revealed a much steeper decline in pairwise kinskip, F ij , in the plastid genome than in the nuclear genome. For the plastid genome, kinship values were on average 0.528 in pairwise comparisons up to 4 m and were still significantly different from the permuted values at differences up to at least 256 m. For the nuclear genome, kinship coefficients were on average 0.133 in pairwise comparisons up to 4 m and were not significantly different from permuted values at distances above eight meters.  The Sp statistic, describing the decline in average kinship values with distance, was estimated to Sp = 0.0264 for plastid haplotype data, and Sp = 0.0026 for the nuclear microsatellite data (Online resource 8).

Affinity to European members of Dactylorhiza
The molecular data obtained in the present study agree with a position of D. foliosa close to D. maculata subsp. maculata. All samples examined had a plastid haplotype corresponding to that of haplotype group II of Devos et al. (2003). Such haplotypes dominate in western European D. maculata subsp. maculata, but are more uncommon in D. maculata subsp. fuchsii (Ståhlberg and Hedrén 2010). In ITS, all plants of D. foliosa were fixed for the ITS allele designated as ITS-I in Pillon et al. (2007). The same ITS allele is strongly dominant in western European D. maculata subsp. maculata, but is rare in subsp. fuchsii (Pillon et al. 2007;Ståhlberg and Hedrén 2010). Of the four nuclear loci analysed, ms13 was fixed for an allele of 90 bp in all populations. This allele is also strongly dominant in western European D. maculata subsp. maculata and virtually absent from subsp. fuchsii (Hedrén and Ståhlberg, unpublished). Still, the diploid D. foliosa cannot be derived directly from the tetraploid D. maculata subsp. maculata, but rather from a common group of diploid ancestors to both of these taxa (Brandrud et al. 2020). Dactylorhiza foliosa also differs markedly from the European D. maculata subsp. maculata in general morphology, including the tall stature, the elongate and usually unspotted leaves and in the large pinkish flowers with indistinct markings. It is possible that these characteristics have evolved in relatively short time as a result of directional selection after the arrival of the ancestral population to Madeira. However, D. foliosa is also reminiscent of poorly studied members of the D. maculata group distributed in mountain regions of northwest Africa, including D. maculata subsp. battandieri and D. maculata subsp. maurusia (Landwehr 1977;Delforge 2001;Baumann et al. 2006). Further studies are obviously needed to understand the exact position of D. foliosa within the D. maculata s.l. complex.

Patterns of between-population differentiation in D. foliosa
A compilation of studies mostly based on allozyme data (Frankham 1997) has shown that island populations or island taxa have in most cases lower levels of genetic diversity than related continental counterparts. However, individual studies vary considerably in different directions, and some plants may have equally high amounts of gene diversity in island populations as in continental populations (e.g., García-Verdugo et al. 2015).
The levels of gene diversity observed in D. foliosa were comparable to levels found in other species of Dactylorhiza. The population average for nuclear markers, 0.453, was within the range observed in mainland taxa (0.355-0.668; Hedrén et al. , 2021. The gene diversity obtained for plastid haplotypes, 0.734, was even somewhat higher than that obtained in other taxa (0.259-0.675; op. cit.). Although estimated levels of genetic diversity depend on the degree of variability of the markers used in any particular study, the studies cited here for comparison with D. foliosa were based on partly the same markers, so the comparison should still be relevant.
The global F ST calculated on basis of nuclear markers was 0.018 in D. foliosa. This value is lower than in most other orchids previously analysed for genetic diversity structure. The low F ST value was also in agreement with the lack genetic structuring found by Bayesian clustering in Structure. Average estimates of G ST /F ST in orchids range from 0.087 (Hamrick and Godt 1996), over 0.146 (Phillips et al. 2012), to 0.187 (Forrest et al. 2004). Phillips et al. (2012) compared levels of F ST for rare and common terrestrial orchids and found that the average F ST was significantly lower in the latter category (F ST = 0.279 vs F ST = 0.092, respectively). Dactylorhiza foliosa is characterised as a species that is common within its range (Rankou 2011). Accordingly, the very low F ST value and the absence of structuring in the nuclear genome could be an effect of the relatively high abundance and even distribution of D. foliosa in the laurel forest habitat on Madeira. The laurel forest is more or less continuous on the northern side of the island (Fig. 1), and there appears to be no major barriers to dispersal within the distribution range.
It should also be observed that the total area occupied by D. foliosa is small. The sampled populations were all relatively close to each other and are probably more efficiently connected by gene flow than populations of other studied orchid species, which have larger and patchier distributions.
A third factor likely to contribute to a low F ST in D. foliosa is the pollination system. Like other members of Dactylorhiza, D. foliosa is characterized by a deceptive pollination system (Baumann et al. 2006) and it offers no nectar to visiting insects. Although rewarding orchids have higher fruit set than orchids offering no reward to their pollinators (Gill 1989, Neiland andWilcock 1995), it has been suggested that the pollinators of non-rewarding orchids ensure more efficient cross-pollination as they transport pollen longer distances between plants (Nilsson 1992, Neiland and Wilcock 1995, Jersáková et al., 2006. Empirical data show that orchids characterized by deceptive pollination systems have on average a significantly lower level of population differentiation than orchids offering nectar or other rewards to their pollinators (F ST = 0.25 vs F ST = ca 0.13; Cozzolino and Widmer 2005).
Many orchid species suffer from inbreeding depression following selfing (Smithson 2006). Only a small proportion of orchids are autogamous or asexual. Inbreeding levels are often moderate in orchids, but fixation indices, F IS , may range from 0 to 1. Inbreeding levels in other members of Dactylorhiza range from F IS = 0.15-0.66 (Hedrén & Lorenz 2019). The estimated level of inbreeding was moderate in D. foliosa with mean F IS = 0.222 (Table 1), and somewhat lower still if the two smallest population samples were excluded. As the two smallest populations did hardly differ from other sampled populations in plant density (online resource 1), it is possible that the small number of potential mates, rather than the mean distance between plants, provides the most likely explanation to the observed variation in inbreeding levels.
The populations of D. foliosa were more strongly differentiated in plastid haplotypes than in nuclear microsatellites (F ST = 0.137). However, this value is also quite low in comparison with other orchid species, including members of Dactylorhiza in which it ranges from 0.258 to 0.707 (Hedrén et al. 2021).
Since the plastid genome is exclusively dispersed by seeds in orchids (Corriveau and Coleman 1988), whereas the nuclear genome is dispersed by both seeds and pollen, F ST estimates based on plastid data are expected to be consistently higher than F ST s based on nuclear markers (Petit et al. 2005). Still, orchids have minute seeds that are supposedly easily dispersed long distances by wind (Arditti and Ghani 2000), and the seeds might contribute disproportionally much to gene dispersal, as compared to other plant families with heavier seeds.
The two F ST measures can be combined with the inbreeding coefficient, F IS , to calculate the pollen to seed flow ratio of Ennos (1994), also given as mp/ms in Petit et al. (2005). Whereas orchids often have low values for this ratio, sometimes close to one (Squirrel et al. 2001;Cozzolino et al. 2003), Petit et. al. (2005 calculated a median value of mp/ ms across all investigated seed plants as mp/ms = 17. The pollen to seed flow ratio in D. foliosa was calculated to mp/ms = 7.32. This ratio is slightly higher than in estimates obtained from other members of Dactylorhiza investigated so far (Hedrén et al. 2021). However, the pollen to seed flow ratio might be dependent on geographic scale (McCauley 1997;Hamrick and Trapnell 2011;Hedrén et al. 2021). Whereas pollen dispersal typically declines gradually from the pollen source, depending on the flight patterns of the pollinators, seed dispersal distances typically exhibit a strongly leptokurtic distribution, with large amounts of seeds deposited in the close vicinity of the source plant and only a minor portion dispersed at medium and long distances (Cain et al. 1998). In effect, the relative contribution of pollen dispersal to total gene dispersal is small at short distances, higher within populations and between close populations, but again small as compared to seed dispersal at long distances (McCauley 1997;Hamrick and Trapnell 2011). For example, in the hummingbird-pollinated epiphytic orchid Laelia rubescens in Costa Rica, mp/ms ratios up to 13.67 were recorded for populations separated by less than 2 km, but a ratio of 1.40 between populations at a scale up to 10 km (Trapnell and Hamrick 2004), Likewise, in the terrestrial, bumble bee-pollinated Dactylorhiza umbrosa investigated on the Anatolian-Iranian highland plateau, mp/ms ratios of 8.4 -12.01 were recorded at distances up to 122 km, but ratios of 2.8-4.9 at higher distances (Hedrén et al. 2021). It is possible therefore, that the high value recorded for D. foliosa, with a maximum distance between populations of 26 km, reflects the relatively small sampling area.
Alternatively, the comparatively high pollen to seed flow ratio recorded for D. foliosa could be a consequence of the sheltered forest habitat in which the species is usually found. A forest acts as shelter to strong winds, especially at ground level, and might thus reduce dispersal distances of small, wind-borne seeds. In contrast, as pollen is dispersed by actively flying insects, pollen dispersal in a forest habitat might not be as strongly affected. Other species of Dactylorhiza that have been examined for genetic diversity structure are typically found in wet grasslands or in open fen habitats where seeds are more efficiently dispersed by winds.
There was no significant correlation between pairwise F ST s derived from plastid data and geographic distance between populations, and likewise, there was no significant correlation between pairwise F ST s derived from nuclear data and geographic distance (Table 2). Although the amount of population differentiation as measured by F ST is generally low in orchids, Phillips et al. (2012) found that many common orchids are still characterized by an isolation by distance structure; IBD. However, they also found that IBD was most often detected in studies where the scale of sampling exceeded 250 km, far longer than the 26 km separating the most distant populations sampled of D. foliosa.
I also tested each of the plastid and nuclear microsatellite data sets for the presence of a phylogeographic structure, such that divergent haplotypes or microsatellite alleles in different population are in general more strongly differentiated than divergent haplotypes or alleles present in the same population (Hardy and Vekemans 2002). In D. foliosa there was no evidence for any significant phylogeographic structure in the nuclear microsatellite data, but for the plastid data, N ST , which takes numbers of differences between haplotypes into account, was significantly higher than the permuted value corresponding to an ordinary G ST (0.400 vs. 0.229, p < 0.001 *** ). An indication of a phylogeographic structure in the plastid genome is also given by the parsimony network ( Fig. 3) in which it could be observed that geographically adjacent sampling sites often share the same or related plastid haplotypes with each other; viz. the three subsites at Encumeada, the close sites Quiemadas and Pico das Pedras, and the two easternmost sampling sites Ribeiro Frio and Portela. However, there is also cases where distant populations share the same haplotypes. Such cases might be explained by rare long-distance dispersal events, or simply homoplasy in the data. The interpretation of the combined patterns is that gene dispersal by seeds is not efficient enough to balance out the differentiation that slowly accumulates as result of local mutations within populations of D. foliosa. In contrast, pollen dispersal seems to be effective to counterbalance differentiation by mutation in nuclear microsatellites.

Spatial genetic structure
The genetic differentiation patterns within sampling sites was analysed by spatial autocorrelation and by comparison of mean values of kinship, F ij , for different distance intervals at distances up to 256 m (Fig. 3, Online resource 6). Whereas plants close to each other often shared the same plastid haplotypes and still did so significantly more often at distances up to at least 256 m, the kinship estimate derived from nuclear microsatellites was small already at close distance and kinship values were not significantly different from permuted values at distances above eight metres. Accordingly, pollen dispersal contributes much more to total gene dispersal than seed dispersal also at the local scale.
The Sp value calculated for the nuclear data (Sp = 0.0026) was very low in comparison with most other orchids examined (range 0.0015-0.069 in predominantly outbreeding species; Chung et al. 2011;Hedrén and Lorenz 2019). However, the value found in D. foliosa was still comparable to some populations of the Asian terrestrial species Cymbidium goeringii (Chung and Nason 2007).
The Sp value calculated on basis of plastid haplotype data was 0.0264 in D. foliosa. There are no other studies in which Sp statistics for the plastid genome has been reported for an orchid species. However, in some orchids characterized by uniparental reproduction, and in which pollen dispersal does not contribute to gene dispersal, Sp values have been found to be higher than in the plastid genome in D. foliosa, with values ranging from 0.05 to 0.208 (Hedrén and Lorenz 2019).