Collecting Sites and DNA Amplification

Samples of L. speciosa were collected in 15 forest fragments in the Cuitzeo basin, Michoacán (Table 1, Figure 1). We collected 12 to 20 individual plants from each population; the plants sampled were growing on Quercus deserticola trees separated from each other by at least 30 m. Samples were frozen until DNA extraction. Genomic DNA was extracted from 100 mg of fresh leaf material using the protocol proposed by Lefort and Douglas (1999). Eight nuclear microsatellite loci were selected and amplified by means of multiplex polymerase chain reaction (PCR; Cortés-Palomec, McCauley, & Oyama, 2008). Based on allele size, annealing temperature, and fluorescent labels, the PCR amplification products were sorted into three groups of primers. The first group included the primer pairs for Lspe8 and Lspe12; the second group those for Lspe1, Lspe4, and Lspe7; and the primer pairs for Lsp6 and Lspe10 were included in the third group (Cortés-Palomec et al., 2008). The PCR was performed using the QIAGEN Multiplex PCR kit (QIAGEN) in a 5-µl volume containing 1X Multiplex PCR Master Mix, 2 µM of each primer, dH₂O, and 20 ng of template DNA (Cortés-Palomec et al., 2008). The thermal cycling conditions consisted of 35 cycles of 94℃ for 1 min, annealing for 1 min, extension at 72℃ for 2 min, and final extension at 72℃ for 10 min. Multiplex PCR products were combined with a GeneScan-500 LIZ size standard and the analyses performed on an ABI-PRISM 3100 Avant sequencer (Applied Biosystems). Fragments were analyzed and recorded using the Peak Scanner program 1.0 (Applied Biosystems).

Table 1.

Locality Name, Sample Size, Geographical Coordinates (Degrees), Mean Number of Effective Alleles per Locus (N_e), Mean Observed Heterozygosity (H_O), Mean Expected Heterozygosity (H_E), and the Inbreeding Index (F_IS) for 15 Populations of L. speciosa From Cuitzeo Basin, Michoacán.

Genetic diversity

We tested for the presence of null alleles, large-allele dropout, and errors due to stuttering in the microsatellite data using the MICRO-CHECKER 2.2.3 program (Van Oosterhout, Weetman, & Hutchinson, 2006) with 10² bootstrap simulations and a 95% confidence interval (CI). Deviations from the Hardy-Weinberg equilibrium, such as excess heterozygosity, or F_IS deficit, were tested by means of a Markov-chain approach (with 10³ dememorization steps, 10² batches, 10³ iterations per batch) using the GENEPOP 4.1 software (Raymond & Rousset, 1995). For each of the 15 populations of L. speciosa, we estimated the following genetic diversity parameters: number of alleles per locus (A), observed heterozygosity (H_O), and expected heterozygosity (H_E), using the GENETIX 4.02 software (Belkhir, Borsa, Chikhi, Raufaste, & Bonhomme, 2004).

Genetic Structure and Bayesian Admixture Analysis

To have a clear scenario about the genetic assignment, we run two Bayesian approaches, Geneland and Structure, to test if the genetic assignments are the same or differ in the number of genetic cluster detected. The first approach, Structure, was tested to infer the genetic ancestry of each individual plant (Falush, Stephens, & Pritchard, 2003; Hubisz, Falush, Stephens, & Pritchard, 2009; Pritchard, Stephens, & Donnelly, 2000). STRUCTURE 2.3.4 software (Falush et al., 2003; Hubisz et al., 2009; Pritchard et al., 2000) uses a Bayesian clustering model to determine the proportion of each individuals' ancestry originating from different populations (Evanno, Regnaut, & Goudet, 2005). The optimal number (K) of groups was determined by varying K from 1 to 10 and running the analysis 10 times for each K value to find the maximum posterior likelihood (LnP (D)). Each run was performed using 50⁴ burn-in periods and 10⁶ Markov Chain Monte Carlo (MCMC) repetitions after burn-in. We used an admixture model that allows correlated allele frequencies without any prior information. Following Evanno et al. (2005), we determined the most likely value of K based on the maximum value of ΔK as implemented in the Structure Harvester 0.6.9 software (Earl & von Holdt, 2012). In the second approach, we run the spatial cluster model implemented in the GENELAND package (Guillot, Mortier, & Estoup, 2005) of the R 3.2.3 program (R Development Core Team, 2016). Different sets of parameters (MCMC, thinning and burn-in) were used in different test runs, in order to find the optimal parameters. Following the recommendation of the user's manual, the first step was replicated 10 times, allowing K to vary from 1to 10 clusters and using the Markov chain Monte Carlo (MCMC) repetitions were set at 100,000, thinning was set at 100, and the burn-in period was set at 200, to know the effective number of genetic clusters in the data set. Then, we set ten long runs MCMC iterations were performed, with K fixed to the maximum value of the logarithm of posterior probability of the data (PPD), using 5 × 10⁶ MCMC iterations, and the other parameters unchanged. We calculated the mean logarithm of posterior probability of the data (PPD) for each of the 10 runs and selected one with the highest PPD (Guillot et al., 2005; Guillot, Santos, & Estoup, 2008). To ensure that the run was long enough, we obtained 10 different runs and compared the parameter estimates (K, individual population membership, maps). The pairwise population genetic differentiation R_ST was estimated using the infinite allele model (IAM), performing 10⁴ permutations in the ARLEQUIN 3.5.1.2 software (Excoffier & Lischer, 2010). A hierarchical test of population structure was conducted using the stepwise mutation model (SMM) with analysis of molecular variance (AMOVA) in ARLEQUIN 3.5 (Excoffier & Lischer, 2010). We conducted a hierarchical analysis to estimate genetic variance components between groups, between populations within groups, and within populations, using the assignation of groups obtained by Structure and Geneland. Statistical significance was tested using 10⁴ permutations. The genetic exchange between L. speciosa populations was assessed by estimating the Bayesian-scaled long-term effective population size (Ne) and migration rate (m) with the MIGRATE 3.5.1 software (Beerli & Felsenstein, 2001). For all the analyses, the starting value of the chain was set to 20⁶ visited and 16 recorded genealogies, following a burn-in period of 50³ iterations.

