Do plant populations on distinct inselbergs talk to each other? A case study of genetic connectivity of a bromeliad species in an Ocbil landscape

Abstract Here, we explore the historical and contemporaneous patterns of connectivity among Encholirium horridum populations located on granitic inselbergs in an Ocbil landscape within the Brazilian Atlantic Forest, using both nuclear and chloroplast microsatellite markers. Beyond to assess the E. horridum population genetic structure, we built species distribution models across four periods (current conditions, mid‐Holocene, Last Glacial Maximum [LGM], and Last Interglacial) and inferred putative dispersal corridors using a least‐cost path analysis to elucidate biogeographic patterns. Overall, high and significant genetic divergence was estimated among populations for both nuclear and plastid DNA (ΦST (n) = 0.463 and ΦST (plastid) = 0.961, respectively, p < .001). For nuclear genome, almost total absence of genetic admixture among populations and very low migration rates were evident, corroborating with the very low estimates of immigration and emigration rates observed among E. horridum populations. Based on the cpDNA results, putative dispersal routes in Sugar Loaf Land across cycles of climatic fluctuations in the Quaternary period revealed that the populations’ connectivity changed little during those events. Genetic analyses highlighted the low genetic connectivity and long‐term persistence of populations, and the founder effect and genetic drift seemed to have been very important processes that shaped the current diversity and genetic structure observed in both genomes. The genetic singularity of each population clearly shows the need for in situ conservation of all of them.

In coastal Atlantic Forest of eastern Brazil, there is a particular region named Sugar Loaf Land (de Paula, Forzza, Neri, Bueno, & Porembski, 2016), where a humid tropical climate combined with occasional tectonic activity in the past has created a "sea of hills" landscape (Salgado, Bueno, Diniz, & Marent, 2015). This area is inserted basically in the Araçuaí orogen (Alkmim, 2015), which has an internal zone dominated by high-grade gneiss and granitoid inselbergs that formed between 625 and 490 Ma (Gradim et al., 2014;Varajão & Alkmim, 2015). The exposure of these features could have started sometime around the beginning of the Neogene (c. 20-2.6 Ma), and more recent tectonic motions might have contributed to their final sculpture (Varajão & Alkmim, 2015).
Inselbergs are a special class of residual landforms with a landscape configuration featuring a contrast of prominent elevations (height > 100 m) with surrounding plains (Lima & Corrêa-Gomes, 2015). They occur isolated or in groups, separated by a few or many kilometers, forming inselberg landscapes. Ecologically, they are characterized by a range of harsh conditions (Lüttge, 2008), such as a high degree of insolation, high evaporation rates, very restricted local soil occurrence, and low water availability (Porembski, 2007;Szarzynski, 2000). These edaphic-climatic characteristics promote the occurrence of specialized vegetation comprising a great number of endemic species (Porembski, 2007). Specifically, Sugar Loaf Land comprises a specialized Bromeliaceae flora, influenced by climatic factors and regional species pools, and seems to have played a pivotal role in the speciation of many rupicolous members (de Paula et al., 2016). The isolated nature of these terrestrial islands is one of the characteristics that make them interesting models for studies of the genetic connectivity of populations (Porembski, 2007;Porembski & Barthlott, 2000). Strong genetic structure has been repeatedly verified for species in this environment, which indicates an association between low gene flow and high isolation among populations (Barbará, Martinelli, Fay, Mayo, & Lexer, 2007;Boisselier-Dubayle, Leblois, Samadi, Lambourdière, & Sarthou, 2010;Byrne & Hopper, 2008;Millar, Coates, & Byrne, 2013;Palma-Silva et al., 2011;Pinheiro et al., 2011Pinheiro et al., , 2014Tapper et al., 2014a,b).
A few years ago Hopper (2009) proposed the Ocbil theory, a series of hypotheses that attempt to explain the evolution and ecology of biota on very old, climatically buffered, infertile landscapes (Ocbils), in contrast to Yodfels (i.e., young, often disturbed, fertile landscapes).
The author focused on the Southwest Australian Floristic Region, the Greater Cape Floristic Region of South Africa, and the Pantepui region of the Guyana Shield in South America as three of the most significant areas on Earth with Ocbils. However, this theory was recently reviewed and new likely Ocbil regions were pointed out by authors, including inselberg systems, in terrains where Ocbils are interspersed among Yodfels (Hopper, Silveira, & Fiedler, 2016).
is an endemic species of granitic inselbergs of Sugar Loaf Land (Forzza, 2005;de Paula et al., 2016), which makes it a good target species for evaluating the history of genetic connectivity among populations in this Ocbil landscape.
This study used eight species-specific nuclear microsatellite markers (nSSR, simple sequence repeat) and five chloroplast microsatellite markers (cpSSR) to investigate the genetic diversity and structure of E. horridum populations in an Ocbil landscape. The specific aims were to estimate levels of genetic variation and population differentiation, identify spatial population structure and dispersal barriers across populations, and elucidate historical and contemporaneous patterns of connectivity among E. horridum populations.

