Beyond isolation by distance: What best explains functional connectivity among populations of three sympatric plant species in an ancient terrestrial island system?

Understanding how landscape features affect gene flow is critical to connectivity conservation and restoration management. Here, we examined the relationship between functional connectivity (gene flow) and structural connectivity (area and spatial configuration of habitats) in three co‐occurring short‐range plant taxa in an ancient terrestrial island system.


| INTRODUC TI ON
Human activities affect most of the Earth's terrestrial systems (Hobbs, Higgs, & Harris, 2009) and have degraded billions of hectares (Gibbs & Salmon, 2015;Nkonya et al., 2016). These activities can result in habitat fragmentation and a reduction in population sizes of the plant and animal species affected. Habitat fragmentation often has adverse effects including the disruption of gene flow that may ultimately lead to the loss of allelic diversity and increased genetic divergence among populations (Ellstrand, 1992;Slatkin, 1987;Young, Boyle, & Brown, 1996). Fragmentation can also have adverse impacts on mating systems by a reduction in available mates and, in plants, increased selfing for self-compatible species. This leads to reduced levels of reproductive success and fitness costs to progeny via inbreeding depression (Aguilar, Ashworth, Galetto, & Aizen, 2006;Jacquemyn, DeMeester, Longejans, & Honnay, 2012). Therefore, information on connectivity among populations of impacted species is critical to mitigate disturbances. Such information should be collected prior to activities that cause habitat fragmentation but is rarely collected at all.
Landscape genetics provides the basis for examining the relationship between structural and functional connectivity. This allows insight into the evolution of plant populations such as how they diverge and ultimately speciate (Manel, Scwartz, Luikart, & Taberlet, 2003;Storfer et al., 2007). Information on the drivers of genetic structure in species is also important from a conservation perspective because it increases understanding of dispersal and the impacts of habitat fragmentation. This information can inform the development of management strategies required to maintain (or restore) connectivity and population viability (Avise, 2000). However, studies of landscape effects on gene flow in plants are commonly conducted on a single species, and findings should not be generalized to other species in the same location (Richardson, Brady, Wang, & Spear, 2016).
We examine spatial genetic patterns in three co-occurring, conservation priority plant species associated with banded iron formations (BIF) of an ancient terrestrial island system in south-western Australia (SWA). Such discrete rocky habitats are found globally (Porembski & Barthlott, 2000;Jacobi & Carmo, 2008;Gibson, Yates, & Dillon, 2010;Gibson, Meissener, Markey, & Thompson, 2012) and cover around 3% of the Earth's surface (Guillot & Hattori, 2013). Ironstone ranges, including the Canga of Brazil and the BIFs of SWA, are ancient features and hotspots for plant species diversity and endemism (Gibson et al., 2012(Gibson et al., , 2010Jacobi, Carmo, Vincent, & Stehman, 2007;Yates, Gibson, Petit, Dillon, & Palmer, 2011). The high species endemism characteristic of many of the world's terrestrial island systems provides a model for studying evolutionary patterns and processes (Byrne et al., 2019). However, there have been no studies on multiple plant species using a landscape genetics framework in these systems. Acquiring information on genetic connectivity is important because in both Brazil and SWA, ironstone ranges are mined for iron ore (Ye, 2008). This creates a challenge to conserve biodiversity in these areas while exploiting their ore reserves. Thus, information on dispersal could be used to define and assess mining impact, and assist restoration.
In this study, we have two specific aims and associated hypotheses: 1. To compare spatial genetic structure in three sympatric, shortrange BIF endemic plant taxa with different pollination systems, using Bayesian clustering. We hypothesized that the three study species will have different genetic patterns despite their co-occurrence and narrow geographical ranges and that each will require different management recommendations; and 2. To determine whether the variation in gene flow of the plant taxa studied can best be explained by the distance between populations (homogenous) or by connected corridors of suitable habitat (heterogeneous). We hypothesized that functional connectivity would not be spatially homogenous and could thus be significantly better explained by suitable habitat between populations than by geographical distance alone.