To identify geographic and genetic breaks between L. speciosa populations, we used the Monmonier's maximum difference algorithm with the BARRIER version 2.2 software (Manni, Guerard, & Heyer, 2004). BARRIER creates a map of the sampling locations from their geographical coordinates. Barriers are then represented on the map by identifying the maximum values within the population-pairwise genetic distance matrix. We used a pairwise matrix of average square distances (Goldstein, Ruiz-Linares, Cavalli-Sforza, & Feldman, 1995; Slatkin, 1995) estimated for the 15 populations of L. speciosa. Resampling random subsets of individuals within populations with the MSA program (Dieringer & Schotterer, 2003) provided 100 bootstrap replicate distances that were used to achieve statistical significance for the predicted barriers.

Changes in population size

We used the BOTTLENECK 1.2 software (Piry, Luikart, & Cornuet, 1999) to detect population bottlenecks, using the three genetic groups previously identified with Structure. Recent bottlenecks can be identified as a population where rare alleles are the first to be lost, then the mean number of alleles per locus is reduced accordingly. In contrast, heterozygosity is less affected, producing a transient excess in heterozygosity relative to that expected based on the resulting number of alleles (Cornuet & Luikart, 1996; Luikart & Cornuet, 1998). For this test, we used 90% stepwise and 10% multistep mutations and performed 10⁴ iterations of the Wilcoxon's signed-rank test with the SMM, the IAM, and the two-phase mutation (TPM) model.

Additionally, we estimated the effective size of 15 populations of L. speciosa using the LDNe software (Waples & Do, 2008). This program implements the bias-correction method developed by Waples (2006) to obtain Ne from a sample of S individuals. We set Pcrit = 0.02 (i.e., alleles with a frequency <0.02 are excluded), which generally provides a good balance between accuracy and bias (Waples & Do, 2008). CIs for Ne were calculated with the chi-square approximation implemented in LDNe (Waples & Do, 2008).

Genetic Diversity

No evidence of null alleles was found in any of the sample-loci combinations. The tests for errors due to stuttering and large-allele dropout yielded negative results in all cases. Values of the genetic diversity parameters estimated in 15 populations of L. speciosa ranged from moderated to high in Zacazo, Indapa, and Caurio (Ne = 6.87–9.62, H_O = 0.635–0.714); Laguna, Laguni, Potrer2, Felipe, Potrer1, and Huiram (Ne = 6.87–7.62, H_O = 0.602–0.682); and in Olvido, Robles, Aguila, Coapas, Coapas2, and Correo (Ne = 6.25–7.0, H_O = 0.534–0.589), see Table 1. Wright's inbreeding coefficient within populations (F_IS) showed positive values in all cases, denoting a heterozygote deficiency (Table 1); F_IS values were low in Zacazo, Indapa, and Caurio (F_IS = 0.009–0.279); slightly high in Laguna, Laguni, Potrer2, Felipe, Potrer 1, and Huiram (F_IS = 0.124–0.215); and high in Olvido, Robles, Aguila, Coapas, Coapas2, and Correo (F_IS = 0.154–0.285). The gene exchange levels detected were moderate to high in all the populations examined (Table 2).