| Habitat and study species
Encholirium horridum occurs on inselbergs located from southern Bahia State to the northern Rio de Janeiro State, mainly in the states of Espírito Santo and Minas Gerais, in southeastern Brazil (Hmeljevski, Reis, & Forzza, 2015). The climate of the area of occurrence can be basically classified as Tropical wet-dry (Aw; Köppen-Geiger system), with an average temperature and annual precipitation of c. 24°C and 1,100 mm, respectively (pt.climate-data.org). The climate is characterized by distinct seasons, with most of the precipitation occurring in "summer" while the driest months are in "winter," which is when the species flowers.
This bromeliad species has abundant populations with thousands of individuals. A study on fine-scale spatial structure (SGS) and paternity analysis in one E. horridum population revealed weak SGS and restricted gene flow with pollination events that did not go beyond 45.5 m (Hmeljevski et al., 2015). Besides this, clonal growth was found to be a rare event, indicating that the species is monocarpic (Hmeljevski et al., 2015). This plant species has a cryptic self-incompatibility system and is pollinated by vertebrates (Hmeljevski, Wolowski, Forzza & Freitas 2017). Although abundant populations are observed for this bromeliad, it has been categorized as Endangered in the The Red Book of Brazilian Flora . This is mainly because of the degradation of inselbergs from mining granite, as Espírito Santo State is one of the main Brazilian centers of mining and processing ornamental rocks (Menezes & Sampaio, 2012). Other threats include predation of rosettes and immature fruits by goats and cattle and the invasion of exotic grasses that results in frequent fires (Hmeljevski et al., 2015).

| Sampling and DNA extraction
We collected leaf samples from 526 reproductive individuals of 11 natural populations covering the entire geographic range of E. horridum ( Figure 1, Table 1). On average, c. 48 and 20 individuals were sampled for nuclear (nDNA) and plastid (cpDNA) genome characterization, respectively (Table 2). Due to the relief's slope, the sampling strategy was to collect from the largest area within populations where it was possible to walk. Weak SGS has been reported for this species (Hmeljevski et al., 2015) so individuals from all E. horridum populations were sampled at intervals of at least 10 m to avoid potential sampling of relatives. Leaf material was stored over silica gel at −20°C until DNA extraction. DNA extraction was performed using a NucleoSpin Plant II (Macherey-Nagel, Düren, Germany) extraction kit.
Locus amplification was performed individually, and the fragments F I G U R E 1 Geographic distribution of Encholirium horridum L.B.Sm. sampled populations and respective plastid network (A). Pie charts reflect the frequency of occurrence of each haplotype in each population. Haplotype colors correspond to those shown in networks. In the statistical Median-Joining network, the haplotype frequencies are proportional to circle size. The number of mutations required to explain transitions among haplotypes is indicated along the lines connecting the haplotypes by cross hatches. Small black circles correspond to unsampled or extinct intermediate haplotypes were subsequently pool-plexed for allele sizing. Allele scoring was made in GeneMapper Software 4.0 (Applied Biosystems).

