Life history and past demography maintain genetic structure, outcrossing rate, contemporary pollen gene flow of an understory herb in a highly fragmented rainforest

Study system

Aphelandra aurantiaca (Scheidw.) Lindl. is an understory herb of neotropical rainforests from southern Mexico to Bolivia. It is a self-compatible species that bears inflorescences with yellow floral buds that turn red when flowers open and produce nectar. In the rainforest of Los Tuxtlas in Mexico, this species is pollinated by the hummingbird Phaethornis longirostris (I Ramírez-Lucho, P Suárez-Montes & J Núñez-Farfán, pers. obs., 2013) altough it is also visited by butterflies and bumble-bees (Calvo-Irabién, 1989; Islas Luna, 1995). Seed dispersal is ballistic, ranging from 1 to 8.5 m (modal value of 1.5 m) from the maternal plant. This herb also exhibit vegetative reproduction by stolons. Its life span ranges from 13 to 18 years (Calvo-Irabién, 1989). Reduction in plant abundance of A. aurantiaca is related to the forest regeneration cycle where light is the most variable abiotic factor; the species inhabits both shaded forest understory and forest light-gaps (Calvo-Irabién, 1989; Calvo-Irabién, 1997).

Study site and data collection

The study was carried out at Los Tuxtlas Biosphere Reserve in southern Mexico, which constitutes the northernmost distributional limit of tropical rainforest in the Americas (Dirzo & Miranda, 1991). The region has lost more than 90% of its original forest cover in the past fifty years. Nowadays, the current landscape is composed of areas used for human settlement (1.27%), roads (0.78%), water bodies (1.92%), cattle ranching and crops (42.82%), riparian strips (4.29%), live fences (3.28%), isolated trees (1.03%), secondary vegetation of rainforest (0.53%), and fragments of rainforest (23.24%). Rainforest fragments are relatively small (<100 ha), surrounded by grassland and located in lowlands or restricted to the top of the mountains, in glens or areas of difficult access. At higher elevations (>600 m a.s.l.) cloud forest (4.61%) and secondary cloud forests (0.25%) is the predominant vegetation (Dirzo & Garcia, 1992; CONANP, 2011; see Salazar Arteaga, 2015) (Fig. 1).

Twelve populations of A. aurantiaca were sampled throughout Los Tuxtlas rainforest (Table 1). Preserved areas covered by rainforest are mostly surrounded by a matrix of pasture lands used for cattle ranching. Because the study area is highly fragmented, the choice of sampling sites was based on size area, accessibility, and the possibility of getting large sample sizes for genetic analyses (∼30 individuals). Forest fragment size best explains the differences in composition and plant structure (Arroyo-Rodríguez & Mandujano, 2006) therefore, we classified fragments as small (≤10 ha), medium (20–120 ha) and large and undisturbed fragments (>640 ha) (Fig. 1, Table S1). Unlike medium and large forest fragments, forest structure of small fragments is characterized by the absence of large primary trees in the canopy and lower abundance of palms and some herb species in the understory, but high abundance of shade intolerant secondary species (see Arroyo-Rodríguez & Mandujano, 2006). The large, undisturbed and continuous forest is in the core of Los Tuxtlas Biosphere Reserve; however it is partially separated by a deforested area produced by an illegal invasion (see Fig. 1). We established four sites within the large fragment as different populations based on their distance from deforested areas, and geographic distance between sampled sites (from 4. 39 to 14 km).

In each fragment size, we collected young leaf tissue from 29 to 38 adult plants. Individuals were selected as reproductive if they had an inflorescence/infructescence, or scars of these on the stem. We collected leaf tissue from individuals separated at least three meters apart. In total, we sampled individuals of A. aurantiaca in five small (n = 153), three medium (n = 97) and four large fragments (n = 138) (Table 1).

Microsatellites analysis and PCR amplification