Table 2.

Levels of Genetic Exchange Estimated With the Program MIGRATE.

Genetic Structure and Bayesian Admixture Analysis

We applied two complementary Bayesian clustering algorithms, namely Structure (Falush et al., 2003; Pritchard et al., 2000) and Geneland (Guillot et al., 2005), to infer population structure (i.e., a number of clusters, K) and to assign individuals probabilistically to populations based on individual multilocus. Both approaches shown similar results over 10 replicated runs tested varying K from 1 to 10 to get the highest probability. In terms of the number of genetic clusters, both approaches, according to the maximum posterior likelihood (LnP (D)) and the maximum ΔK value, showed that K = 3 is the optimum number of genetic groups, see Figure 2. Also, Structure and Geneland approaches confirm in the majority of the assignment of each individual to each genetic cluster (Figures 3, 4, and 5). For instance, in the bar plot and the pie charts in Figures 3, 4, and 5 show the distribution of ancestry proportions in each collection site. Cluster 1 (e.g., red in Structure and cluster 1 Geneland) is widespread across the Cuitzeo basin and includes the Zacazo, Indapa, Laguna, Potrer2, Olvido, and Caurio populations, except for Laguna that shift for Laguni in Geneland. Cluster 2 (e.g., green in Structure and cluster 3 Geneland) includes the populations of Robles, Felipe, Aguila, and Potrer1, which are consistently structured across the landscape. Cluster 3 (e.g., blue in Structure and cluster 2 in Geneland) has a wider distribution and includes the populations of Correo, Coapas, Coapas2, Laguni, and Huiram (Figures 3, 4, and 5) except for Laguni that was changed for Laguna in Geneland. Geneland and Structure also differs in the presence of admixture (not noticeable by Geneland), as an evidence of gene exchange, particularly in the populations of Indapa, Caurio, Laguna, Potrer1, and Potrer2, which have ancestry coefficient (Q) values ranging between 0.732 and 0.848 (Figures 3 and 4).

Figure 2.

Estimated number of populations from Structure and Geneland analyses. Structure values of K plotted against K, the peak indicates the most probable number of genetic groups given the data using Structure Harvester program. Geneland plot of the number of populations simulated from the posterior distribution obtained with GENELAND.

Figure 3.

Genetic assignment of individuals and populations according to the Bayesian method implemented in the program STRUCTURE. Each thin horizontal line represents an individual and the proportion of each color is the proportion of ancestry derived from each of the three main genetic groups (K 3) inferred. Green, red, and blue genotypes are representing the genetic ancestry groups corresponding to L. speciosa populations. Populations are separated by black lines.

Figure 4.

Each pie chart represents the proportions in each population of the three genetic groups as assigned by the program STRUCTURE. Green, red, and blue genotypes are representing the genetic ancestry groups corresponding to L. speciosa populations. Genetic discontinuities (bold lines B-1 and B2) obtained with Monmonier's maximum difference algorithm on genetic distances derived from microsatellite allele frequencies.

Figure 5.

Maps of Geneland individual assignments to clusters for K = 3 (scale units in latitude, longitude). The three plots represent the assignment of pixels to each cluster: cluster 1, cluster 2, and cluster 3. The highest membership values are in light yellow and the level curves illustrate the spatial changes in assignment values. The plot is based on the highest probability run at that value of K.

The AMOVA was performed using both genetic assignation from Structure and Geneland. Results showed that most of the variation (95.7 and 95.6%) occurred within L. speciosa populations; the variation between groups accounted for 3.29% and 3.51%, and the variation between populations within groups accounted for 1.02% and 0.87% (Table 3). Similarly, there was a moderate pairwise genetic differentiation between some populations, as indicated by the R_ST values observed (Table 4): Aguila/Olvido (R_ST = 0.149), Coapas/Aguila (R_ST = 0.143), Caurio/Coapas2 (R_ST = 0.082), Laguna/Olvido (R_ST = 0.107), and Caurio/Correo (R_ST = 0.082). Finally, a low differentiation was found between the following population pairs: Indapa/Caurio (R_ST = 0.006), Zacazo/Coapas (R_ST = 0.014), and Laguni/Felipe (R_ST = 0.028), see Table 4.

Table 3.

Analysis of Molecular Variance (AMOVA) Performed on the nSSR Data and Using R_ST for the Three Group Genetic Clusters Obtained by Means of STRUCTURE and by GENELAND for Populations of L. speciosa.

