Niche specificity influences gene flow across fine‐scale habitat mosaics in Succulent Karoo plants

While the tempo of diversification in biodiversity hotspots has received much attention, the spatial scale of diversification has often been overlooked. Addressing this deficiency requires understanding the drivers of population divergence and the spatial scales at which they operate in species‐rich clades and ecosystems. South Africa's Succulent Karoo (SK) hotspot provides an excellent system for such research, being both compact (ca. 110,000 km2) and home to spectacular in‐situ radiations, such as the ruschioid Aizoaceae. Here we use GBS to document genetic structure in two co‐occurring ruschioid species, at both coarse (>10 km) and fine (<500 m) spatial scales. Where Ruschia burtoniae shows strong between‐population genetic differentiation and no gene flow, Conophytum calculus shows weak differentiation, with high levels of admixture suggesting recent or ongoing gene flow. Community analysis and transplant experiments reveal that R. burtoniae occupies a narrow, low‐pH edaphic niche, and at scales of a few hundred metres, areas of elevated genetic turnover correspond to patches of edaphically unsuitable habitat. In contrast, C. calculus occupies a broader niche and exhibits isolation‐by‐distance without a habitat effect. We suggest that edaphic specialisation, coupled with highly restricted seed and pollen dispersal in heterogeneous landscapes, has played a major role in driving rapid diversification at small spatial scales in this system. However, the contrasting patterns in our study species show that these factors do not influence all organisms uniformly, being strongly modulated by lineage‐specific traits that influence both the spatial scale of gene flow and habitat specificity.