We amplified 11 polymorphic microsatellites specifically developed for A. aurantiaca to determine its genetic structure (Suárez-Montes, Tapia-López & Núñez-Farfán, 2015) (Table S2). DNA extraction, amplification, and laboratory setup are detailed in Suárez-Montes, Tapia-López & Núñez-Farfán (2015). We genotyped all sampled individuals and scored alleles using the software GeneMarker V.2.4.0 (SoftGenetics, State College, PA, USA).

Genetic diversity analyses

We used MICRO-CHECKER software (Van Oosterhout et al., 2004) to detect null alleles. Deviations from Hardy–Weinberg equilibrium (HWE) and linkage disequilibrium were tested using GenePop 4.2 (Rousset, 2008) and FSTAT v.2.9.3.2 (Goudet, 1995), respectively. We estimated descriptive statistics of genetic diversity including allelic richness (A) and expected and observed heterozygosity (He and Ho) using Arlequin v.3.5.1.3 (Excoffier, Laval & Schneider, 2005). The private allelic richness (P(1)) was calculated by rarefaction with H_P−RARE (Kalinowski, 2005). Mean inbreeding coefficient (F_IS) was estimated using GENETIX v.4.05 (Belkhir et al., 1996–2004), based on 10,000 permutations. The hypothesis of isolation-by-distance was tested by a Mantel test based on 999 replicates using the ade4 package (Dray & Dufour, 2007) in R 3.1.3 (R Development Core Team, 2015).

Population size and demographic history analyses

We estimated effective population sizes (N_e) using a linkage disequilibrium method in NeEstimator V2 (Do et al., 2014). We employed two approaches to detect changes in population size. The first approach is based on the detection of heterozygosity excess or deficiency in a very recent period of time (2N_e–4N_e generations) using the program Bottleneck 1.2.02 (Piry, Luikart & Cornuet, 1999). Populations recently bottlenecked would lose rare alleles faster than heterozygosity, resulting in an apparent heterozygosity excess in comparison with a population at equilibrium (Heq) (Cornuet & Luikart, 1997; Piry, Luikart & Cornuet, 1999).

The second approach employs coalescent simulations in an approximate Bayesian computation (ABC) framework to infer past demographic history, as implemented in the software DIYABC 2.0 (Cornuet et al., 2008). ABC chooses a demographic scenario that best fits the observed data by running simulations constrained by the specifications of the model (e.g., a demographic bottleneck). It approaches the posterior probability distributions of parameters by selecting the simulated datasets with the smallest Euclidian distances to the observed data, as measured by summary statistics (Cornuet et al., 2008; Cornuet et al., 2014). We compared four demographic scenarios: two with constant N_e, one with a decline, and one with an expansion. We ran one million simulations for each scenario. The parameter settings and priors are shown in Table S3. The change time in N_e was set at 10–10,000 before present, assuming generation time of one year for the species (Calvo-Irabién, 1989). We used a generalized stepwise mutation model with a mutation rate of 10⁻⁴–10⁻³. The summary statistics were: mean number of alleles, mean heterozygosity, and mean allelic size variance. We assessed the fit of the models to the data by principal components analysis (PCA), implemented in DIYABC. We estimated the posterior probabilities using a logistic regression approach on the first 1% simulations closest to the observed dataset. To check the confidence of model choice we estimated type I and type II error rates by simulating 500 pseudo-observed data sets (Cornuet et al., 2008; Cornuet et al., 2014).

Genetic structure and clustering analyses

We estimated Rst and a hierarchical analysis of molecular variance (AMOVA) using Arlequin v.3.5.1.3 program (Excoffier, Laval & Schneider, 2005). Population structure was explored with both model based (Structure) (Pritchard, 2010) and distance based approaches (DAPC, Pop Graph and NetStruct) (Jombart, Devillard & Balloux, 2010; Greenbaum, Templeton & Bar-David, 2016).