Table 4.

Pairwise Population Genetic Differentiations Using the Matrix of Slatkin Linearized R_ST (SMM) Among Localities of L. speciosa in Cuitzeo Basin, Michoacán, Over 10⁴ Replicates in the Arlequin 3.5 Software (Schneider et al., 2005).

The analysis performed using 100 bootstrap replicates of the average square distance matrices revealed two barriers (with over 50% bootstrap support) between the 15 populations of L. speciosa studied (Figure 4). The most significant barrier, with 95% bootstrap support, is a complex break separating the southwestern part from the central part of the Cuitzeo basin. This indicates that, despite the admixture revealed by the Structure analysis, some populations have become isolated such as the Potrer1, Potrer2, Olvido, Coapas, Coapas 2, Laguni, Laguna, and Aguila (Figure 4). The second barrier, with 89% bootstrap support, sets apart the Caurio population (located in the northwest part of the basin) from the rest of the L. speciosa populations.

Changes in population size

Results from the analyses to test for evidence of recent bottlenecks (excess heterozygosity) using the IAM, the TPM, and the SMM are shown in Table 5. No significant results were found with the IAM and TPM. However, the SMM showed significant (p = .001, .003, .001) evidence of mild bottlenecks in populations of the three genetic groups, particularly in loci Lspe 4, 6, and 14, and also significant evidence of excess heterozygosity in the L. speciosa populations suggesting that they have recently gone through a bottleneck. The estimates of effective population size (N_e) obtained using LDNe for the L. speciosa populations showed that the third genetic group had the highest value (N_e = 441 individuals), followed by the second (N_e = 349) and the first (N_e = 122) groups. All N_e estimates had a high Jackknife support and a good CI (Table 5).

Table 5.

Bottleneck Analysis for L. speciosa Populations in the Cuitzeo Basin Using Wilcoxon Rank Test Under Infinite Allele, Stepwise Mutation, and Two-Phase Model.

Genetic Diversity

Habitat fragmentation, coupled with the massive extraction of individuals from natural populations, has reduced the genetic diversity and population size of these plants, with the ensuing reduction in population connectivity (Honnay et al., 2005; Jump & Peñuelas, 2006). In this study, we examined the genetic structure of populations of L. speciosa, an orchid species that is threatened by the indiscriminate extraction of specimens for trading in local markets. Surprisingly, the L. speciosa populations examined showed moderate to high levels of genetic diversity, as determined from nuclear microsatellites. Recent studies based on microsatellite variation have also found high genetic diversity (H_O: 0.0–0.728) in the rare orchid Isotria medeoloides (Stone et al., 2012), the critically endangered spider orchid Caladenia huegelii (H_O: 0.587–0.766; Swarts et al., 2009), and in two Ophrys species (H_O: 0.76–0.91 to 0.77–0.92; Mant, Peakall, & Schiestl, 2005; Soliva & Widmer, 2003). Most of these species share similar characteristics with L. speciosa, such as being epiphytic, long-lived, cross-pollinated, and long-distance seed dispersal. Furthermore, widespread species such as L. speciosa tend to exhibit a higher genetic diversity compared with threatened or endemic species (Li & Jin, 2007; Mathiasen, Rovere, & Premoll, 2007; Qian, Wang, & Tian, 2013; Yu, Yang, Sun, & Liu, 2011). This may be the result of historical patterns of genetic variation in L. speciosa, for example, from a formerly wide-ranging distribution to becoming rare only recently. These patterns suggest that the massive extraction of individuals of this species has not yet had a noticeable effect on gene frequencies and heterozygosity or in the alleles that cause genetic drift in populations (Aguilar, Quesada, Ashworth, Herrerías-Diego, & Lobo, 2008; Schaberg et al., 2008). Similarly, Wright's inbreeding coefficients within populations (F_IS) were positive in all cases, denoting heterozygote deficiency (see Table 1). Fragmentation may also lead to inbreeding due to either increased self-fertilization or mating between related individuals within remnant fragments (Chung et al., 2014; Honnay et al., 2005; Mathiasen et al., 2007; Swarts & Dixon, 2009). In this case, inbreeding seems to have been due to the activity of pollinators, as they were reported visiting nearby forest fragments more often than distant ones (Bacles, Lowe, & Ennos, 2004; Smith-Ramírez & Armesto, 2003). The ecological features of L. speciosa could also contribute to an increased gene flow via pollinators or long-distance dispersal between populations and suggest that isolated trees might serve as stepping-stones for gene flow between populations (Barrett & Kohn, 1991; Herrera-Arroyo et al., 2013; Soliva & Widmer, 2003; Tremblay & Ackerman, 2001).