| Study system and field sampling
The BIF ranges of the Yilgarn Craton in SWA are isolated, ancient ranges within a predominantly flat landscape. Our study site is the Helena and Aurora Ranges (HAR), which are topographically complex, low-altitude BIFs rising about 700 m above sea level (asl) and extending over ~13 km in length (surface area 52 km 2 ) (Di Virgilio, Wardell-Johnson, Robinson, Temple-Smith, & Hesford, 2018 banded ironstone formations, circuit theory, gene flow, habitat conservation, IBD, inselbergs,   landscape connectivity, landscape genetics, restoration planning Threatened and Priority flora taxa, including endemic species, and restricted plant communities not represented in secure conservation reserves (Gibson et al., 2012(Gibson et al., , 2010. Our three study taxa are These taxa inhabit different niches within the HAR; A. adinophylla is the most widely distributed and found on the slopes and adjacent plains surrounding ironstone ridges (Figure 1d; Maslin, 1999); L. bungalbin and T. aphylla subsp. aphylla are confined to cliff tops and steep stony slopes (Figure 1e; Barrett, 2007;Butcher, 2007), although the former prefers south facing, shady aspects. The presence of an elaiosome on seed of each of the three taxa suggests
Each marker was amplified in a 6 µl reaction volume containing PCR buffer, Bioline Immolase DNA polymerase and dNTPs based on recommendations by Bioline, 1.5 mM MgCl 2 , 0.06 µM of M13labelled forward locus-specific primer, 0.13 µM of reverse locusspecific primer, 0.13 µM of fluorescently labelled M13 primer (as described by Schuelke, 2000) and 15 ng gDNA. The following PCR conditions were used: 94°C for 5 min followed by 11 cycles at 94°C for 30 s, 60°C for 45 s (dropping 0.5°C per cycle) and 72°C for 45 s; followed by 30 cycles at 94°C for 30 s, 55°C for 45 s and 72°C for Prior to analysis, data sets were checked for clonality using gen- -Haond & Belkhir, 2006), as per Millar, Byrne, and Coates (2010). Replicate multilocus genotypes found likely to have arisen by asexual reproduction were removed from subsequent analyses. We tested for linkage disequilibrium (LD) among loci using fstat v.2.9.3.2 (Goudet, 2002). Sequential Bonferroni corrections were applied to alpha values in the determination of significance to correct for multiple comparisons of LD (Rice, 1989). Departure from Hardy-Weinberg equilibrium was assessed for each locus and population by chi-square tests in genalex v.6.5 (Peakall & Smouse, 2012 Standard measures of genetic variation including observed and expected heterozygosity (H O , H E ) and private alleles (PA) (alleles found in only one population) were calculated using GenoDive (Meirmans & Van Tienderen, 2004). Allelic richness (AR) was also calculated and rarefied to the smallest population size using HP-Rare (Kalinowski, 2005) (n = 13 for A. adinophylla, n = 18 for L. bungalbin and n = 19 for T. aphylla subsp. aphylla).

| Population differentiation
Pairwise genetic differentiation F ST (Wright, 1965), a measure of past gene flow among populations, was calculated in fstat v.2.9.3.2 (Goudet, 2002) and used as a response variable. Gene flow among populations was also assessed by D EST estimated in smogd v.1.2.5 (Crawford, 2010). D EST is an alternative measure of allelic differentiation among populations not biased by the genetic diversity of the populations (Jost, 2008).

| Population size
Changes in population size can have a strong influence on genetic patterns, particularly on rates of genetic drift. To examine whether past changes in population sizes account for any genetic differences between species, populations were assessed for past reductions in population size (over 10s to 1000s of years) using bottleneck v.1.2.02 (Cornuet & Luikart, 1997). Of the three available tests, the Wilcoxon signed rank test was applied, because (a) the sign test has low statistical power; and (b) the standardized difference test requires data from 20 or more loci. We used the two-phase mutation model (TPM), which is intermediate between the stepwise mutation model (SMM) and the infinite allele model (IAM), because few microsatellite loci follow the strict (one-step) SMM (Di Rienzo et al., 1994). We ran the TPM simulation as 90% one-step mutations and 10% multistep changes.

| Bayesian clustering
We used a Bayesian clustering model to examine the level of population structure in each species. structure v.2.3.4 software (Pritchard, Stephens, & Donnelly, 2000) assigns individuals probabilistically to user-defined K populations to achieve Hardy-Weinberg and linkage disequilibrium within populations. structure was run using the admixture model, assuming correlated allele frequencies, and population information was used for all individuals. structure was run with 250,000 Markov Chain Monte Carlo iterations after a burn-in period of 100,000 iterations. Preliminary data analysis using PCoA of pairwise F ST values suggested <7 genetic clusters for each species.
Therefore, we modelled with K = 1 to K = 8, with 10 iterations of each K. Structure Harvester (Earl & vonHoldt, 2012) was used to infer an optimal K based on the method of Evanno, Regnaut, and Goudet (2005). The 10 runs of the optimal value of K were summarized using CLUMPP (Jakobsson & Rosenberg, 2007) with the Greedy algorithm and graphically displayed.

| Landscape variables
We sought to examine the influence of distance between plant populations (IBD) and habitat suitability on the functional connectivity (genetic variation) of our taxa. The IBD model was developed as a raster where every pixel had a cost of one to represent a completely flat landscape and thus ignores the effect of terrain on gene flow (Noguerales, Cordero, & Ortego, 2016). Species distribution models (SDM; habitat suitability) were generated using MaxEnt species distribution modelling software (Phillips, Anderson, & Schapire, 2006) using default parameters.
Inputs into the SDM were derived from a 2-m digital elevation model, constructed from the last returns of LiDAR data, comprising measures of morphometry (topographic position index-TPI; topographic roughness index-TRI), hydrology (topographic wetness index-TWI; SAGA wetness index-SWI) and annual solar radiation (ASR). All metrics were uncorrelated. The TPI compares the elevation of each cell in a digital elevation model to the mean elevation of a specified neighbourhood of 10 × 10m around that cell.
Positive TPI values represent locations that are higher than the average of their neighbourhood window (e.g., ridges), negative values are lower (e.g., valleys), and flat areas are close to zero (Guisan, Weiss, & Weiss, 1999). TRI was used to quantify relief heterogeneity, with higher values indicating greater landscape complexity (Riley, Gloria, & Elliot, 1999). The hydrological indices used were the topographic wetness index (TWI) and the SAGA wetness index Model accuracy was determined using a 20% subset where absence data were generated in equal proportions for each species using a random point generator and summarized by calculating the area under the curve (AUC) of respective receiver operating characteristic (ROC) curves (Fielding & Bell, 1997). We used three summary statistics to compare SDMs. Niche overlap, which considers the pairwise similarity of the suitability values, was computed using Schoener's D (Schoener, 1968;Warren, Glor, & Turelli, 2008).
It ranges from 0 (no overlap) to 1 (complete niche overlap). Range overlap was calculated as the number of grid cells in which both species are predicted to occur above a suitability of 0.5, divided by the minimum number of grid cells in which either species in predicted to be present. Area overlap was calculated as the amount of suitable habitat from one species that is present in the suitable habitat of a second species.

| Circuitscape modelling
The two models (IBD and SDM) for each taxa were used as input into Circuitscape 4.0 (McRae, Shah, & Mohapatra, 2013) to model connectivity of multiple pathways and thus represent the empirical costs of gene flow through a flat landscape and one that incorporates the landscape features of the study area. All inputs to Circuitscape were set to represent resistance, which required inversion of the SDM, as high habitat suitability was assumed to have low resistance (Nowakowski, DeWoody, Fagan, Willoughby, & Donnelly, 2015). Current was set to flow in any direction (eight neighbours). Models were output as ASCII grids for visualization, and pairwise resistance distance matrices for all combinations of all focal nodes (population centroids) were extracted for statistical analyses.

| Statistical analyses: Relationship between resistance distance and genetic distance
The explanatory power of the two pairwise distance matrices (IBD and SDM) was tested against both response variables (F ST and D EST ) using multiple matrix regressions with randomization (MMRR) with the ecodist 2.0.1 package (Goslee & Urban, 2007) in r 3.3.0 (R Core Team, 2016). The response matrix was permuted 100,000 times to determine the significance of regression coefficients (Lichstein, 2006). Noting that all variables were significant at α = 0.05, we subsequently partitioned the variance into its two components (IBD, SDM) using the varpart function in VEGAN (Oksanen et al., 2013).

| Genetic structure
We did not find any identical multilocus genotypes, and there was no evidence of linked loci after Bonferroni corrections. All loci were polymorphic in all populations, and there was no consistent departure from Hardy-Weinberg equilibrium for any locus across all populations (see Table S2.1-3 in Appendix S2). Replicate PCRs produced the same allele scores as the originals in ~ 99% of comparisons.
From MICROCHECKER, the average frequency of null alleles in A. adinophylla was < 5% across the whole data set but two loci (accur26 and accur5) had higher frequencies (10 and 8%, respectively). To determine the overall affect that these loci had on the results, all analyses were performed with these loci both present and absent. Excluding these loci reduced the overall average null frequency to only 1.7%, and neither F ST nor diversity metrics changed greatly. As such, all subsequent analyses were performed on the observed allele frequencies. For L. bungalbin and T. aphylla subsp.
aphylla, the average frequency of null alleles was < 5% across the whole data set.
Genetic diversity was high for the three species (see Table S2.4-6 in Appendix S2) but varied among populations. However, there was F I G U R E 2 Gene flow potential across the Helena and Aurora Ranges using circuit theory. Relationships between genetic differentiation and its most explanatory variable are paired with maps of connectivity between populations. These were isolation by distance for Acacia adinophylla (a, b); and habitat suitability for Tetratheca aphylla subsp. aphylla (c, d). The most explanatory variable for Lepidosperma bungalbin was inconclusive. We have chosen to show isolation by distance (e, f) as it is the most visually intuitive. Sampled population centroids are shown as white-filled triangles. Pie charts show the proportion of membership values from structure analysis for each population. Colours in pie charts indicate the mean proportion of assignment to clusters for all individuals within a population, where optimal number of clusters was defined by the Evanno et al. (2005) Table S2.7-9 in Appendix S2). Mean D EST was 0.187 ± 0.045. For T. aphylla subsp. aphylla, mean F ST was 0.031 ± 0.005 and mean D EST was 0.078 ± 0.027.
The STRUCTURE analysis and the method of Evanno et al. (2005) identified K = 2 as the optimal number of clusters to the data for A. adinophylla and T. aphylla subsp. aphylla (Figure 2). Geographical patterns of population clustering were generally similar. Populations in the south-west and north-east formed distinct clusters, with more admixed ancestry for centrally located populations. In contrast, analyses supported K = 3 as the optimal number of clusters in L. bungalbin and individuals were strongly assigned to each of the three multipopulation genetic clusters (Figure 2). Populations located in the south-west corner of the HAR formed a distinct cluster, populations in the central part of the range formed a second cluster and populations in the north-east of the species range formed a third cluster.

| Landscape genetic analysis
SDMs of the three species are shown in Supporting Information (see Figure S1.1 in Appendix S1), and a composite model is shown in Figure 3. All models were highly accurate (AUC = 0.96 -0.99).
We found that topographic control on hydrology (the SWI) was the most important variable for L. bungalbin (82%) and T. aphylla subsp.
The TRI was most important for A. adinophylla (55%). Solar radiation had some importance on the distribution of L. bungalbin (11%), which showed a preference for shade. Species response curves were presented in Robinson, Virgilio, Temple-Smith, Hesford, and Wardell-Johnson (2019).
A. adinophylla occupied the largest niche space of 2,332 ha (Table 1; Figure 3) and was positioned lower on the ranges than the other two species, exhibiting minimal niche (D = 0.27-0.35) or range overlap (0.07-0.10) between them (  Figure 3). Although considerably smaller in area (272 ha), and far more fragmented, the range of L. bungalbin overlapped the range of A. adinophylla on the southern side of the HAR (Figure 3), as reflected from the range overlap coefficient of 0.92 (251 ha; Table 1).
Coefficients of determination were significant for both variables for all three taxa using either single or multiple regression and required constrained ordination to examine the shared variance.
Variance partitioning identified that, except for a small unique contribution (R 2 = 0.07), the cost surface derived from the SDM was largely nested within IBD for A. adinophylla (Figure 4a- which was best explained by habitat suitability (R 2 = 0.56-0.59;

| D ISCUSS I ON
The persistence of disjunct populations of rare plant species impacted by anthropogenic activity may depend on gene flow between other populations to limit the loss of genetic variation through drift (Young et al., 1996). Thus, the development of strategies to maintain gene flow should identify structural features (e.g., habitat suitability) that facilitate functional connectivity. We studied the spatial genetic patterns and effect of structural connectivity on gene flow among populations of three short-range BIF endemic plant taxa in an ancient terrestrial island system in SWA. For each, we asked whether isolation by distance (IBD) and/or habitat suitability (SDM) best explained gene flow. We found different levels of population structure in these species and that they differed in response to these two explanatory variables.

| Population structure
structure analysis suggested that the number and location of genetic clusters were similar for A. adinophylla and T. aphylla subsp.
aphylla. In these taxa, geographically proximate populations were more similar than geographically distant ones and there were more mixed proportions of membership for centrally located populations, suggesting some influence of IBD. In contrast, structure clearly segregated populations of L. bungalbin in central, south-western and north-eastern sections of the HAR, and individuals were relatively strongly allocated to one of the three genetic clusters. In agreement, all of the structure models showed a divide between the north-eastern and south-western populations. IBD is the most likely explanation but another possibility is isolation by adaptation (IBA) whereby gene flow among populations is reduced by local genetic adaptation to different ecological characteristics (Orsini, Vanoverbeke, Swillen, Mergeay, & Meester, 2013). Nonetheless, we found no habitat differences between sites on either side of this divide or current disjunctions in the species distribution that could also explain the location of geographical clusters.
There are no direct estimates of pollen or seed dispersal distances for A. adinophylla, L. bungalbin and T. aphylla subsp. aphylla.
However, studies of congeners on BIF ranges in the region found genetic structure (examined using F ST and structure) at similar spatial scales to our study in A. woodmaniorum (Millar, Coates, & Byrne, 2013), A. karina (Funnekotter, Millar, Krauss, & Nevill, 2019) and L. sp. Mt Caudan (Binks, Millar, & Byrne, 2015), but much higher structure in T. paynterae subsp. paynterae (Butcher et al., 2011;Butcher, McNee, & Krauss, 2009 (Meirmans & Hedrick, 2011), and an alternative measure of genetic divergence, D EST , which is unaffected by genetic diversity levels (Jost, 2008;Meirmans & Hedrick, 2011), was also estimated. D EST indicated higher levels of differentiation for the study species than F ST , and a greater range of values, particularly for L. bungalbin, which also had the highest H E . Levels of population differentiation are comparable to the few studies of short-range Acacia and Lepidosperma species on BIFs (e.g., Millar et al., 2013;Binks et al., 2015;Funnekotter et al., 2019). However, T. aphylla subsp. aphylla had a much lower level of differentiation than T. paynterae subsp. paynterae (Butcher et al., 2011. Why spatial genetic structure was lower than expected in T. aphylla subsp. aphylla is difficult to explain given the lack of direct pollinator observations in this taxon. However, seed and most pollen dispersal is likely to be far more extensive than in T. paynterae subsp. paynterae. This finding indicates differences in processes that structure genetic variation in congeners with limited ranges and similar habitats. Furthermore, it supports the conclusions of previous studies in the region which have shown that species on or associated with BIFs show a broad range of spatial genetic structure that is not predictable based on life history traits alone (reviewed in Byrne et al., 2019).
Spatial genetic structure and population divergence are difficult to predict because they can arise from a variety of factors, including adaptation and range disjunctions, but gene flow through pollen and seed dispersal are key determinants in their establishment (Wright, 1943). However, seed dispersal is unlikely to be the major factor in the maintenance of connectivity among populations of our study species, many of which are separated by at least 100s of metres from the closest neighbouring population. All three species have seeds adapted for ant-mediated seed dispersal and ants are thought to disperse seeds over relatively small distances compared to other dispersal modes (Davidson & Morton, 1981;Gomez & Espadaler, 2013), although a recent study in the region suggests that seeds are often dispersed 10s of metres by ants (Pascov et al., 2015).
It is possible that some seeds may be dispersed by other modes (water, wind or other animals), but this is unlikely. Firstly, given a semi-arid climate and complex topography, it is questionable how seed could be dispersed long distances by water across ridges on the HAR. Secondly, seed morphology suggests that wind is unlikely to play a role in seed dispersal. Thirdly, seeds are likely to be quickly collected by ants, often within 24 hr of gravity dispersal (Majer, 1980(Majer, , 1984. Passerine species may be responsible for low levels of

| Landscape genetic analysis
We hypothesized that functional connectivity would be better explained by suitable habitat between populations than by geographical distance alone. This was an accurate postulate for T. aphylla subsp.
aphylla where connectivity derived from suitable habitat, which was almost entirely predicted from modelled water flow (SWI), explained up to 59% of the variation and geographical distance made no unique contribution. Accordingly, connectivity among populations of T. aphylla subsp. aphylla would most effectively be maintained through the retention of corridors of suitable habitat.
In contrast, gene flow among the study populations of A. adinophylla was best explained by geographical distance with a minor contribution from habitat suitability, which was otherwise nested within IBD. This species distribution was predicted predominantly from terrain ruggedness and water flow and occupies the greatest amount of suitable habitat of the three studied (2,332 ha vs. 272-1029 ha).
The relative abundance of unfragmented suitable habitat appears to have resulted in an overestimation of connected pathways, consequently reducing its explanatory value in the linear regression models.
Furthermore, whilst habitat may be suitable for A. adinophylla, it does not consider impedance from competing species, especially given that rugged habitats often exhibit higher species richness than physically simpler ones (Di Virgilio et al., 2018). Nonetheless, given that IBD was the most explanatory variable, modification of suitable habitat may therefore have less of an effect on gene flow for this species.
Relationships between the explanatory variables and L. bungalbin were poorer than for the other taxa and inconclusive. Suitable habitat was mainly derived from water flow and solar radiation (shade), and its range was nested within that of T. aphylla subsp. aphylla.
Given the strong response of genetic distance to SDM resistance for T. aphylla subsp. aphylla, along with its relatively small range (272 ha), we expected a similar response for L. bungalbin. This proved to be the case when regressed against F ST (R 2 = 0.43) but converted to IBD when regressed against D EST (R 2 = 0.42), confounding conclusions.
We know little about the pollinators of most plant species on BIFs in SWA even though their foraging characteristics and preferences for habitat are likely to be important determinants of connectivity.
Whilst there are no studies of the mating system of L. bungalbin, the species is likely to be wind-pollinated (Barrett, 2013). Therefore, pollinator behaviour is not important for connectivity. Rather, the distance and direction that wind can disperse pollen is relevant.
Therefore, we suggest that further research to model wind direction during flowering would provide useful insight.

| Management implications
We identified divergent population structures and different factors that explain functional connectivity in the three narrow range study species. There is long-standing debate on the suitability of using neutral markers for delineation of conservation units (e.g., Patkeau, 1999;Funk, McKay, Hohenlohe, & Allendorf, 2012;Hoelzel, Bruford, & Fleischer, 2019 aphylla, which both have lower levels of genetic differentiation than L. bungalbin. Our findings emphasize that connectivity is complex and difficult to predict. Whenever possible, multispecies studies of genetic structure and gene flow should be conducted to allow a better understanding of the effects of proposed developments on the future viability of gene flow and populations, and accordingly, whether species require separate management plans (Rouget, Cowling, Lombard, Knight, & Kerley, 2006). Our findings also highlight that direct studies of both pollen and seed dispersal should be undertaken to inform conservation planning. Findings from this study could also be used to guide restoration planning and monitoring. For example, the estimates of historical gene flow could be set as baselines and targets for connectivity following restoration (Proft, Jones, Johnson, & Burridge, 2018;Ritchie, Dyer, Nevill, Sinclair, & Krauss, 2019 The potential of any development activity to minimize impact on surrounding biodiversity will depend on the effort expended to inform the environmental impact assessment process, ideally through the collection of multiple sources of data (here, genetic, remote sensing and field-derived plant density) for the whole site, and for multiple species. Undertaking these activities in partnership with industry end users will inform stakeholder planning and help mitigate adverse effects on surrounding biodiversity.

ACK N OWLED G EM ENTS
This research was funded by Mineral Resources Ltd. We thank D.

DATA AVA I L A B I L I T Y S TAT E M E N T
All genotype data is available from the Dryad data repository. https :// doi.org/10.5061/dryad.h0t0dk5