The Bayesian clustering algorithm implemented in Structure V2.3 (Pritchard, 2010) estimates the probability of genotypes being distributed into K number of clusters. Simulations were run using correlated allele frequencies, under admixture ancestry models, conducting a burn-in of 10⁶, MCMC iterations of 10⁶, and K varying from 1 to 12. The total number of clusters (K) was inferred with the Evanno ΔK method (Evanno, Regnaut & Goudet, 2005) using STRUCTURE HARVESTER (Earl & vonholdt, 2012) and CLUMPAK pipeline (Cluster Markov Packager across K) (Kopelman et al., 2015) to visualize bar plots.

NetStruct 1.2 program (Greenbaum, Templeton & Bar-David, 2016) infers genetic structure using network theory. The equivalent of a genetic population structure is the community partition of a network constructed with individuals as nodes and edges (paired connections of nodes), defined by using a similarity measure. Clustering is done by locating groups of nodes (community) that are strongly connected within the group but weakly connected to nodes outside the group. The strength association (SA) measures how strongly the individuals are related to the community at which they were assigned to. The strength association distribution (SAD) analysis examines the distribution of SA values of different communities and provides information about relative gene flow. A narrow SAD indicates low gene flow, while left-skewed SAD suggests constant moderate gene flow; recent migrants will display low SA values, increasing the variance and left–skewness of the distribution. Finally, NetStruct evaluates the statistical significance of community partitions using permutation tests (Greenbaum, Templeton & Bar-David, 2016). We used the Fast Greedy algorithm with a medium threshold of 0.24 and 999 permutations for modularity significance test. We characterized the SAD of communities by the Coefficient of Variation (CV), as a measure of dispersion.

Membership probability analysis

To evaluate the membership probability of individuals to their population of origin we used the discriminant analysis of principal components (DAPC), a multivariate analysis implemented in the adegenet package (Jombart, 2008; Jombart, Devillard & Balloux, 2010) in R 3.1.3 (R Development Core Team, 2015). Based on the retained discriminant functions, the analysis derives membership probabilities for each individual of original source populations (Jombart & Collins, 2015). We followed adegenet directions for alpha scores optimization; we performed the analysis with 50 PCs retained. We also evaluated admixed individuals with no more than 0.5 probability of membership to any population.

Spatial structure and genetic barriers analyses

To visualize the spatial genetic structure and the connectivity across the landscape we used Population Graphs implemented in gstudio (Dyer & Nason, 2004) and popgraph packages (Dyer, 2009; Dyer, 2014) in R 2.15.3 (R Development Core Team, 2013). Population Graphs is a graph-theoretical approach where the total genetic variation is decomposed into a geometric interpretation of components within and among the strata, and then modified using conditional covariance to the minimal topological configuration. The genetic structure within population variance is represented as nodes that are connected by edges whose magnitude is proportional to their interpopulation variance (Dyer, 2015). Nodes and edges were then mapped on their spatial coordinates. We tested isolation across graph distance (IBGD) through physical and graph distances using Graph software of GeneticStudio software (Dyer, 2009; Dyer, 2014). To identify long distance dispersal events or restricted gene flow we compared pairwise physical distances with their corresponding pairwise edge lengths using Graph software (Dyer, 2009). Edges whose length is significantly longer than expected indicate long distance dispersal events while edges with length shorter than expected indicate limited dispersal across the landscape (Dyer, 2015).

Recent migration analyses

Recent bidirectional migration rates (in the last 2–5 generations) were estimated for paired populations using the program BayesAss v 3.0 (Wilson & Rannala, 2003). Because we expected reduced connectivity for smaller and more isolated fragments, we also estimated migration rates between different fragment sizes. BayesAss does not depend on Hardy–Weinberg equilibrium and estimates m as the fraction of individuals in population i that are migrants derived from population j (per generation) (Wilson & Rannala, 2003). To check for consistency, we performed 10 runs with a different random seed number, and calculated their respective Bayesian deviance (Meirmans, 2014). The run with the lowest deviance value was used to select the best-fitting model. Each run consisted of three million iterations for the chain with an initial burn-in period of one million iterations, and interval between samples for 2,000 chains of MCMC.