| INTRODUC TI ON
Biodiversity hotspots are areas characterized by high concentrations (per unit area) of species richness and taxonomic endemism (Guilhaumon et al., 2008;Mazel et al., 2014;Mittermeier et al., 2011;Myers et al., 2000; but see Orme et al., 2005). Therefore, notwithstanding the critical role of long-term climatic stability in facilitating the local accumulation of endemic species (Harrison & Noss, 2017), interrogating the processes that promote diversification at small spatial scales is central to understanding the richness of these hotspots. Speciation requires either strong divergent selection or a reduction in gene flow between populations occurs (Coyne, 1992;Coyne & Orr, 2004;Slatkin, 1985;cf. Nosil, 2008), and so the challenge of understanding what drives speciation in hotspots reduces to understanding the factors-both environmental and organismalthat promote ecological divergence and/or restrict gene movement at small spatial scales.
Environmental heterogeneity influences gene flow both by presenting barriers of unsuitable habitat that impede between-population dispersal and so result in population isolation, and by presenting contrasting selective environments which simultaneously stimulate adaptation to local conditions and select against gene flow between the locally-adapted forms (Rundle & Nosil, 2005;Sobel et al., 2010).
The relative influence of these mechanisms varies from one system to the next, depending on the magnitude and spatial configuration of landscape heterogeneity. In oceanic archipelagos and high-elevation mountain systems, for example, population isolation by areas of unsuitable habitat (e.g., sea water or low-elevation habitat) commonly provides the initial impetus for differentiation leading to the formation of new species across broad spatial scales (Gillespie & Roderick, 2002;Grant, 1996;Steinbauer et al., 2016). In rainforest systems, by comparison, local adaptation across ecotones may be more important in powering speciation, the strength of selection being sufficient to counter gene movement (Brousseau et al., 2020;Fine et al., 2013;Kirschel et al., 2011;Misiewicz & Fine, 2014;Misiewicz et al., 2020;Moritz et al., 2000;Smith et al., 1997). The latter process can also operate at fine spatial scales, underpinning the "budding" speciation pattern observed in small, sympatric or parapatric populations of monkeyflower (Grossenbacher et al., 2014) and in twelve plant families in the California Floristic Province (Anacker & Strauss, 2014).Importantly, regardless of the mechanism involved, heterogeneity-driven divergence is not expected to yield a pattern of isolation by distance, the latter being more consistent with intrinsically limited short-range dispersal in a continuously-distributed population (Wright, 1943).
The spatial scale of gene flow is, of course, also influenced by intrinsic organismal features, particularly those influencing the spatial scale of gene transmission and the degree of niche specificity.
Seed and pollen dispersal distances in plants are influenced primarily by floral and fruit/propagule traits that determine pollinator/disperser type, as well as by other traits, such as plant height and seed mass (Boucher et al., 2017;Dick et al., 2008;Givnish, 2010;Schurr et al., 2005;Thomson et al., 2011;Vittoz & Engler, 2007). As such, traits affecting seed and pollen movement determine the scale of environmental heterogeneity that is relevant to gene flow and population genetic differentiation. This is further influenced by breeding system (outcrossing versus selfing), plant life span, plant habit (annual versus perennial; woody versus herbaceous), geographic range extent, and rare and unpredictable long-distance seed dispersal events (Hamrick & Godt, 1996;Hamrick et al., 1992;Nathan & Muller-Landau, 2000). Similarly, niche specificity is influenced by specific traits (e.g., plant habit, plant size, and environmental stress tolerance; Boucher et al., 2017;Boulangeat et al., 2012) which in turn determine the scale of environmental heterogeneity at which immigrant inviability limits gene flow. To sum up, the spatial scale of gene flow is influenced by organismal traits, and this undermines our expectation of a simple relationship between landscape heterogeneity and the spatial scale of speciation, which should instead be modulated by lineage-specific idiosyncrasies (Fine & Baraloto, 2016).
In the Cape biodiversity hotspot (i.e., Cape Floristic Region, CFR, sensu Goldblatt, 1978) of South Africa, complex topography, and the habitat heterogeneity that it entails, appears to underpin the exceptional species richness and local endemism of its flora (Goldblatt, 1978(Goldblatt, , 1997Goldblatt & Manning, 2002;Linder, 1985, landscape is the presence of a complex system of quartzitic gravel "islands", commonly termed quartz fields, which range from ca. 1->100 m across (Schmiedel & Jürgens, 1999) and are separated from each other by a sharply-contrasting intervening matrix of red sands (Schmiedel, 2002;Schmiedel et al., 2015). Since the high albedo of the quartz pebbles ameliorates the desiccating effect of the hot, dry summers, these quartz fields provide refuge to many species, in particular a charismatic suite of drought-intolerant, dwarf-succulent plants (Musil, Schmiedel, & Midgley, 2005, 2009Schmiedel & Jürgens, 2004), and the plant communities they support are distinct in composition and growth form from those of the intervening red sands (Schmiedel & Jürgens, 1999;Schmiedel et al., 2015).
Overall, the Knersvlakte quartz fields host 63 locally-endemic plant species, of which 37 belong to Aizoaceae (table 27 in Schmiedel, 2002). This richness is in part facilitated by the presence of environmental heterogeneity within quartz fields, particularly in edaphic properties. Particularly important is pH, which varies from neutral or slightly acidic to strongly acidic, with changes often being spatially abrupt Schmiedel & Jürgens, 1999;Schmiedel et al., 2015). This heterogeneity underlies striking changes in local species composition over tens of metres, producing a mosaic of distinct plant communities mostly dominated by Aizoaceae (Schmiedel, 2002;Schmiedel et al., 2015). These features make the Knersvlakte an excellent system for studying both coarse and finescale genetic patterns in plants and how they relate to environmental heterogeneity (Jelinski, 1997).
In this paper, we examine spatial patterns of genetic variation in two quartz field-endemic Aizoaceae species, Ruschia burtoniae L.
Bolus and Conophytum calculus (A. Berger) N.E.Br., in order to test the hypothesis that fine-scale spatial edaphic heterogeneity promotes fine-scale genetic differentiation in quartz field succulent species, and to elucidate the mechanism through which heterogeneity limits gene flow. For this purpose, we quantify patterns of genetic diversity at two spatial scales using single nucleotide polymorphism (SNP) variation obtained via genotyping-by-sequencing (GBS; Elshire et al., 2011). Using plant community and soil pH data, plus a common-garden comparison of each species' growth on different soils, we also characterize each species' edaphic preferences and, using this information, assess the extent to which genetic differentiation in each is a consequence of gene flow limitation associated with (a) barriers of edaphically-unsuitable habitat; or (b) divergent adaptation to soils of contrasting pH, or both.
In the absence of adaptive divergence, strong genetic structure would point to a primary role for edaphic barriers. We also test  for isolation by distance effects expected if genetic structuring is influenced primarily by intrinsic, and not extrinsic, limitations to gene flow.

| Study system
The study was conducted in the Knersvlakte region of the Succulent Karoo biodiversity hotspot (Figure 1). To broaden the scope of our findings, we selected species with contrasting morphological and ecological traits, which represent the broader variability that resides within Ruschioideae. Where R. burtoniae is an erect, succulent-leafed shrub growing to ca. 60 cm tall, having small, open, pink flowers borne on a many-branching cyme, and fruit capsules with complex seed-retention structures, C. calculus is a miniature succulent growing to no more than 10 cm tall with tubular, yellow, night-scented flowers borne singly in the axil of each leaf pair, and fruit capsules lacking complex seed-retention structures. Both species are largely confined to the Knersvlakte: R. burtoniae extends marginally into the hills to the north and east (Klak & Buys, 2012), while C. calculus is endemic excluding an outlying population located some 200 km to the north-east (Hammer, 2002) which is recognized as taxonomically distinct (subspecies vanzylii). Klak, Bruyns, et al. (2013) found Ruschia to be polyphyletic and that R. burtoniae falls outside the core Ruschia clade, being placed instead within a four-species clade alongside Ruschia phylicoides and the two species of Scopelogena. This clade is typified by a lack of closing bodies in the capsules and an unusual preference (in Aizoaceae) for relatively low-nutrient soils (Klak, 2000;Klak, Buys, et al., 2013;Klak & Buys, 2012) but is otherwise typical of the other woody ruschioid shrubs that dominate the SK flora. Conophytum, on the other hand, is a large, monophyletic genus with 108 species and 165 taxa (Hammer & Young, 2017;Klak, Bruyns, et al., 2013;Powell, 2016), of which 60% are endemic to the SK region and 28% have a total range size <10 km 2 (Young & Desmet, 2016). All are miniature leaf succulents and lack complex seed retention structures, but floral traits vary greatly in terms of morphology (from short and open to long and tubular), colour (white, yellow, pink and magenta) and circadian phenology (70% day-flowering, 30% night-flowering; Hammer & Young, 2017;Liede & Hammer, 1990). We selected four quartz field sites situated across the Knersvlakte at which both study species occurred ( Figure 1). These were: a site in the Groot Graafwater area of the central Knersvlakte ("GG"); a northern site in the Kareeberg area ("KB"); a southern site on the farm Quaggaskop ("QK"); and a south-western site in the Moedverloren area ("MV"). All sites except QK were within the Knersvlakte Nature Reserve. The exact locations of each site are withheld to combat poaching and can be obtained on request. The distances between GG, the central site, and the other sites ranged from 17 to 28 km, with the peripheral sites being 24-42 km apart. Each site spanned an area of less than 1 km 2 and had largely unbroken quartz cover.

| Sampling design
At each site, variation in community composition and associated soil chemistry was sampled in October 2014 using a series of circular plots, each of 5 m 2 area. The number of plots varied between sites, with 37, 35, 51 and 29 plots at GG, KB, MV and QK, respectively.
Plot locations were nonrandom, and were chosen to represent the diversity of quartz field plant communities present at each site, with C. calculus occurring in 13, 11, 13 and 7 plots, and R. burtoniae occurring in 26, 19, 19 and eight plots, at GG, KB, MV and QK, respectively. For each plot, a list of all Aizoaceae taxa (rank dependent on taxonomic uncertainty; see Schmiedel, 2002) present was noted and a soil sample was taken, combining soil from two to three points within the plot, at a depth of ≤10 cm. These samples were submitted to the Analytical Services facility at Elsenburg (Western Cape Department of Agriculture) for pH determination. This was done by sifting the soil to remove grains >1 mm and suspending in KCl, followed by measurement of pH using a standard pH meter.
For assessment of genetic patterns, 23-24 individuals of C. calculus and R. burtoniae were sampled from each site, giving a total of 95 samples per species across all four sites. Individual plants were sampled from across the species' distribution at each site, taking care wherever possible to sample individuals >5 m apart. Samples comprised leafy cuttings with leaves in good condition.

| Analysis of community and soil data
In order to characterize the ecological niches of the two focal species in a community context, we explored the relationship of community composition to soil pH across the four sites. Unless otherwise stated, community analysis was conducted in r v. 3.5.3 (R Core Team, 2018) using vegan v. 2.5-4 (Oksanen et al., 2019). Using only species present in ≥10 plots across all four sites in order to avoid biases introduced by rare species, between-plot Jaccard distances, based on species presence/absence data, were calculated using "vegdist". Thereafter, "varpart" was used to determine the proportion of variance in the matrix of Jaccard distances explained uniquely by pH (irrespective of site) and by site (irrespective of pH) using adjusted R 2 values. To visualise and assess the significance of the effect of pH on community composition, irrespective of site, we implemented the constrained ordination procedure introduced by Legendre and Anderson (1999) known as distance-based redundancy analysis (dbRDA) using the function "capscale". We specified site as a conditioning variable to ensure that the constrained axis (CAP1) corresponded to variance in community composition explained by pH after "partialling out" the effect of site, and that the unconstrained axes (denoted MDS1, MDS2, etc. and corresponding to orthogonal axes of variance in community composition not explained by pH) were also independent of site-specific variance. Significance was assessed using the randomisation procedure implemented in the function "ANOVA.
cca", with 999 randomisations. The niche characteristics of C. calculus and R. burtoniae were assessed visually using biplots of CAP1 and MDS1 showing both plot and species mean values. As a complementary analysis, we also assessed species co-occurrence in a pairwise manner and ignoring site and pH, using the method of Veech (2013) implemented in cooccur v. 1.3 (Griffith et al., 2016). This provided a probabilistic method to test whether pairs of species co-occur significantly more or less frequently than expected by chance.

| DNA extraction and sequencing
DNAs were extracted from freshly field-sampled leaf material. For C. calculus, we used the Qiagen DNEasy Plant Mini Kit (Qiagen, Hilden, Germany), following the manufacturer's recommendations, with minor modifications. Fresh leaves weighing ca. 1 g had their cuticle removed using a scalpel blade and were crushed using pestle and mortar, with gradual addition of the ligation buffer (Buffer AP1) to the point where the mixture could be poured easily. Approximately 0.1 g of PVP was also added to reduce the polyphenol content. For R. burtoniae, ca. 2 g of fresh material was ground into a fine powder using liquid nitrogen for use in the large-scale extraction protocol of Healey et al. (2014). This protocol modifies the widely used CTAB method of Doyle & Doyle (1987). Where necessary, samples for both species were cleaned and concentrated using the Zymo Genomic Clean and Concentrator kit (Zymo Research).
Samples were sent to the Cornell Institute of Genomic Diversity for GBS (Ithaca, New York) to be processed using a workflow which generates reduced representation libraries for Illumina sequencing following the laboratory protocol outlined in Elshire et al. (2011).
EcoT22I (restriction site ATGCAT) was used as the restriction enzyme for digestion as it produced fragments of mostly <500 bp with no evidence of repetitive DNA amplified by PCR. The enzyme PstI (CTGCAG) was also tested but was deemed inferior as it produced fragments ranging from 160 to 1,000 bp. The 95 samples of each species were run, alongside a blank, on separate lanes of a 100 bp single-end Illumina HiSeq 2000 instrument at the Cornell Core Laboratories Centre.

| Read alignment, SNP calling and SNP filtering
Read parsing, alignment, SNP calling and SNP filtering were all conducted separately for the two species. Raw reads were parsed, quality-filtered and trimmed using Trimmomatic (Bolger et al., 2014), and demultiplexed using the gbs-snp-crop pipeline v. 1.1 (Melo et al., 2016). Read alignment and variant calling were conducted following the ddocent pipeline v. 2.2.19 (Puritz et al., 2014). Reads were  (Rimmer et al., 2014), all with the default parameters. For each of the resultant variant call format (VCF) files a first round of filtering removed individuals with >60% missing data (five individuals of R. burtoniae and one of C. calculus), flagged variants, indels, nonbiallelic SNPs, SNPs with a minor allele count of one, genotypes with read depth <5, and SNPs with >10% missing data. For these steps we used VCFtools (Danecek et al., 2011) and vcflib v.
1.0.0 (Garrison, 2012). Of the remaining SNPs, we kept only those called by all three variant callers using bcftools v. 1.9 (http://samto ols.github.io/bcftools), using the freebayes VCF as the template for downstream analyses. Next, SNPs in putative linkage disequilibrium (LD) were removed using the sliding-window "indep-pairwise" function in plink v. 1.90b4 (Purcell et al., 2007), with pairwise LD R 2 calculated for SNPs within 50 bp windows, SNPs pruned until no pairs with R 2 > .2 remained, and the window shifted five SNPs after each step. Finally, loci present in putative repeat regions of the S. steenbokensis genome (SA Schlebusch, unpublished data) were removed using bedtools v. 2.28 (Quinlan, 2014).

| Analysis of genetic data
All analyses, unless otherwise indicated, were conducted in r v. 3.5.3 (R Core Team, 2018). VCF files were read into R using vcfr v. 1.9.0 (Knaus & Grünwald, 2017) and converted to genlight and genind objects native to adegenet v. 2.1.2 (Jombart & Ahmed, 2011) for analysis. Global between-site genetic differentiation and bootstrapped 95% confidence intervals were estimated using the F ST estimator of Weir and Cockerham (1984), which is appropriate for large SNP data sets (Willing et al., 2012), using the functions "wc" and "boot.vc" in hierfstat v. 0.4.22 (Goudet, 2005).
We implemented two methods to assess genetic structure at the between-site scale. Firstly, a principal component analysis was done using the "glPCA" function in adegenet. Secondly, individual assignment analysis was conducted using the sparse non-negative matrix factorization (sNMF) algorithms developed by Frichot et al. (2014), which are analogous to those used in other genetic clustering (i.e., "STRUCTURE-like"; Pritchard et al., 2000) programs, to produce individual ancestry coefficient estimates for the number of clusters (K) ranging from 2 to 8, with 100 runs for each K, using the R package lea v. 3.0.0 (Frichot & François, 2015). VCF file conversion to STRUCTURE format was done using pgdspider v. 2.1.1.5 (Lischer & Excoffier, 2011). The entropy-based cross-validation technique (implemented in lea) and PCA eigenvector plots were used to aid the choice of most likely K for each species, as recommended by Frichot and François (2015). The resulting set of Q-matrices of ancestry proportion estimates were summarised using the clumpak v.
1.1 online server (http://clump ak.tau.ac.il/; Kopelman et al., 2015), which identifies major and minor modes based on the degree of consensus between runs within each K. The resulting estimated ancestry coefficients were plotted using the R package pophelper v. 2.3.0 (Francis, 2017). Mantel tests (Diniz-Filho et al., 2013;Mantel, 1967) were used to test for isolation-by-distance (IBD) between the four study sites. For this purpose, Weir and Cockerham's F ST was used to estimate pairwise genetic differentiation between all site pairs (with 95% confidence intervals), using the package stampp v. 1.5.1 (Pembleton et al., 2013).
In addition, the function "genleastcost" in popgenreport v. 3.0.4 (Adamack & Gruber, 2014)  Two approaches were used to investigate spatial genetic structure between individuals within sites. Firstly, the correlation between log-transformed Euclidean geographic and genetic distance for each species at each site was assessed using Mantel tests calculated in vegan, with 9,999 permutations. Pairwise genetic distance matrices were estimated using "gd.kosman" in popgenreport v. 3.0.4 (Adamack & Gruber, 2014), which determines Kosman's genetic distance (Kosman & Leonard, 2005). Secondly, we implemented the "Mapping Averaged Pairwise Information" (MAPI) approach in mapi v. 1.0.1 (Piry et al., 2016) to visualise and test two-dimensional genetic divergence patterns for each species at each site, again using Kosman's genetic distance. Briefly, the MAPI approach entails the following: (a) a hexagonal grid is drawn so as to encompass the full geographic extent of all sampled individuals; (b) each pairwise genetic distance has its value assigned to an ellipsis connecting the pair; and (c) each grid cell is assigned the mean value of all ellipses that intersect it (this is denoted m w ). Then, to test the null hypothesis that m w values are randomly distributed in space, a null distribution of m w values is generated by randomly permuting the sampling localities 5,000 times, each time repeating the MAPI procedure. Using this null, the m w of each cell is assigned a probability of being either higher or lower than expected, and each p-value is subjected to the false discovery rate procedure of Benjamini and Yekutieli (2001) to account for positive spatial autocorrelation. Since the eccentricity of the ellipses may affect the outcome of the analysis, the above procedure was repeated for eccentricities of 0.8, 0.9, 0.975 and 0.99.
The results were qualitatively similar for different eccentricities (see Supporting Information) so we continued with the default value of 0.975.
To assess whether areas of unsuitable edaphic habitat, as determined on the basis of pH, corresponded to areas of high genetic turnover, genetic distance values were extracted from the MAPI plots for each point with community and pH data using raster v.
3.0.12 (Hijmans, 2020), and the relationship between pH and genetic distance was assessed using separate linear models for each species at each site.

| Common garden experiment
We conducted a common garden experiment to assess: (a) the po- We followed the approach of Paine et al. (2012) in using (mostly) nonlinear models to estimate separate growth curves for each group of interest, with 95% confidence intervals used to determine statistical significance. For each species, we fitted independent models for each soil type, relating size (i.e., height or diameter) measurement day, with pot and individual included as random effects to account for block effects and repeated measurements, respectively. Models were fitted using nlme v. 3.1-143 (Pinheiro et al., 2019), with logistic, Gompertz and linear link functions. The best-fitting models were chosen using pseudo-R 2 , calculated following Paine et al. (2012) for the logistic and Gompertz fits, and following Nakagawa et al. (2017) for the linear fits, in the package mumin v. 1.43.15 (Barton, 2019). In most cases, growth curves were asymptotic, for which logistic and Gompertz fits performed almost identically; therefore, logistic fits were preferred arbitrarily. In these cases, asymptotes were allowed to vary randomly by individual and pot. Where linear fits worked best, intercepts were allowed to vary randomly. We followed Paine 1.4.0 (Bengtsson, 2020). Some iterations resulted in models that did not converge, and these were excluded.

| Edaphic variation and community composition
Variance partitioning showed that both pH and site independently explained a substantial proportion of the total variance in Knersvlakte succulent community composition, though the effect of pH was much greater (pH|site: adj. R-squared = 0.151; site|pH: adj. R-squared = .0632). A nonparametric ANOVA of the dbRDA model suggested that pH was highly significant as a main effect after conditioning for the effect of site (F 1,147 = 25.9, p = .001) implying a constant effect of pH on community composition at all sites. Figure 2 shows the results of the dbRDA (see Table S1 for abbreviation key).
The constrained axis (CAP1), corresponding to variation in community composition that correlated with pH, explained 14.96% of the total variance in community composition after conditioning out the effect of site, while the first unconstrained axis (MDS1), representing variation not explained by site or pH, explained an additional 18.54%.
Ruschia burtoniae displayed the strongest association with low pH soils of any species, with the lowest CAP1 score (−1.72), a mean pH of 3.98 ± 0.50 SD and a pH range of 3.1 to 6.0 across plots in which it occurred (n = 72). In contrast, C. calculus displayed a weaker affinity for low pH soils, with a CAP1 score of -0.78, a mean pH of 4.22 ± 0.67 SD, and a pH range of 3.4 to 7.0 (n = 44). A two-sided t test revealed a significant pH difference between plots occupied by R. burtoniae and C. calculus (t = −2.01, df = 72.3, p = .0483). Ruschia burtoniae also had more unusual and F I G U R E 2 Ordination biplot of the distance-based redundancy analysis (dbRDA) of Aizoaceae community composition across the four study sites, with variation between sites conditioned out. The plot is scaled relative to species using Hill scaling. The proportion of variance in the first constrained axis (CAP1) explained by pH, and its direction, are indicated by the arrow. Colouring of points indicates plot pH. The first unconstrained axis (MDS1) represents variation in community composition not associated with site or pH. The percentage of variance explained by each axis is indicated in the axis titles. The two focal species are highlighted in larger text, and density kernels are drawn based on plots inhabited by Conophytum calculus (green shading) and Ruschia burtoniae (blue shading). See Table S1 for the key to abbreviated species names  (Table S2). Nonetheless, C. calculus did occur with R. burtoniae more frequently than expected by chance (observed = 28 plots, expected = 20.8, p = .008). Altogether these results indicate that habitat suitability decreases with increasing pH for both focal species, but that R. burtoniae is limited to a narrower and more extreme (low) pH range.  (Figure 3 and Figure S2) indicating that the four sites represent genetically distinct populations. The cross-validation

F I G U R E 3 (a) Principal component analysis biplots and (b) sNMF individual assignment results
showing estimated ancestry coefficients for the major mode identified by clumpak (summarised over 100 sNMF runs) at K = 4 for Ruschia burtoniae and Conophytum calculus. In the PCA plots (a), the proportion of variance explained by each axis is shown in parentheses, points represent individuals, shapes represent sites (as in b) and 95% confidence ellipses are drawn for each site. In the assignment plots (b), colours represent different genetic clusters, bars represent individuals (grouped by site), and the proportion of each colour in each bar represents the estimated ancestry coefficient for that cluster for that individual. For clarity, individuals are sorted by their estimated ancestry coefficient within each site. Note that clumpak placed all 100 runs at K = 4 under the major mode for R. burtoniae, whereas 65/100 runs at K = 4 were placed under the major mode for C. calculus. Note also that K = 4 was the best-supported K for R. burtoniae but K = 1 was best for C. calculus; we show K = 4 only for comparison [Colour figure can be viewed at wileyonlinelibrary.com] • • by entropy criterion analysis showed the lowest entropy at K = 4 ( Figure S4), and PCA eigenvectors showed a drop in explained variance from the first to fourth component axes, followed by a plateau ( Figure S5), also supporting four distinct clusters. In contrast, global F ST was 0.0090 (0.0078-0.0102 95% CI; p = .001) for C. calculus, which showed very limited spatial genetic structure, with both the assignment analysis and PCA recovering MV as somewhat distinct from the other three sites, which were not reliably distinguishable from each other (Figure 3 and Figure S3). For C. calculus, entropy showed a steady increase from K = 1 to K = 8, while PCA eigenvectors showed a steep drop in explained variance from the first to second axes ( Figures S4 and S5), reflecting the distinctness of MV and a lack of further structure in the data.
Evidence of between-site isolation-by-distance (IBD) also differed markedly between the species (Figure 4). Although differentiation was low in C. calculus, it showed a strong and significant pattern of IBD between sites. The stronger IBD effect for on-quartz distance was caused primarily by the MV-QK pair, which has an on-quartz distance almost twice that of its straight-line distance (see Figure 1).
In contrast, R. burtoniae showed no IBD pattern.
At the fine-scale (i.e., within sites, with between-individual distances <500 m) spatial genetic structure was remarkably strong in both species but considerably more pronounced in R. burtoniae.
Mantel tests (Table 1) Figure S6), with the MAPI analysis ( Figure 5) demonstrating that genetically diverged clusters of genetically similar individuals are separated by unsuitable, high pH habitat. In contrast, this pattern was never observed for C. calculus.

| Common garden experiment
The number of individual plants with repeated growth measurements in each soil type ranged from 10 to 19 for R. burtoniae (all with four repeated measurements), and from 10 to 16 for C. calculus (three to four repeated measurements). Growth profiles ( Figure 6) were explained well by logistic fits for all soil types for R. burtoniae C. calculus (matrix soil) relative to that in home soil ( Figure 6 and Figure S7). In addition, there was some evidence of a contracted growth period in QK soil relative to home soil. However, no differences in growth profile were observed in MV and KB soil. By contrast, C. calculus showed no evidence of reduced growth rate in any soil relative to home soil, with seedlings unexpectedly performing worse in home soil than in alternative soil ( Figure 6). The hypothesis that fine-scale genetic differentiation in the SK flora is linked to fine-scale patchiness of edaphic habitats appears to hold for R. burtoniae. Our results support previous work showing that pH strongly affects plant community composition in the Knersvlakte (Schmiedel & Jürgens, 1999), and the transplant experiments show that R. burtoniae seedlings experience reduced growth in unsuitable soils. Reduced seedling growth during the winter rains is likely to affect survival probability over the hot, dry summer as competition with adult plants intensifies (Milton, 1995;Musil et al., 2005Musil et al., , 2009, increasing the likelihood that unsuitable soils present a barrier to movement. Furthermore, because our common garden experimental approach controls for other potential soil-specific factors (e.g., herbivory and droughting) that may further reduce performance in unsuitable soils, the results probably more closely reflect the effect of edaphic variation alone.

| D ISCUSS I ON
Community and common garden transplant data indicate that R. burtoniae prefers strongly-leached, low-pH soils, a preference shared by the closely-related Scopelogena bruynsii and R. phylicoides (Klak, Bruyns, et al., 2013), both of which associate with fynbos vegetation growing on low-nutrient, acidic sands (Klak, Bruyns, et al., 2013;Klak & Buys, 2012). Since soils of this type, on the Knersvlakte, are typically associated with higher-lying areas, which are inselberg remnants of fluvial erosion processes that have dissected the Knersvlakte landscape since the Early Cretaceous (Kounov et al., 2008), their contemporary distribution is patchy.
The strong edaphic specificity of R. burtoniae explains both why genetic variation in this species is clustered at the fine spatial scale (i.e., within-site, between individuals <500 m apart), with genetic turnover at this scale being greatest across areas of unsuitable edaphic habitat, and why genetic differentiation shows no relationship to geographic distance at the coarse spatial scale (i.e., across sites, 17-42 km). The spatial scale of gene flow is not, however, a function of edaphic barriers alone, with seed and pollen dispersal distances influencing whether gene movement across such barriers is possible (Wiens, 1976). Our results support the hypothesis that the complex seed retention structures of R. burtoniae (and most Ruschioideae) restrict long-distance seed dispersal (Ihlenfeldt, 1994;Parolin, 2001Parolin, , 2006 and, in line with expectations for spring-flowering SK Aizoaceae that employ a generalist pollination syndrome (Ihlenfeldt, 1994;Scodanibbio, 2002), they suggest that pollination in R. burtoniae is localised.
In contrast to R. burtoniae, the edaphic heterogeneity hypothesis does not explain genetic structure in C. calculus. Instead, geographic distance appears to be the primary determinant of genetic divergence in this species. Whether the isolation-by-distance pattern observed in C. calculus is a product of ongoing gene flow or recent range expansion (Peter & Slatkin, 2013), both scenarios imply greater dispersal ability in C. calculus than in R. burtoniae. In part, this may be due to the weaker edaphic specificity of C. calculus compared with R. burtoniae, which may facilitate this species' movement both within and between quartz islands. Although our common garden results suggest that C. calculus is capable of growing in matrix soils, however, it is unlikely to do so in the wild as the windblown nature of these sands would result in seedling burial, and the lack of quartz would make them susceptible to overheating (Musil et al., 2005(Musil et al., , 2009Schmiedel & Jürgens, 2004 (Jürgens, 2004;Jürgens & Witt, 2014) are likely to be strong flyers and may be forced to search farther afield for nectar sources in late summer when little else is flowering (Hammer, 2002;Liede & Hammer, 1990). The flowers of this species are also scented F I G U R E 5 Maps to show spatial variation in pH and genetic distance. First column: maps of each site showing outlines of the MAPI grids for Ruschia burtoniae (blue) and Conophytum calculus (green), pH sampling points (pale blue points), interpolated pH values (colour shading), and areas of higher-than-expected genetic distance for R. burtoniae (hatched regions). Second and third columns: results of the MAPI procedure showing inferred mean pairwise genetic distances, genetic sampling locations (points), and areas of higher or lower than expected genetic turnover (black). GD, Kosman's genetic distance. Note the different scales for genetic distance for the two species [Colour figure can be viewed at wileyonlinelibrary.com] 100 m  (Jürgens, 2004) which is likely to enhance their attractiveness to pollinators (Jürgens & Witt, 2014;Parachnowitsch et al., 2012). While a comparison of maternally versus biparentally inherited markers would provide insight into the importance of seed versus pollen-mediated gene flow, we were unable to confidently recover adequate chloroplast markers from our data set to do this.
In addition to its effect as a barrier to seed and pollen dispersal, edaphic heterogeneity may promote genetic divergence via local adaptation, and the associated genetic hitchhiking/pleiotropy and selection against immigrants that this entails (Kawecki & Ebert, 2004;Rice & Hostert, 1993). Whereas C. calculus showed no evidence of reduced growth in soils from other sites, R. burtoniae exhibited reduced growth in QK relative to home soil which might be interpreted as evidence of local adaptation, though confirmation of this idea would require reciprocal transplants showing that QK plants perform poorly on GG soil. Importantly, however, no reduction in growth was observed in soil from the other two sites (KB and MV), suggesting that genetic differentiation in R. burtoniae may be driven more by landscape barriers to gene flow than by local edaphic adaptation. This inference is corroborated by the relative genetic similarity of GG to QK (compared to that of GG to KB and MV; Figures 3   and 4), the only site to show any potential signal of local adaptation.
While we acknowledge that soil variation represents just one axis of selective heterogeneity between sites, with local climate differences potentially being a stimulus for local adaptation, we consider the potential influence of climate to be limited given the relative uniformity of the landscape and the proximity of our sampling sites.
Nevertheless, the contrast between C. calculus and R. burtoniae suggests that reduced gene flow goes hand in hand with local adaptation. This accords well with the work of Ellis et al. (2007) showing that in the Knersvlakte-endemic dwarf succulent A. pearsonii, landscape barriers to gene flow are a prerequisite for genetic differentiation even in the context of strong divergent selection between sites.
Although R. burtoniae and Argyroderma show highly localized genetic differentiation, C. calculus reveals that these patterns do not automatically generalize to the entire SK flora, nor even its Aizoaceae component. This is because the genetic mobility of individual species is idiosyncratic, being dependent on a suite of species-and lineage-specific traits that dictate habitat specificity, seed dispersibility and pollen mobility. Given its minute stature, the relatively weak genetic differentiation shown by C. calculus, compared to that shown by the taller R. burtoniae, is somewhat surprising given that plant size has been identified as a potential driver of diversification rate/density (Boucher et al., 2017(Boucher et al., , 2020. However, the potential influence of other traits means that this relationship will often be noisy (Boucher et al., 2017) and certainly it does not hold here.
Given that minute stature has evolved several times independently in Aizoaceae (Klak, Bruyns, et al., 2013) and is prevalent in other important SK families, this observation brings into focus the question of whether and how minutism might drive diversification.
Inferring speciation mechanism from patterns of population-level genetic differentiation is not straightforward (Wiens, 2004). Given their genetic and even ecological distinctness, the four populations of R. burtoniae sampled in this study might represent distinct evolutionary species (De Queiroz, 2007;Wiley, 1978). In our experience, however, these populations show no discernible morphological divergence. Since most plant species, including those that make up the remarkable SK flora, have been differentiated on the basis of morphological distinctiveness (Bickford et al., 2007;Ellis et al., 2014), it is unclear how our understanding of cryptic population differentiation in R. burtoniae advances our understanding of the floristic radiation of the SK generally. One possibility is that genetic differentiation represents an important step towards the reproductive isolation of populations (Coyne & Orr, 1989;Moyle et al., 2004;Scopece et al., 2010) which are then free to diverge morphologically, via genetic drift and/or selection (Rieseberg & Willis, 2007), to the point where taxonomists would recognise them as separate species (Treurnicht et al., 2017). Another is that the process of population differentiation described here for R. burtoniae is atypical for SK plants, with most differentiation being adaptive and accompanied by morphological divergence. For example, it may be that erosional fragmentation of the Knersvlakte landscape (Kounov et al., 2008) fosters a more neutral mode of differentiation than is typical elsewhere in the SK.  (Liede & Hammer, 1990) and in being tetraploid (de Vos, 1947) where the majority of Conophytum, including C. calculus, are diploid (Hammer & Liede, 1992). Both C. breve and C. flavum typically occur on quartzitic substrates to the north of the Knersvlakte where they are broadly sympatric (Hammer & Young, 2017;Young & Rodgerson, 2017). This pattern is consistent with the idea that both allopatric/budding speciation (mediated by long-distance seed dispersal) and sympatric speciation (mediated by shifts in floral traits and/or ploidy) have contributed to diversification in Conophytum (Hammer & Liede, 1992;Powell, 2016;Powell et al., 2019). Given that there is floral divergence despite evidence of extensive gene flow between C. calculus populations, we encourage further research on speciation process in Conophytum, e.g., investigating the extent of gene flow between sister taxa with different pollination syndromes.
In conclusion, understanding the uneven global distribution of biodiversity requires that we understand the factors that promote differentiation at fine spatial scales. By focusing on patterns of genetic and ecological divergence in response to eco-spatial landscape features, we hope to gain insight into the early stages F I G U R E 6 Plots to show change in size over time for Ruschia burtoniae (left column) and Conophytum calculus (right column) seedlings germinated and grown in different soils in a greenhouse, with plots separated into within-and between-site comparisons (top and bottom rows, respectively). Within-site comparisons contrast growth in home (i.e., occupied) quartz soils with that in alternative (i.e., unoccupied, high pH) quartz soil and matrix (i.e., nonquartz) soil. Between-site comparisons contrast growth in home soils with that in occupied quartz soils from other sampling sites. Raw values per individual are shown as narrow dashed lines (coloured by soil type). Predictions of (non-) linear mixed-effects models are shown as thick lines with shaded areas representing 95% confidence intervals determined by bootstrapping Although our data corroborate existing studies Ellis et al., , 2007 in identifying soil heterogeneity as a potentially critical driver of fine-scale differentiation, we find that the scale of differentiation is influenced by other factors, notably organismal traits influencing niche specificity and seed and pollen dispersal distances. Future work examining additional species, and integrating genetic, community and trait data, is needed to assess the generality of the processes and patterns described here. In addition, there is a need for comparative studies that relate processes taking place at the population level to diversification pattern at the phylogenetic scale.

ACK N OWLED G EM ENTS
We thank the staff of the Knersvlakte CapeNature field office for providing accommodation and field assistance; Steven Hammer,