Some of the populations of L. speciosa exhibited evidence of having recently gone through a mild-bottleneck, but only with the SMM model. Based on the well-documented history of L. speciosa in Michoacán (Miranda, 1997), we hypothesized that these populations have been subject to massive extraction and habitat fragmentation for at least 100 years. Studies on conservation genetics have put forward a rule to define the effective population size necessary to prevent genetic damage (Barrett & Kohn, 1991). An effective population size of 50 is considered the minimum necessary to maintain sufficient genetic variability, while 500 individuals are required to offset effective drift (Barrett & Kohn, 1991; Brzosko, Wróblewska, & Talalaj, 2004; Honnay, et al., 2005; Jacquemyn et al., 2012; Tremblay & Ackerman, 2001). Despite the massive decline in abundance caused by drastic environmental change, habitat destruction, and illegal collection, the populations of L. speciosa still have moderate to high values of effective population size. The wide distribution range and the cross-pollinated syndrome of L. speciosa seem to have been partly responsible for its high genetic diversity and moderate-to-high effective population size.

Genetic Structure and Bayesian Admixture Analysis

The AMOVA showed that most of the variation in L. speciosa occurs within populations. This is consistent with the results obtained by Da Cruz et al. (2011); based on the genetic structure of plants, Li and Ge (2006) suggested that cross-pollinated plants tend to have a greater diversity within populations and very little between them.

The distribution of ancestry proportions clearly showed a substantial connectivity between populations across the Cuitzeo basin. Nevertheless, we identified breaks in the continuity of gene flow such as the first barrier, which separates populations located in the southwestern from those in the central part of the Cuitzeo basin. We also observed some populations that have become isolated by a complex barrier; this was also evident in the map of oak forests remaining in this region (Figure 4). Geographical restrictions and a low gene flow between populations were also indicated by the second barrier that separate populations located in the northwestern from those in the central part of the basin. This pattern reflects the effect of habitat fragmentation on the patchy distribution of populations across the Cuitzeo basin. Pine and oak trees have almost disappeared from some localities of this basin due to their value as timber for construction or as a source of charcoal, an important economic activity in the region (Aguilar, Guilardi, Vega, Skutsch, & Oyama, 2013).

Small forest fragments that might still remain between larger populations would possibly function as vegetation corridors where isolated trees could function as stepping stones for gene flow between populations, as is the case in Quercus castanea in the same basin (Herrera-Arroyo et al., 2013). The history of L. speciosa in the Cuitzeo basin points to the fact that both fragmented and continuous populations have been recently subject to moderate reductions in effective population size, which have been large enough to affect allelic diversity (Barrett & Kohn, 1991; Lowe et al., 2005). Such reductions may have been caused by episodic environmental fluctuations (e.g., drought or geographical restrictions) or resulted from the recent founding of some of the populations (Chung et al., 2014; Forrest, Hollingsworth, Hollingsworth, Sydes, & Bateman, 2004; Murren, 2003). One likely explanation could be the effect of historic long-distance dispersal, which can maintain diversity through colonization. Thus, diversity is maintained at the regional rather than at the population level, with a moderate differentiation between populations (Bialozyt, Ziegenhagen, & Petit, 2006; Trapnell & Hamrick, 2004). Exponential population growth in newly colonized sites can also preserve genetic diversity near the leading edge, as rare alleles are less likely to be lost to drift (Excoffier, Foll, & Petit, 2009). Loveless and Hamrick (1984) suggest that genetic differentiation is the result of the ability of the species to disperse pollen and seeds; this holds true in orchids as they are wind dispersed (Chung et al., 2007; Gustafsson & Sjögren-Gulve, 2002; Li & Ge, 2006) and are capable of long-distance dispersal (Ackerman & Ward, 1999; Arditti & Ghani, 2000; Neiland & Wilcock, 1998). Therefore, it is quite possible for orchids to disperse across considerable distances of up to several kilometers (Flores-Palacios & García-Franco, 2003; Sharma, Clements, & Jones, 2000; Trapnell & Hamrick, 2004). Furthermore, Bombus pennsylvanicus, B. sonorus, and B. ephippiatus have been documented as potential pollinators of L. speciosa (Van der Pijl & Dodson, 1996). Other studies have concluded that, in orchids, population genetic structure is strongly influenced by the behavior of pollinators (Chung et al., 2014; Cozzolino & Widmer, 2005; Neiland & Wilcock, 1998; Qian et al., 2013; Sharma, Jones, & French, 2003; Tremblay, Ackerman, Zimmerman, & Calvo, 2005).