To detect possible first-generation migrants and their source population we used GENECLASS2 (Piry et al., 2004). We used the likelihood-base estimator L_home using Paetkau et al.’s, 2004 algorithm. L_home, is the likelihood of the individual genotype within the population where it was sampled. GENECLASS2 assumes Hardy–Weinberg equilibrium and species, sexual reproduction. The probability of each individual to be encountered in a given population was estimated using 100,000 simulated individuals with a threshold value of 0.01.

Mating system and contemporary gene flow

We established a 50 × 20 m plots (200 m²) within four fragments, two medium size fragments and two populations from large fragments. We did not include small fragments due to the lack of enough samples of maternal families (mother and progeny). Plant density was characterized in each population by counting the number of flowering individuals. In each plot, we collected leaf tissue and mature infructescences of maternal plants (n = 33). We germinated seeds and collected the leaf tissue from the emergent seedlings for DNA extraction. The number of maternal families varied between populations and each family between 2 and 11 individual plants (Table 2). To estimate outcrossing rate (t) and pollen pool structure (Φ_FT) we used six highly polymorphic unlinked microsatellite loci (0432, 5409, 1233, 4536, 4483, and 1808) (Table S2) (Suárez-Montes, Tapia-López & Núñez-Farfán, 2015) and analyzed all members of each maternal family (6–10 families per population).

We estimated parental inbreeding coefficient (F), multilocus outcrossing rate (t_m), single-locus outcrossing rate (t_s), and biparental inbreeding due to mating among relatives (t_m − t_s) using MLTR (Ritland, 2002). Standard errors were derived from 1,000 bootstraps. To estimate the confidence interval (CI at 95%) we used the estimated MLTR mean ± 1.96 × 1 sd. We also estimated the CI for F, t_m,  and t_s by the bootstraping in MLTR (Ritland, 2002). Progeny inbreeding was calculated with GENETIX v4.05, based on 10,000 permutations (Belkhir et al., 1996–2004).

We estimated the differentiation of pollen pools (Φ_FT) sampled by different maternal families using the TwoGener method (Austerlitz & Smouse, 2002) as implemented in GenALEx 6.502 (Peakall & Smouse, 2012). This model assumes uniform individual distribution and accurate density estimates. If the analysis reveals high Φ_FT then pollen dispersal is restricted, suggesting that different mothers are sampling pollen from, at least partially, non-overlapping sets of fathers. To avoid overestimation of Φ_FT values due to parental inbreeding (F_p) and selfing (s), we used the formula

${{\mathrm{\Phi }}^{\prime }}_{FT}=\frac{{\mathrm{\Phi }}_{FT}}{1+F_p}$, and given s as 1 − t_m where t_m is the multilocus outcrossing rate, Φ′_FT transforms to ${\mathrm{\Phi }}^{″}=\frac{{2{\mathrm{\Phi }}^{\prime }}_{FT}-S^2}{2{\left(1-S\right)}^2}$ for the selfing rate. Finally the pollen dispersal distance (δ), the effective number of pollen donors $\left(N_{ep}=\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{\Phi }}_{FT}^{″}\right)$, and the effective pollination neighborhood area, $\left(A_{ep}\mathrm{=}\frac{N_{ep}}{d}\right)$, where d is the density of reproductive plants (Austerlitz & Smouse, 2001a; Austerlitz & Smouse, 2001b).

Genetic diversity and inbreeding

Overall, levels of mean population genetic diversity of A. aurantiaca at Los Tuxtlas rainforest were high (N_a = 5.5, H_o = 0.57, H_e = 0.61) and similar among fragments. High genetic diversity values are expected for long-lived perennial plants (H_o = 0.63, H_e = 0.68; Nybom, 2004). Inbreeding coefficient values were low (F_IS = 0.097, p = 0.00), ranging −0.066–0.29 for populations, and statistically similar between fragments of different size (range: 0.05 to 0.114), but slightly lower for large fragments (Table 1). Mantel’s test did not detect a relationship between genetic and geographical distances of all populations (r =  − 0.28, p = 0.94). Geographic distances between populations are given in Table S4.

Effective population size and bottleneck analysis