| Prior to genetic analysis
As initial analyses, gametic disequilibrium and null allele presence were tested before including the microsatellite set in the study. Loci that are included in analyses despite gross violations of these assumptions or high rates of error could lead to inaccurate and biased genetic estimates (Selkoe and Toonen, 2006). For nDNA data, linkage disequilibrium was tested for each E. horridum population with the program FSTAT (Goudet, 1995). In order to avoid false positives, the p-values (=.05) obtained were adjusted applying a sequential Bonferroni correction for multiple comparisons (Rice, 1989

| Genetic diversity analysis
For the nSSR data, genetic diversity parameters were estimated by the number of alleles (A), the allelic richness (R), the number of private alleles (A P ), the number of rare alleles (A R , alleles with frequency < 0.05), and expected heterozygosity (H e ) in Hardy-Weinberg equilibrium (HWE). Analyses were run using the program FSTAT (Goudet, 1995), except for the private alleles that were estimated in GeneAlEx6.5 (Peakall & Smouse, 2006, and allelic richness that were calculated in the R package, diveRsity (Keenan, McGinnity, Cross, Crozier, & Prodohl, 2013). In order to make the allelic richness independent from sample size, the rarefaction method was used to standardize it to the smallest sample size. The inbreeding coefficient (F IS ) was estimated considering the presence and absence of null alleles, and its significance (determined by 10,000 permutations across loci) was tested using the program SPAGeDi 1.5 (Hardy & Vekemans, 2002).
For cpDNA data, alleles at cpSSR were combined into haplotypes for each individual. From this, we characterized each population by the number of haplotypes detected and haplotypic diversity (h), calculated with the software CONTRIB (Petit, El Mousadik, & Pons, 1998).
Phylogeographical signal was also investigated for both classes of genetic markers by comparing measures of differentiation based on ordered versus unordered alleles (Hardy, Charbonnel, Fréville, & Heuertz, 2003;Pons and Petit, 1996). This reveals the importance of mutation relative to other causes of genetic differentiation (i.e., gene flow and divergence time) (Hardy et al., 2003). For nSSR, the global F ST and R ST values were compared to examine whether mutations were a significant cause of population differentiation using the R ST permutation procedure described by Hardy et al. (2003). Populations that have been isolated for a significant period of time would have accumulated mutations that contribute to their differentiation, and if mutation at the microsatellite loci has contributed significantly to this differentiation, the observed value of R ST should be higher than the permuted R ST (pR ST ). Conversely, if populations were isolated for a relatively short period of time (in comparison with the mutation rate) or if gene flow has been important, R ST would not be significantly different from F ST (Hardy et al., 2003). The test was implemented in the program SPAGeDi using 10,000 permutations across allele size. For cpSSRs, the global F ST and N ST indices of amongpopulation differentiation were computed in the program SPAGeDI to test whether N ST > F ST . For the estimation of N ST , the distance between haplotypes was calculated as the sum of their absolute length differences across all loci. We performed 10,000 permutations of rows and columns of the distance matrix between haplotypes. Such a significant relationship suggests that a phylogeographic signal occurred when distinct haplotypes found in the same population were related more closely (i.e., separated by fewer mutational events) on average than distinct haplotypes sampled in different populations (Duminil et al., 2010).
For the nSSRs and cpSSRs, we estimated the pairwise genetic distance and verified the isolation by distance among populations using SPAGeDi, and Mantel's test with 30,000 randomizations was performed using the program IBDWS 3.23 (Jensen, Bohonak, & Kelley, 2005). Bayesian clustering was also performed, in STRUCTURE 2.3.4, to assign individuals to genetic clusters (K) and to estimate admixture proportions (Q) for each individual (Pritchard, Stephens, & Donnelly, 2000). The analysis was performed under the assumption that the allele frequencies in different populations can be correlated with one another and that alleles carried at a particular locus by a particular individual originated in some unknown population (admixture model). Analyses were carried out with a burn-in period of 20,000 and run lengths of 200,000, and 10 iterations per K, ranging from 1 to 12, to confirm stabilization of the summary statistics (Pritchard et al., 2000). We first ran the dataset without considering sampling localities, followed by analysis considering this information. In order to estimate the appropriate number of populations, ΔK was estimated as an ad hoc quantity related to the second-order rate of change of the log probability of data with respect to the number of clusters, as proposed by Evanno et al. (2005). The most likely K number graphics and barplot figures were generated with the CLuMpaK server on the web (Kopelman, Mayzel, Jakobsson, Rosenberg, & Mayrose, 2015; http://clumpak.tau.ac.il/).

| Quantification of migration among populations
For the nDNA data, short-(the past one to three generations) and long-term gene flow (average over the past n generations, where n = the number of generations the populations have been at equilibrium) was estimated to examine the connectivity patterns of the population. A Bayesian approach was used to determine recent migration rates (m), which was implemented in BAYESASS 1.3 (Wilson & Rannala, 2003). Samples were run for 2.0 × 10 7 interactions with a burn-in period of 10 7 generations and sampled every 2,000 interactions. To estimate long-term gene flow, we used the program MIGRATE 3.2.6 (Beerli & Felsenstein, 1999. This program
Environmental data (19 standard BIOCLIM variables) were obtained for all geographical coordinates, at 1 km spatial resolution, from the WorldClim database (Hijmans, Cameron, Parra, Jones, & Jarvis, 2005), and cropped to the range of Brazil. To avoid overparameterization of SDM due to redundant variables (Dormann et al., 2013), the correlations between bioclimatic variables were assessed and those with presumed reduced biological relevance (r > .9) were eliminated, which left 14 for subsequent analysis (Table S1)

| Visualizing dispersal corridors
To understand how environmental conditions and genetic data contributed to shaping the distribution of E. horridum in Sugar Loaf Land across the time periods, we inferred corridors of highest dispersal probability, as described in Chan et al. (2011). According to Neto et al. (2015), during the beginning of the last Quaternary glaciation (LIG, 115,000-70,000 years ago), Brazil experienced a period of maximum cold and aridity that had strong repercussions on the vegetation.
Climatic fluctuations were frequent between 70 and 22 kyr BP, with alternating drier and wetter periods. The final cool and arid period began approximately 22,000 years ago (LGM) and lasted until approximately 14 kyr BP, when rainfall and temperatures increased again, which favored the return of forests. The warmest interglacial period appears to have occurred in the mid-Holocene (5,600-2,500 years ago), and after a short warm period, the temperatures decreased again and the climate became drier.
From the cpDNA haplotype network (Figure 1), we generated population connectivity maps by summing the least-cost path (LCP) among haplotypes from different localities. The LCP is based on the idea that resistance to movement can be mapped by assigning each cell in a map a relative "weighted distance" or "cost" of moving across that cell (Chan et al., 2011). The cell "cost" is determined by the habitat characteristics of the cell and is evaluated by calculating, for each cell, the cumulative weighted distance between the cell in question and two designated source areas. This analysis results in a map that shows the relative linkage value across the landscape (where routes through the landscape encounter more or fewer landscape barriers) between the two source areas (Singleton, Gaines, & Lehmkuhl, 2002). The combined interpretation of genetic population connectivity results and the LCP maps highlight areas and corridors that may have been a source population or historical refuge (Chan et al., 2011).
Using SDMtoolbox 1.1a (Brown, 2014) in ArcMap 10.2.2 (ESRI 2011), with the "Create Friction Layer > Invert SDM tool," we first inverted species distribution predictions that had previously been generated under current and palaeoclimatic models (mid-Holocene, LGM, and LIG). Inverted species distribution models produce friction layers with highest values in cells where probabilities of species occurrence are lowest (i.e., highest resistance to colonization). Next, the corridors that minimize the cost of dispersal between sampling localities by following paths of lowest friction were calculated. For these calculations, the "Least-Cost Corridors and Paths > Pairwise: All Sites tool" was used (default settings).

| RESULTS
In addition to the analysis of the independent segregation hypothesis that showed 6.5% of locus combinations significantly deviated from p = .05, no pairs of loci were found to be in significant genotypic disequilibrium after Bonferroni's correction (p < .0018). The presence of null alleles was detected for all nuclear loci and populations ranging from 0 to 4 and 0 to 6, when considering the number of loci per population and the number of populations per locus, respectively (Table S2).

| Nuclear and plastid genetic diversity
We identified a total of 193 alleles for the set of nSSRs. The mean number of alleles was 44 per population and ranged from 10 to 95 (Table 2). The observed values of allelic richness for the different populations varied from 1.24 to 11.05 (Table 2). The population with the highest allelic richness value was VP (11.05), followed by AD (8.97), CC (8.41), and MU (7.69). Private and rare alleles were detected in high frequency for the set of populations, which were 37% and 70%, respectively. Mean genetic diversity at the species level was equal to H e = 0.464, but with values highly variable among populations: 0.010 in CA to 0.833 in VP ( Table 2). The mean inbreeding coefficient per population was 0.206, reducing to 0.036 when estimated with exclusion of null alleles (Table 2).
For the cpSSRs, we identified from two to five alleles per locus, whose genetic combinations resulted in 11 different haplotypes (H1-H11, Table 2 and Figure 1). Although the high estimation of total haplotypic diversity (h = 0.951), the mean haplotypic diversity was low (h = 0.064) due to the high frequency of fixed haplotypes in populations (Table 2, Figure 1). In the haplotype network, it was possible to identify two haplogroups, one with eight haplotypes and the other with three haplotypes (Figure 1a). H3 identified in EC showing the higher number of mutational steps in relation to the other haplotypes ( Figure 1A).

| Distribution of genetic variability and phylogeographic pattern
High and significant genetic divergences were estimated among pop- High and significant isolation by distance for nSSRs were detected using F ST values (Z = 185.55, r = .4404, p < .01; Figure 2). However, no IBD was detected for cpDNA genome (Z = 6,971.24, r = −.0268, p = .4247). Bayesian clustering of groups of individuals revealed 12 and 11 genetic clusters (i.e., populations), considering no a priori population localities and with posterior information of the identity of populations, respectively (Figure 3, Fig. S1). The almost total absence of genetic admixture among populations was evident, especially in the K = 11 barplot (Figure 3).

| Estimates of gene flow
Values of immigration and emigration rates among populations based on the nSSRs were very low (Table 3), corroborating the strong genetic structure revealed from previous analyses (AMOVA and Bayesian clustering). Both short-and long-term estimates were not statistically significant among populations (Table S3). No pairs of populations showed significant asymmetric short-term estimates of gene flow, as seen by the overlapping of 95% confidence intervals. In fact, short-term estimate values were quite similar among the population set, varying from 0.0053 to 0.0074 (SD ± 0.0006). Estimates of longterm gene flow were more variable, ranging from 0.0254 to 0.0001 (±0.047), and only four population pairs showed asymmetric estimates (Table S3). Short-and long-term bidirectional estimates of gene flow across population pairs were not significantly correlated (Spearman's ρ = +.001, p = .95).

| Distribution and dispersal corridors across time periods
Based on the cpDNA results, putative dispersal routes in Sugar Loaf Land across cycles of interglaciation and glaciations events in the Quaternary period revealed that E. horridum population connectivity changed little during climatic fluctuations (Figure 4). Actual climatic conditions seem to be the most propitious to species genetic exchange. Areas with higher environmental suitability and presumable genetic connectivity correspond to inselberg landscapes in Espírito Santo State. These areas, characterized by high inselberg F I G U R E 2 Mantel's test between F ST /(1-F ST ) and Euclidian geographic distance log (Km) F I G U R E 3 Proportion of multilocus genotype of Encholirium horridum L.B.Sm. for K = 12 and K = 11 genetic groups detected by the software STRUCTURE. Distinct colors represent distinct genetic groups. See Table 1

| Patterns of genetic diversity and gene flow among inselbergs
Low gene flow, prolonged isolation, and persistence seems to be a common pattern to populations located in inselberg landscapes around the world (Barbará et al., 2007;Boisselier-Dubayle et al., 2010;Byrne & Hopper, 2008;Millar et al., 2013;Pinheiro et al., 2011Pinheiro et al., , 2014Tapper et al., 2014a,b). Random genetic drift and inbreeding have been suggested as main drivers of plant diversification and evolution on terrestrial islands due to limited connectivity among disjunct populations (Millar et al., 2013;Pinheiro et al., 2011Pinheiro et al., , 2014Tapper et al., 2014a). For example, in extreme cases of isolation, within-inselberg interspecific gene flow between congeneric species has been found to be higher than intraspecific connectivity between different inselbergs (Barbará et al., 2007;Palma-Silva et al., 2011). Also, persistence of species on granite outcrops indicates that these outcrops must provide heterogeneity and refugial opportunities (Tapper et al., 2014b) as long as there is a certain degree of connectivity among populations.
For instance, genetic diversity indicated that gene flow was sufficient to largely counteract any negative genetic effects of inbreeding and random genetic drift in inselberg populations of Acacia woodmaniorum (Millar et al., 2013).
Encholirium horridum showed, on average, median genetic diversity but demonstrated a great range in values among populations (Table 2). We observed a tendency for peripheral populations, such as GT, MA, and CA, to have a lower diversity. These populations are located on geographically isolated inselbergs in the landscape. Decline in genetic diversity accompanied by increased differentiation among populations from the center of a species' geographic distribution to its periphery is a frequent pattern across plants and animals (Eckert, Samis, & Lougheed, 2008). Isolated populations may experience declines in their sizes resulting in a loss of genetic diversity due to the effects of genetic drift and suffer from the Allee effects (Forsyth, 2003;Groom, 1998). Furthermore, the landscape morphology of northwestern Espírito Santo and western Minas Gerais is characterized by a concentration of inselbergs in chain formations. In these landscapes, E. horridum populations are found geographically closer to one another, potentially allowing more genetic connectivity and this is precisely where populations with the highest genetic diversity are located.
In face of the detection of high genetic diversity and private alleles in some E. horridum populations, and that strong genetic structure was detected in all of the analyses, we suggest there has been long-term persistence of the populations across time. The stability of potential dispersal routes also reinforces this pattern. Although no clear phylogeographic structure was evident in the results, considering the almost complete fixation of haplotypes, we also suggest that seed dispersal among inselbergs is an occasional phenomenon and the founder effect seems to be an important process in the formation of new populations. Furthermore, the relationships between AMOVA values estimated for both genomes was noteworthy (Φ ST(n) = 0.463, Φ ST(plastid) = 0.961; p < .001). Hamilton and Miller (2002)   inselbergs Pinheiro et al., 2014;Tapper et al., 2014b).

| Biogeographic patterns on Atlantic Forest granitic inselbergs
In Brazil, during the LGM, the retraction of seasonally dry tropical forests occurred and, consequently, Cerrado areas expanded until the mid-Holocene (Werneck et al., , 2012

| Encholirium horridum in the light of the Ocbil theory
Ocbils are defined across three continuously varying axes of landscape age, disturbance regimes, especially relative climatic buffering, and nutrient impoverishment of soils (Hopper, 2009;Hopper et al., 2016). Considering the key hypotheses of Ocbil theory (see Hopper, 2009), we can highlight some interesting points. The strong genetic structure found in both nuclear and plastidial genomes is in consonance with the first prediction, which enunciates the reduced dispersibility and strong differentiation of populations. Restricted pollen flow by paternity analysis was also detected in E. horridum (Hmeljevski et al., 2015). The second hypothesis, which predicts the elevated persistence of lineages and long-lived individuals, is supported by the evidence of the long persistence of populations across time based on genetic data and potential dispersal routes. Still, besides the absence of clonality in E. horridum, individuals are perennial and have slow growth (Hmeljevski et al., 2015). Other results, such as mean median genetic diversity of populations-probably purging of deleterious recessive genes (see Hmeljevski et al., 2015)-pollination mediated by vertebrates (Hmeljevski, 2013), and higher importance of pollen dispersal for genetic exchange than dispersal by seeds, are indicative of the occurrence of the James Effect (i.e., the retention of heterozygosity by selection mechanisms; third hypothesis). Lastly, the possible invasion of Atlantic Forest granitic inselbergs by routes from drier surrounding environments is in accordance with the fourth (Semiarid Cradle Hypothesis) prediction of the Ocbil theory. Therefore, based on this set of evidence, E. horridum is an Ocbil organism.
In the current scenario, genetic analyses highlighted the low genetic connectivity and long-term persistence of E. horridum popu- scholarships. We would like to thank also the Associate Editor and two anonymous reviewers for comments and suggestions.

CONFLICT OF INTEREST
None declared.