Effective population sizes (N_e) are relatively large in all fragment sizes (Table 1). Average effective population size of large fragments (N_e = 205.7 (CI [123–480.6])) is not significantly higher than in small fragments (N_e = 119.2 (CI [86.4–176.7])). However, effective population size tends to decrease with reduction of fragment size. The lowest estimated N_e corresponds to one small fragment (19.2 at fragment 3 SM) while the highest value (infinite N_e estimate) corresponds to a large fragment (9 Lrg).

Under the Stepwise Mutation Model (SMM) and Two-Phase Mutation Model (TPM), the bottleneck test failed to detect a recent genetic bottleneck. The proportion of heterozygotes observed (H_O) was lower than expected (Heq), suggesting an absence of recent genetic bottlenecks in fragments of all sizes (all p < 0.001). All fragments sizes exhibited significant allele deficiency (Table S5).

Demographic history (ABC)

Population expansion was the best scenario and had the highest posterior probability (p = 0.99, 95% CI [0.9997–0.9998]). The PCA representation exhibited a good recovery of the posterior predictive distribution and the observed data (Fig. S1). Under this model, we found evidence of a small population with an effective population size of approximately 1,180 individuals that expanded to a present effective population size of approximately 68,600 individuals. We estimated that this expansion occurred approximately 2690 years before present (Table S3). We found a type I error rate of 0.03, and a type II error rate of 0.02, indicating a statistical strength of 97% and a high degree of confidence for the population expansion scenario.

Population differentiation

Populations showed a lack of genetic differentiation despite isolation by fragmentation. The hierarchical AMOVA analysis indicated weak genetic differentiation among all fragments (R_st = 0.037, p = 0.00), with most genetic variance (96.2%) within populations, while only 3.8% of the variance was among populations within fragment size classes.

Bayesian statistical modeling for clustering implemented in STRUCTURE showed the most likely number of clusters at K = 2 (LnP = −11100.43) and K = 3 (LnP = −10965.43) (Fig. S2). For K = 2, cluster I includes almost all sampled populations while cluster II was composed of only two populations (3 SM and 4 SM). For K = 3, cluster I is also composed by almost all populations (1 SM, 2 SM, 5 SM, 6 Med, 7 Med, 9 Lrg, 11 Lrg, and 12 Lrg), while cluster II and III included two main populations (8 Med, 10 Lrg and 3 SM, 4 SM, respectively) (Fig. S2).

As in the STRUCTURE analysis, the constructed network in NetStruct also detected three groups or communities with significant community partitions (p < 0.05). The network indicates that these communities are dispersed over the landscape without a clear pattern (Fig. S3A). The assignment of individuals to communities showed that all are composed of individuals from all populations. However, community I was mainly composed by populations 10 Lrg, 8 Med, 6 Med, and 5 SM; community II includes 7 Med, 1 SM, 3 SM, 4 SM and 11 Lrg populations; and community III is composed by 9 Lrg, 12 Lrg and 2 SM. Moreover, the mean of SAD was low and similar for the three detected communities (I: 0.0012, II: 0.0010 and III: 0.0013), all with a wide left-skewed distribution suggesting moderate strength of association and constant gene flow levels (Fig. S3B). Community II showed slightly lower association without a very wide skewed tail in comparison with the other communities, suggesting gene flow in the past. The coefficient of variation was high for all communities, but the lowest was for community II.

DAPC showed that membership probabilities were higher for individuals in their home population, ranging from 46% to 77% (Table S6 and Fig. S4). Population 3 SM, and 4 SM had the highest membership probabilities for individuals in their home, suggesting higher isolation. DAPC detected 89 admixed individuals: 33 from small fragments, 18 from medium fragments, and 38 for large fragments. The highest number of admixed individuals belongs to population 11 Lrg (15) and population 5 SM (12). The lowest value was for population 3 SM with one admixed individual.

Spatial genetic structure, connectivity and barriers

The Population Graph consisted of 12 populations connected by 24 edges that exhibited a significant conditional covariance. The topology showed that all populations formed a single interconnected network, indicating gene dispersal (Fig. 2). No IBD (Isolation by distance) among populations was detected (Mantel Z = 67.8, p = 0.805). Conditional genetic distances (cGD) are shown in Table S4. We found extended edges (between 1 SM–2 SM, 1 SM–3 SM, 1 SM–11 Lrg, 3 SM–4 SM, and 5 SM–6 Med) whose lengths were longer than expected from the spatial distances, indicating long distance dispersal (Fig. 2). We also found compressed edges (between 2 SM–11 Lrg, 3 SM–12 Lrg, 4 SM–7 Med, 4 SM–8 Med, 4 SM–10 Lrg, 7 Med–8 Med, and 11 Lrg–12 Lrg) whose lengths were shorter than expected from the spatial distances, suggesting a reduced permeability of landscape.

Migration rates and first generation migrants

The BayesAss analysis showed that most gene flow (in the last 2–5 generations) occurred from large to both medium and small fragments (Table 3). However, the large–medium rate was higher than the large-small rate, suggesting higher connectivity between large and medium fragments. GeneClass2 identified 15 putative first-generation migrants out of a total of 388 individuals (p < 0.01) (Table S6). These results also indicated considerable gene dispersal from the largest fragment. We found that nine out of the 14 migrants derive from the large fragment, five reside in small fragments, two in medium fragments and two in large fragments. Three derived from medium fragments and reside in smaller populations. Two out of 14 migrants derived from small fragments, one resides in a medium fragment and the other in a large fragment (Table S6).

Mating system, pollen structure and pollen movement

Plant density of A. aurantiaca at Los Tuxtlas is variable among populations. Densities ranged from 0.42 to 3.18 individuals/m²; similar values have been reported by Calvo-Irabién (1989), who found a range from 0.38 to 3.32 individuals/m². The density of reproductive individuals in our study ranged from 0.18 to 1.6 individuals/m².

Aphelandra aurantiaca has a mixed mating system (t_m = 0.81), altough mating system estimators of each population are not related to fragment size (Table 2). There were no significant differences among populations in multilocus outcrossing rate, single outcrossing rate, and biparental inbreeding. Bootstrapped confidence intervals showed a significantly lower multilocus outcrossing rate in one population (11 Lrg). The AMOVA of gametes indicated that most genetic variance was contained within mothers (70–90%). The correlation of paternity estimates(Φ_FT), were high and significantly higher than zero for all populations (Table 2). Hence, the correlation of paternity estimates,Φ_FT, indicated restricted pollen dispersal (part of the offspring are full-sibs), with a relatively low number of effective pollen donors (N_ep), and short pollen distance movement in this species (Table 2). Some populations showed a high (45%) or a low (17–19%) fraction of siblings sharing the same father. The number of effective pollen donors (N_ep) was relatively low (1.92, on average) (Table 2). Average inbreeding values [F_is(s.d.)] were slightly lower for adult maternal plants [0.03(0.15)] than for seedlings [0.05(0.04)]. Furthermore, we detected biparental (or uniparental) inbreeding in all populations (t_m − t_s ranged from 0.02 to 0.24), that did not differ between populations. The effective pollination neighborhood ranged from less than 0.6 to 14.4 m², whereas pollen distance movement ranged from half a meter up to 1.1 m (Table 2). We suggest caution in the interpretation of TwoGener results due to uniform individual distribution assumptions.

Genetic structure and habitat fragmentation

Habitat fragmentation has not produced genetic differentiation or immediate reductions in genetic diversity of A. aurantiaca despite dozens of generations since the onset of fragmentation in the region of Los Tuxtlas. Regardless of fragment size, populations possess private alleles and high genetic diversity (H_e = 0.61), similar to those of long-lived perennial herbs (Nybom, 2004), but higher compared to other understory Acanthaceae plants [e.gr., Graptophyllum reticularum (H_e = 0.31), Graptophyllum ilicifolium (H_e = 0.43) and Ruellia nudiflora (I = 0.26) (Shapcott, 2007; Vargas-Mendoza et al., 2015). Effective population size estimates (N_e) showed that most populations are effectively large, suggesting that habitat fragmentation has not as yet reduced N_e enough to detect an impact on genetic diversity.

The weak genetic structure detected (R_st = 0.037) is supported by a number of analyses, suggesting that A. aurantiaca populations have remained genetically connected. Most genetic variation is contained among individuals within populations (96%) rather than between populations (3%) in fragments of different size. A similar result was obtained for other understory tropical herbs in the same region (Cuartas-Hernández & Núñez-Farfán, 2006; Suárez-Montes, Fornoni & Núñez-Farfán, 2011). Genetic clustering analyses revealed that most populations shared genetic information of one cluster or community, which explain the low genetic differentiation. However, small populations (3 SM and 4 SM) form a genetic group, suggesting ongoing isolation.

The high genetic diversity, large effective population size, and low genetic differentiation found in populations of A. aurantiaca could be related to historical processes of populations rather than with the present landscape configuration. Also, life history characteristics of A. aurantiaca, such as its mating system, generation time, vegetative reproduction, and overlapping generations, might help to diminish the impact of genetic drift, maintaining large effective population size, and buffering the loss of genetic diversity due to habitat fragmentation (Weidema, Magnussen & Philipp, 2000; Hailer, Helander & Folkestad, 2006; Breed, Christmas & Lowe, 2014; Pellegrino, Bellusci & Palermo, 2015).

Genetic structure, gene flow and past demographic change

Although low genetic differentiation and occasional long-dispersal events of A. aurantiaca were detected, we also found evidence of restricted gene flow. The contrasting results of restricted ecological dispersal of A. aurantiaca over short distances and low genetic structure could indicate a lack of population equilibrium under current demographic conditions. When populations suffer frequent extinction and re-colonization processes, low Fst values are expected if colonist individuals are drawn from distant populations. Besides, since dispersal could be highly variable through time, direct measures of dispersal could miss long distance dispersal events (Coyne et al., 1982; Slatkin, 1985; Slatkin, 1994; Whitlock & McCauley, 1990).

Historical data suggest ancient contraction-expansion of Los Tuxtlas rainforest. Contractions occurred during a period of low temperatures and humidity (from 20,000 to 12,000 years ago during the last glacial maximum (LGM)) followed by subsequent vegetation expansion events (Graham, 1975; Toledo, 1982; Haffer & Prance, 2001; Gutiérrez-Rodríguez, Ornelas & Rodríguez-Gómez, 2011). In A. aurantiaca, ABC analyses suggest a plausible past population expansion scenario at Los Tuxtlas around the end of the LGM, when warmer climatic conditions established. More recently, contraction-recolonization of Los Tuxtlas rainforest could also be related with volcanic activity (during the last 153 years ago) (Martin Del Pozzo, 1997; Guevara & Laborde, 2012), and with forest fragmentation due to human activities (only during the last 42 years) (Dirzo & Garcia, 1992). Further evidence indicates no relationship of geographic and genetic distances among populations of A. aurantiaca, suggesting a relatively recent origin. There is also evidence of ancient and recent population expansion for an abundant palm of the understory of Los Tuxtlas (J Juárez–Ramírez, 2015, unpublished data; Martínez-Ramos et al., 2016). Therefore, the low genetic structure of A. aurantiaca could be due to different historical processes at Los Tuxtlas rather than recent habitat fragmentation.

Current gene flow and habitat fragmentation

The highest rates of migrants and first-generation migrants derived from the largest and medium fragments, underscore the importance of relatively large forest patches to prevent genetic isolation. Moreover, pollinators may use a series of different fragment sizes to forage, helping to maintain connectivity across the landscape (Llorens et al., 2012; Volpe et al., 2014). Specifically, hummingbirds are effective pollinators that can fly across relatively large areas during their foraging routes, carrying pollen grains to individuals’ located far apart (Stouffer & Bierregaard, 1995; Kraemer, 2001). However, it is necessary to conduct specific studies to assess the effect of Los Tuxtlas forest fragmentation on the abundance and behavior of the hummingbird pollinator Phaethornis longirostris, and their consequences on reproductive output of A. aurantiaca.

Although both pollen and seed dispersal are relevant to the pattern of genetic structure in A. aurantiaca, analyses do not allow us to determine which process is the most important contributor to gene flow. Natural gene flow often follows a leptokurtic distribution, implying that most genes move over short distances and only a small fraction move over long distances. Pollen dispersal kernels are often short, resulting in self-pollination or gene exchange among closely related individuals (Betts et al., 2014; Ellstrand, 2014). However, even a small number of long distance migration events can suffice to reduce F_st. Unfortunately, these rare events are difficult to detect in field studies (Nathan et al., 2003; Mona et al., 2014). For A. aurantiaca, it is likely that current events of long distance dispersal contribute to maintain the landscape connectivity preventing genetic differentiation and increasing local genetic diversity. In contrast, restricted pollen (0.41–1.09 m) and seed dispersal (1.5 m; Calvo-Irabién, 1989) promote substructure within populations.

Compressed edges suggest reduced gene permeability among populations that are geographically close, even in the large fragment (Fig. 1). Within the large fragment only population 12 Lrg shows contact with three other populations, but one of them is a compressed edge (11 Lrg–12 Lrg). This could be related to the physical barrier imposed by the “Vigia” hill (ca 600 m a.s.l.) within the preserve, reducing gene dispersal. Thus, topography (elevation) of the landscape should be considered in future studies as a factor affecting gene flow.

Mating system and fragmentation effects

The species’ mating system is an important factor that affects the distribution of genetic diversity. Aphelandra aurantiaca is predominantly outcrosser (t_m = 0.81). The description of A. aurantiaca as a selfing species with a mixed mating system agrees with values found in other hummingbird-pollinated plants in the Neotropics (Wolowski et al., 2013). However, we found that its mating system is predominatly of outcrossing. A mixed-mating system can combine advantages of both reproductive strategies: outcrossing promotes genetic diversity when pollinators are abundant, while self-fertilization may ensure reproduction when pollinators are scarce or absent (Goodwillie, Kalisz & Eckert, 2005; Ruan & Teixeira da Silva, 2012).

Self compatible tropical herbs do not necessarily suffer of inbreeding because they may possess breeding system traits that promote outcrossing (McDade, 1985). In general, A. aurantiaca showed no signs of inbreeding, although some populations exhibited inbreeding (Table 1). This finding could be a consequence of the Wahlund effect caused by genetic structure within populations (Murren, 2003) due to limited seed/pollen dispersal within populations. For A. aurantiaca at Los Tuxtlas fragmented rainforest, maintaining a mixed mating system with a high outcrossing rate appears to help preventing the loss of genetic variation, as would be theoretically expected for smaller populations.

Forest fragmentation does not seem to affect the contemporary pollen dynamics of A. aurantiaca, as in other herbs at Los Tuxtlas (Cuartas-Hernández, Núñez-Farfán & Smouse, 2010). The pollen pool structure was variable; in some populations it was more restricted than in others. The effective pollination neighborhood estimated resulted smaller than 14 m², and plants received pollen from neighbors located, on average, within a radius of 1.1 m. This finding agrees with an estimation of pollen movement using fluorescent dyes (reported by Calvo-Irabién, 1989). The effective number of pollen donors of A. aurantiaca is relatively low (range 1.0–2.8), suggesting that the potential pollen donors contribute little to N_ep. The moderate biparental inbreeding in A. aurantiaca could be explained by the limited seed dispersal and mating among close relatives, which may be due to hummingbirds’ moving among relatively close plants (P Suárez-Montes, pers. obs., 2014). High plant density may contribute also to increase inbreeding and shorten pollen dispersal distance in A. aurantiaca. However, assessing whether high plant density reduces pollen dispersal makes necessary an extensive sampling of populations with different densities. We suggest caution when interpreting results of TwoGener given that it assumes uniform individual distribution.