Evolutionary stability, landscape heterogeneity, and human land‐usage shape population genetic connectivity in the Cape Floristic Region biodiversity hotspot

Abstract As human‐induced change eliminates natural habitats, it impacts genetic diversity and population connectivity for local biodiversity. The South African Cape Floristic Region (CFR) is the most diverse extratropical area for plant biodiversity, and much of its habitat is protected as a UNESCO World Heritage site. There has long been great interest in explaining the underlying factors driving this unique diversity, especially as much of the CFR is endangered by urbanization and other anthropogenic activity. Here, we use a population and landscape genetic analysis of SNP data from the CFR endemic plant Leucadendron salignum or “common sunshine conebush” as a model to address the evolutionary and environmental factors shaping the vast CFR diversity. We found that high population structure, along with relatively deeper and older genealogies, is characteristic of the southwestern CFR, whereas low population structure and more recent lineage coalescence depict the eastern CFR. Population network analyses show genetic connectivity is facilitated in areas of lower elevation and higher seasonal precipitation. These population genetic signatures corroborate CFR species‐level patterns consistent with high Pleistocene biome stability and landscape heterogeneity in the southwest, but with coincident instability in the east. Finally, we also find evidence of human land‐usage as a significant gene flow barrier, especially in severely threatened lowlands where genetic connectivity has been historically the highest. These results help identify areas where conservation plans can prioritize protecting high genetic diversity threatened by contemporary human activities within this unique cultural UNESCO site.


| INTRODUC TI ON
Loss of biodiversity through human-induced habitat alteration has large impacts in reducing the sustainability of ecosystem services (Balmford & Bond, 2005;Flombaum & Sala, 2008;Sala et al., 2000). The impact of urbanization on biodiversity has received great attention (Alberti, 2015;Johnson & Munshi-South, 2017;Miles et al., 2019;Rivkin et al., 2018), yet anthropogenic activities in general have likely disturbed all natural areas (Kareiva et al., 2007). Abrupt changes due to human land-usage, including deforestation and agriculture, impact biodiversity long after recovery efforts (Jung et al., 2019;Moreno-Mateos et al., 2017;Newbold et al., 2015). Management efforts indicate that local factors are the best predictors of biodiversity, but that multi-scale landscape approaches are needed to marry historical and contemporary patterns with land-usage change (Auffret et al., 2018;Gonthier et al., 2014).
The details of evolutionary population demography, such as gene flow and population structure across temporal and spatial scales, are key to understanding anthropogenic impact on longterm biodiversity. As human-fragmented landscapes reduce gene flow, species can be challenged by lower genetic diversity within populations and increased genetic differentiation among them (Keyghobadi, 2007;McKinney, 2006), which over long-term can lead to increased inbreeding and local extinction (Charlesworth & Charlesworth, 1987;Frankham, 2005). While microevolutionary studies should directly incorporate anthropogenic impacts in modeling biodiversity dynamics Des Roches et al., 2020), how patterns of gene flow across both natural and anthropogenic features, especially across highly heterogeneous environments, influences conservation models is still a burgeoning field of research. Recent work emphasizes that challenges stem from the unknown effects of land-usage from large agricultural development to small but dense urban areas, the patchy habitat structure of the landscape, the entangled historical and contemporary effects, and interactions among variables (Ballare & Jha, 2020;Centeno-Cuadros et al., 2017;Fusco et al., 2020;Gonzalez-Serna et al., 2019;Miles, Johnson, et al., 2018;Thatte et al., 2020). Accounting for these challenges is vital to our biodiversity hotspots, where heterogeneity is often at its highest and information on population genetic connectivity is crucial for management (Myers et al., 2000).
The Cape Floristic Region (CFR), found in the southwestern tip of South Africa with ~9000 plant species in only 90,000 km 2 (Cowling et al., 1996;Goldblatt, 1997;Linder, 2003), 70% of which are endemic, is the most diverse extratropical area for plant biodiversity (Latimer et al., 2005;Linder, 2005;Rundel et al., 2016). The CFR is one of several Mediterranean Biomes, where high plant species diversity is typical , yet the southwestern CFR alone is an outlier, with endemism equal to or higher than the neotropics (Latimer et al., 2005;Rundel et al., 2016). The high levels of endemism are associated mostly with fynbos ("fine bush"), which is the dominant and most speciose vegetation in the CFR (Cowling et al., 1992), covering half of the geographic area with >80% of its species (Goldblatt, 1978;Linder, 2003).
The CFR is a global and local conservation priority (Rouget et al., 2003), with its status as a UNESCO World Heritage Site due to its ancient cultural history, and evidence of some of the earliest development of symbolic human behavior that likely evolved as a result of the resource-rich area (Marean, 2016;Marean et al., 2014). Today, the combination of a rich plant diversity and an equally rich cultural diversity in the CFR presents unique applications for biomedicine (Van Wyk, 2011). Unfortunately, human population growth in the CFR is the second highest among its Mediterranean Biome peers (Underwood et al., 2009), with agricultural and urban expansion devastating low-lying coastal areas and threatening inland mountain refugia (Bomhard et al., 2005). By 2030, it is expected that 75% of the South African human population will live in urban areas (Cooperative Governance & Traditional Affairs, 2016); thus, conservation efforts within the CFR need to link anthropogenic activity with eco-evolutionary dynamics Cilliers et al., 2014;Graham, 2017;Rebelo et al., 2011).
Anthropogenic impacts aside, CFR diversity has long been explained by a suite of environmental traits that interact on a landscape consisting of a highly heterogeneous array of microregions Latimer et al., 2005Latimer et al., , 2006Linder, 2003). Soil fertility varies from nutrient-poor to mineral-rich and has been hypothesized as a barrier among vegetation types on fine scales (Campbell, 1983;Linder, 2003). Mountain ranges, which create elevation gradients changing from sea-level to over 2200 m within a few kilometers, have been identified as isolating barriers creating potential refugia (Verboom et al., 2009, and plant species richness has been positively correlated with altitude (Cowling & Lombard, 2002).
Overall, CFR plant species diversity declines moving west to east (noted as "Levyn's Law"; Cowling et al., 2017) and has been associated with higher winter-rainfall in the west to less rainfall seasonality in the east (Cowling et al., 1992(Cowling et al., , 1997(Cowling et al., , 1998Cowling & Lombard, 2002).
In teasing apart what biotic and abiotic factors impact CFR diversity, there has been much debate on the evolutionary explanations. Some have suggested the high CFR diversity has been the result of historically high speciation rates, followed by re-occurring extinction and colonization (Linder, 2005), a signature of instability within this heterogeneous landscape. Others have postulated that biome stability and topographic heterogeneity maintain species diversity (Colville et al., 2020;Cowling et al., 2005;Potts et al., 2015).
The endemic fynbos plant Leucadendron salignum, a shrub in the Proteaceae, is a suitable model for evaluating how neutral population processes are influenced by factors in the CFR . Within the Proteaceae-the most prominent flowering fynbos group-the Leucadendron genus is considered basal, with its radiation potentially predating the origin of fynbos vegetation (Barker et al., 2004). L. salignum or "common sunshine conebush" is a uniquely biogeographically diverse CFR species (Barker et al., 2004), found across a range of soil types, rainfall regimes, and elevations, and inhabiting every area of fynbos. It is an evergreen diploid, dioecious plant, whose chloroplast (cp) DNA is maternally inherited (Pharmawati et al., 2004), with cpDNA and nuclear (nu) DNA reflecting seed and pollen dispersal, respectively. It is insect-and wind-pollinated, with weakly developed serotiny (retention of seeds born in cones), and survives fire by resprouting from a lignotuber as well as serotiny for seed dispersal (Hattingh & Giliomee, 1989;Welsford et al., 2014Welsford et al., , 2016). An historically abundant species, L. salignum has seen recent declines in lowland areas like other CFR species, whereas the lack of human interest in agriculture and development in montane fynbos retains large population sizes for the species in these areas (Rebelo et al., 2019).
Here, we sample L. salignum as a model to characterize molecular evolutionary and ecological patterns of CFR population genetic connectivity across spatial and temporal scales. This is the first study to use next-generation DNA sequencing and broad-and fine-scaled landscape genetic analyses to test overall hypotheses of CFR evolutionary stability, as well as effects of human land-usage, in addressing several questions. First, do population genetic patterns mirror species-level patterns of long-term biome stability in the southwestern CFR (which predicts greater and older evolutionary genetic divergence in this region), in contrast to patterns of biome instability in the eastern CFR (which predicts more recent evolutionary divergence in this region)? Second, do CFR population genetic patterns also mirror species-level patterns with evidence of local fine-scale environmental factors, such as seasonal rainfall and elevation, shaping patterns of population genetic connectivity? Lastly, do we find evidence that anthropogenic activity, such as human land-usage, has negatively impacted population genetic connectivity? Addressing these questions in the L. salignum model will help reveal how CFR population analyses shed light on CFR species diversity, as well as provide applied management plans with a means to identify where CFR genetic connectivity is most threatened.

| Study design and sampling
We employed a constructed random sampling of 306 L. salignum individuals from 51 locales (each <0.01 km 2 ) distributed across its biogeographical range (Figure 1; Table S1). Locales were identified from publicly available map data of vegetation types within fynbos habitat (SANBI, 2018), where L. salignum is found (Rebelo et al., 2019). We prioritized the sampling of CFR locales across climate regimes, elevation, vegetation types, and ecozones from dense mountain to coastal fynbos, as well as human land-usage from agriculture to development ( Figure S1).
We generated datasets of chloroplast (cp) and nuclear (nu) DNA to reflect the L. salignum dispersal of diploid embryos as seeds and haploid pollen, respectively. This approach is necessary to determine whether factors on the landscape impact gene flow differentially for seed and pollen dispersal (Holderegger & Giulio, 2010).
For example, phylogenetic analyses of CFR plants that show higher pollen dispersal tend to evolve with lower seed dispersal (Tonnabel et al., 2014). Thus, as previously noted, pollen dispersal in L. salignum is facilitated by wind and insects (Hattingh & Giliomee, 1989;Welsford et al., 2014Welsford et al., , 2016, and thus, it may be the case that while nuDNA (pollen-related) patterns reflect high gene flow, that cpDNA (seed-related) patterns are consistent with lower gene flow. With this multi-tiered approach, the sampling of many locales was our focus for analyses of population connectivity to resolve spatial and environmental heterogeneity (Dyer, 2015;Dyer et al., 2010Dyer et al., , 2012. In this respect, our cpDNA dataset sampled six individuals from each of the 51 locales, which results in 306 cpDNA alleles/sequences (i.e., 306 haploid samples), whereas our nuDNA dataset sampled four individuals from each locale, which results in 408 alleles/sequences (i.e., 204 diploid samples). This genetic sampling is the largest to date of any plant or animal in the CFR and is unusual in comparison to even other land-usage population-based studies outside this region (Miles et al., 2019).

| Genetic marker collection
Sampled leaves were immediately placed in sealable bags of silica gel in the field and subsequently stored at room temperature to rapidly dry them for preservation. Genomic DNA was extracted using the DNeasy Plant Mini Kit (QIAGEN) according to the manufacturer's protocol. We designed a relatively simple and affordable approach that targeted hundreds of putatively neutrally evolving cpDNA and nuDNA fragments using degenerate PCR primers to generate SNP data (Data S1). We used the Nextera XT DNA sample preparation kit (Illumina, Inc.) and Illumina MiSeq (250 bp, paired end) next-generation technologies to barcode and collect full sequence data. The resulting fasta files were aligned using the MUSCLE default parameters in the program Geneious v. 4.8.4 (Kearse et al., 2012), with postsequence filtering steps to enrich our dataset for single-copy haploid and diploid fragments. Details of the data collection are presented in the Data S1.

| Population structure and connectivity analyses
Unless noted otherwise, population genetic summary statistics and analyses were generated using the PopGenome v. 2.7.5 R package (Pfeifer et al., 2014). Estimates of per nucleotide site diversity (θ) were generated for each of the cpDNA and nuDNA datasets and for each locale using the average number of pairwise nucleotide differences (θπ) and the number of polymorphic sites (θs). We estimated Tajima's (1989) D for each of the cpDNA and nuDNA datasets to reflect the SNP frequency spectrum, which can provide insight into deviations from demographic equilibrium (e.g., population expansion). Estimates of genetic diversity between locales were calculated as pairwise F ST values for each of the cpDNA and nuDNA datasets (using the estimator of Hudson et al., 1992). Following Hudson (2000), we used a permutation analysis with an in-house script in R to test for deviations from panmixia (as in . Specifically, this analysis pooled all locales, and randomly sampled allelic variation using our same sample sizes to reconstitute all locales, after which pairwise F ST values were calculated. This simulation was repeated 1000 times and observed F ST values were compared to these simulated distributions.
These simulations determine the extent to which estimates of genetic diversity between locales are sensitive to the sample sizes of the number of individuals and the number of alleles for each locale for both the cpDNA and nuDNA datasets.
We performed a principal component analysis (PCA) on each of the cpDNA and nuDNA datasets to infer genetic clustering of individuals using the gstudio v. 1.5 R package (Dyer, 2016). We also used the Bayesian clustering program fastSTRUCTURE (Raj et al., 2014) with standard defaults to estimate the potential number of genetic clusters (K) for each of the cpDNA and nuDNA datasets. Samples were coded as individuals and run with 20 iterations for each K cluster. The program CLUMPAK (Kopelman et al., 2015) was used to summarize the output, which was visualized using DISTRUCT (Rosenberg, 2004). Both the PCA and fastSTRUCTURE analysis enable interpretation of genetic variation independent of the locales within which it was sampled, that is, locale-labeling is assigned to individuals only after analysis is performed to identify unbiased clustering of genetic variation.
For analyses of genetic connectivity among locales, we employed a graph theoretical framework (Dyer et al., 2010) using the derived statistical measure of conditional genetic distance as cGD.
The statistic cGD is estimated from the genetic covariance among all locales connected through a population network. Thus, cGD reflects the genetic connections among all locales in the network, in contrast to genetic distance estimators such as F ST that make inferences about gene flow from pairwise dissimilarities only (Dyer et al., 2010).
This genetic covariance measure as cGD was used to generate a population graph or "popgraph" for each of the cpDNA and nuDNA datasets, using the popgraph v. 1.5 R package (Dyer, 2017), where nodes represent sampled locales, and edges represent genetic connections among locales (as in Dyer et al., 2012).
Popgraph topologies can be used to depict parameters utilized in "social networks" to define characteristics of the population genetic network (as in Miles, Johnson, et al., 2018). This social network model evaluates genetic relationships among locales and relative contributions of key "actors" using mathematical graph theory (Wasserman & Faust, 1994). This model visually represents gene flow and "hubs" of connectivity among sampled locales. Social network node-specific parameters including closeness, degree, betweenness, and eigenvector centrality were calculated using the popgraph v. 1.5 R package. Higher "closeness" values indicate higher genetic distance between a locale and other locales in the network, whereas lower "degree" values may indicate a locale is more relatively isolated from other locales in the network (Dyer, 2007). Higher "betweenness" values indicate a locale has relatively more paths that flow through it from other locales. "Centrality" reflects the extent to which a locale drives gene flow as a "hub" and reflects a sum of the other popgraph parameters noted above.
Altogether, we interpret lower values of closeness and higher values of degree, betweenness, and centrality as indicative of areas in the popgraph with higher gene flow compared to the rest of the popgraph.
To identify statistical differences and congruence within and between popgraph topologies (e.g., whether patterns of connections and network parameters differ between cpDNA and nuDNA popgraphs), we employed a permutation analysis using the gstudio v. 1.5 R package (as in Dyer et al., 2012;Miles, Johnson, et al., 2018). Specifically, for each popgraph (cpDNA and nuDNA, respectively), each node (representing each locale) was fixed on the popgraph (constructed as noted above), while we randomly permuted the edges that connect the locale nodes. For each social network node-specific parameter that we estimated (e.g., "closeness"), we simulated 1000 randomized popgraphs, each drawn from a distribution that included the observed parameter estimates within that popgraph. For each test, statistical significance was assessed by determining the probability of each observed network parameter value given our simulated popgraphs generated from each node-specific parameter.

| Landscape and genetic association analyses
We used several approaches to interpret the relationship between genetic and landscape variation among our 51 L. salignum locales.
Isolation-by-distance models predict that as geographic distance between populations increases, so does genetic differentiation between populations (Slatkin, 1993), whereas general isolationby-resistance models (McRae, 2006) predict that features on the landscape, including geographic distance, provide "resistance," either by reducing or even facilitating gene flow (Miles et al., 2019).
Thus, through these association analyses, we may interpret a CFR landscape variable as contributing to overall patterns of gene flow in L. salignum, and with further investigation of those patterns, we can then identify where those patterns increase or decrease population connectivity. One assumption in these analyses is that contemporary landscape measures reflect evolutionary and ecological measures, especially for very recent events, such as those involving anthropogenic activities (Leigh et al., 2019;Miles et al., 2019).
Another consideration is that measures of landscape variables also include variance across time (Anderson et al., 2010). In our case here, variables on elevation, precipitation, and vegetation have been longrecognized as evolutionary stable processes (Cowling et al., 2017), and so we may expect these variables have more power to explain relatively older evolutionary genetic events. While human land-usage has been relatively recent, previous analyses also show that it is possible to detect changes in the CFR that reflect contemporary disturbance (Slingsby et al., 2020).
Our first approach used a generalized linear mixed model (GLMM), which accounts for nonindependence of distance matrices (Clarke et al., 2002), and has been shown to be a powerful approach to test multiple hypotheses and investigate nonindependent features in landscape genetic analyses (Row et al., 2017).
We used the function lme of the nlme R package (Pinheiro et al., 2013) and used the Akaike information criteria (AIC) to inform the best models. In this analysis, the pairwise F ST values from the cpDNA and nuDNA datasets were used as our proxies for genetic distances between all locales and constituted our single dependent variable in the GLMM. We employed four landscape features as independent variables ( Figure S1) that are commonly used to model patterns of CFR plant species richness in South African biomes (Thuiller et al., 2006). These four variables are as follows: (i) elevation (ELE), via WorldClim data (value expressed in meters); (ii) annual rainfall (PPT), via WorldClim data (value expressed in millimeters); (iii) seasonal rainfall concentration (CON), ratio of winter to summer seasonal values (see Schulze, 1997); and (iv) vegetation type (VEG), reflecting different vegetation substrates, soil composition, and fertility as categorized by the SANBI (2018) National Vegetation Map study (resolution of 1 km 2 regions over a 1000 × 2000 grid). As we have previously hypothesized, land-usage change may have the potential to fragment population connectivity. Thus, we added a fifth landscape feature as an independent variable in the GLMM to test hypotheses about the potential impact of land-usage (LAN). We used the GeoTerraImage (2017) Southern Africa land-cover dataset (30 × 30 m raster cell resolution), which reflects both natural and human-altered landscape types, the latter of which includes plantations and cultivated agriculture, as well as urban and rural settlement ( Figure S1).
We acquired raster image files for each of the five independent landscape variables, as noted above, and cropped them to fit the geographic area of the 51 sampled locales (Figures 1 and S1).
Because spatial and landscape resistance distances are inherently correlated, for example, as one increases in elevation between two locales, there is also an increase in geographic distance, we included geographic distance as a covariate within the GLMM to assess the relative independent contributions of each landscape variable.
Geographic distances (GEO) between locales were estimated from latitude and longitude coordinates as Euclidean distances using the PASSaGE program (Rosenberg & Anderson, 2011). For each of the five landscape variables, the mean and variance of the raster pixel values among all pairwise connections of the 51 locales were each calculated using the raster R package (Hijmans & van Etten, 2012). The mean and variance calculated between all locales reflect different spatial properties. For example, the mean of an increase, followed by a decrease, and again an increase in elevation between two locales could be identical to a mean score that reflects "no change" in elevation between two locales. In this way, we can capture this variation in change between two points, which is what potentially influences gene flow between two points on the landscape. In addition, calculating landscape resistance distance as the difference between two locale points misses critical information.
That is, elevation at locale 1 and elevation at locale 2 is identical and thus the difference between these two locale points is "0," but dramatic changes in elevation occur between those two points, as we observe over even short geographic distances in the CFR. These mean and variance statistics were used as landscape resistance measures in separate GLMM analyses of each of the cpDNA and nuDNA datasets to test for relationships between landscape and genetic variables.
Our second approach used a permutation analysis with the gstudio v. 1.5 R package (as above) to identify significant relationships between genetic connectivity (estimated as cGD), and our landscape variables, as in previous analyses (Dyer et al., 2012;Miles, Johnson, et al., 2018). Specifically, for each of the cpDNA and nuDNA datasets, each of the 51 sample locales (or nodes) was fixed on the landscape of the popgraph, which was estimated using cGD, or the genetic covariance among all locales (as above). As above, nodes of the popgraph reflect the locales and connections reflect the cGD measures or edges of the popgraph. We overlaid the connections of the popgraph on each of the raster maps of each of the five landscape variables (ELE, PPT, CON, VEG, LAN) and then calculated the observed mean and variance resistance distances of all the pairwise connections that separate the locales. We then simulated 1000 randomized popgraphs that were constrained by edge density but also by the edge distribution on individual nodes, generating a null distribution of popgraphs. These null distributions for each landscape feature generated from simulations (using the observed popgraph data) were then compared to the observed popgraph for each landscape variable to look for differences in congruence. Statistical significance was assessed by determining the probability of the observed values given our simulated popgraphs of each landscape variable. This analysis randomizes all the genetic, landscape, and spatial variables observed among all locales to simply ask whether specific landscape features are more often associated with genetic distances between locales (i.e., genetic connectivity or gene flow), than expected by chance.

| Estimates of molecular evolutionary relationships
We conducted several analyses to place genetic divergence among locales and lineages in an evolutionary framework. While we are interested in the phylogenetic structure among locales and lineages, we are also interested in generating date estimates for this genetic structure. Although prior phylogenetic tree analyses of the Leucadendron genus have been published (Barker et al., 2004;Hoffmann et al., 2015;Tonnabel et al., 2014), molecular clock divergence estimates have not been previously estimated for this group. Thus, we generated molecular clock divergence estimates for Leucadendron using previously published species nucleotide data (Table S2) and the approach of Valente et al. (2010), which estimated divergence times of the closely related South African Protea genus using a relaxed Bayesian MCMC approach implemented in BEAST (Drummond & Rambaut, 2007). We also collected nucleotide sequence data for a closely related outgroup

| RE SULTS
After multiple levels of postdata filtering, the resulting dataset from our sampling of 51 L. salignum locales included >100 SNPs from 306 haploid cpDNA sequences and >800 SNPs from 408 nuDNA sequences. Statistics reflecting SNP frequency spectra for marker datasets (cpDNA, D = −0.71; nuDNA, D = −1.60) are not significant and show no deviation from demographic equilibrium. Nucleotide sequence diversity for the overall dataset is higher for nuDNA than it is for cpDNA; however, estimates within locales are considerably lower, especially for cpDNA (Table S3), suggesting significant genetic variation is distributed between locales. This is indeed the case as we find significant genetic divergence among many locales ( Figure   S2), and much higher overall structure in cpDNA (F ST = 0.81) compared to nuDNA (F ST = 0.10).
Several key observations emerge from the BEAST analysis of cpDNA (Figure 2), which support evolutionary groups that mirror geographic clustering (Figure 1). We find support for four groups in the western CFR (labeled here as "north," "west," "south," "central"), whereas there is support for a single group in the eastern CFR (labeled as "east"). The estimated coalescence for the population sample at 1.17 Mya (95% CI = 0.73-1.63 Mya) is followed by two highly supported groups, one that includes the west and central groups, and one that include the north, south, and east groups. The west, central, and south groups each have estimated coalescent times coinciding at ~0.50 Mya (0.26-0.78), with the east group emerging more recently from an apparent common ancestor in the southwestern CFR at ~0.31 Mya (0.17-0.47), and the smaller north group emerging more recently at ~0.17 Mya (0.05-0.30). The overall estimate of coalescence of the nuDNA lineages at ~1.7 Mya was consistent with that of the cpDNA tree.
However, unlike the cpDNA tree, there was no significant internal node support for any individual nuDNA groups. To independently confirm these BEAST results, we performed a simple model-free neighbor-joining (NJ) tree analysis using MEGA (Kumar et al., 2016), which supported the BEAST topologies of the cpDNA lineages as well as the nuDNA polytomy.
The fastSTRUCTURE analyses mirror the PCA results. The cpDNA analysis for K = 2-6 (additional K runs resulted in no further clustering) reflects multiple clusters moving geographically west to east (Figures 1 and 4). The only exception is the south group (K = 6), which includes individuals that are geographically distantly separated, but have in common their southwestern coastal locations. Notably, the central group additionally separates into multiple geographic clusters (K = 5), and coincidentally, the cpDNA PCA identified this more subtle genetic separation within the central region (PC5-6 = 8.2%).
On the other hand, the nuDNA fastSTRUCTURE analysis found no genetic clustering (i.e., could not reject a model of K = 1, Figure S3).
In our analysis of the network parameters from the cpDNA popgraph, we found that connections among eastern CFR locales have significantly higher connectivity (i.e., lower "centrality," p = 0.0002) and are more genetically similar (i.e., lower "closeness," p = 0.00009) in contrast to connections anywhere else in the popgraph ( Figure 5; Table S4). The measure of "betweenness" reflects nodes that significantly influence connectivity across the network as a whole. In the cpDNA popgraph, the locales with the highest "betweenness" values are those found at each of higher elevation, within the Cape Peninsula, and also along the southwestern coast (Table S4 and Figure S4). When locales are disconnected from the popgraph network it does not reflect "isolation" but instead signifies that they do not contribute added information to patterns of genetic connectivity among the other locales. This was the case for only two locales in the cpDNA popgraph ( Figure 5). In contrast, the locales of the nuDNA popgraph are almost entirely disconnected; that is, locales equally share nuDNA allelic variation, indicating no distinct pattern of gene flow ( Figure S5). This panmictic population genetic pattern for nuDNA is consistent with the BEAST, PCA, and fastSTRUCTURE results.
The landscape genetic analysis performed in the GLMM shows that geographic distance explains a significant proportion of genetic differentiation (measured as F ST ) in both cpDNA and nuDNA datasets (9% and 15%, respectively), but that landscape variables explain very little variance (Table S5). In contrast, our simulation analyses that assessed the relationship between landscape features and genetic connectivity (measured as cGD) provided significant results ( Figure S6). These results are consistent with the cGD measure as a statistically more powerful proxy of gene flow as it reflects shared variation through the entire network, whereas F ST reflects dissimilarity, and not explicitly connectivity among locales (Dyer et al., 2010). As noted previously, the nuDNA data show no patterns of genetic structure or connectivity; whereas, landscape permutation analyses with the cpDNA popgraph reveal that each of ELE, PPT, CON, and LAN variables were significant ( Figure S6). These analyses find that observed genetic connections on the popgraph occur significantly more often in areas with (i) lower elevations, (ii) higher annual precipitation, and (iii) higher winter-rainfall concentration than predicted by the simulations.
Further investigation into the significant LAN results found that compared to the simulations ( Figure S6), genetic connections on the popgraph were significantly underrepresented in areas where "fynbos thicket" (i.e., a mosaic of Cape fynbos and subtropical F I G U R E 1 Map of southern South Africa with 51 sampled locales of 306 Leucadendron salignum (Table S1 for locale information). Pie charts provide geographic reference to color coding in phylogenetic (Figure 2), population structure (Figures 3 and 4), and population connectivity ( Figure 5) analyses. Shaded area denotes recognized Cape Floristic Region (after Cowling et al., 2017) F I G U R E 2 Consensus BEAST phylogenetic tree topology of 306 Leucadendron salignum chloroplast DNA SNP haplotypes, rooted by outgroup, with estimates of genetic divergence times (see Materials and Methods). Black dots on nodes represent major branches with support >75% F I G U R E 3 Bi-plots of the first four principal components of Leucadendron salignum datasets using (a) chloroplast DNA SNPs from 306 sequences, and (b) nuclear DNA SNPs from 408 sequences. Color-coding reflects the five major phylogenetic tree clusters (Figure 2) F I G U R E 4 fastSTRUCTURE analysis of 306 Leucadendron salignum chloroplast DNA SNP haplotypes from 51 sampled locales (see Figure  1 for geographic sampling, Table S1 for locale information) thicket, Vlok et al., 2003) and "commercial cultivation" (i.e., commercial agriculture) occur ( Figure S1). Thus, these two land-usage types appear to significantly reduce genetic connectivity among locales where these features exist. As previously noted, a GLMM is a valuable approach as it models the potential nonindependence of landscape features; nonetheless, the four landscape variables that showed significant relationships with genetic connectivity are

| D ISCUSS I ON
Using the CFR-wide L. salignum plant as a population and landscape genetic model, we set out to test hypotheses of the evolutionary and environmental factors that drive the rich CFR plant species diversity. Our analyses find that while pollen-related genetic signatures reflect a panmictic population with gene flow largely uninhibited, seed-related genetic signatures reflect high population genetic structure, and demographic patterns of relatively deeper and older genealogies in the southwestern CFR compared to lineages in the northern and eastern CFR. Although we find several environmental factors associated with population genetic connectivity, a striking result was the association of human land-usage with reduced gene flow. We discuss how these patterns are consistent with evolutionary hypotheses, the implications for management strategies, and future directions to evaluate CFR genetic models.

| Contrasting dispersal patterns of genetic differentiation
Our results from nuclear and chloroplast genetic markers are a sharp reminder of how sampling only one genetic mode of inheritance could potentially provide an incomplete picture. Although L. salignum nuDNA shows higher genetic diversity compared to cpDNA (as is expected for biparental, recombining nuDNA markers; Wolfe et al., 1987), in contrast to cpDNA, nuDNA shows very little population structure, genetic connectivity, or phylogeographic resolution. Patterns of cpDNA variation reflect the hypothesis of historically consistent CFR edaphic and climatic conditions (Cowling et al., 2017;Potts et al., 2015), which likely selects for lower seed dispersal in habitat-specialist taxa (Lopez et al., 2008;Olivieri et al., 2015). In contrast, patterns of pollen dispersal in L. salignum unsurprisingly reflect pandemism, consistent with phylogenetic analyses of CFR plants that show higher pollen dispersal tends to evolve with lower seed dispersal (Tonnabel et al., 2014). While pollen dispersal in L. salignum is facilitated by wind and insects (Hattingh & Giliomee, 1989;Welsford et al., 2014Welsford et al., , 2016, the lack of pollen geographic variation here does not suggest that pollinator movement in the CFR is also unaffected by environmental factors. Given the different pollen and seed patterns here, population genetic analyses of CFR pollinators could elucidate the co-evolutionary relationship they have with plant diversity (Eimanifar et al., 2018).

| Signatures of broad-and fine-scale patterns of CFR genetic diversity
Population genetic analyses for L. salignum reveal contrasting patterns of CFR diversity across different spatial and temporal scales.
Population genetic structure overall is significantly high, yet much of this variation is consistently found between southwestern CFR locales in very close geographic proximity. In fact, geographic distance explains a significant, yet small, overall proportion of population genetic structure on broad-scales (9% of cpDNA, and 15% of nuDNA diversity, Table S5). As previously noted (e.g., Cowling et al., 2017;Forest et al., 2018), and predicted here, phylogeographic evolutionary analyses show lineages from the southwestern CFR exhibit higher genetic structure and older estimated ages than lineages found across the most distantly geographically separated locales in the CFR. That is, a simple isolation-by-distance model cannot explain patterns of genetic divergence observed here for either marker type.
Interestingly, although eastern lineages appear to have a derived ancestry from the west, their estimated ages are not necessarily recent. That is, eastern lineages appear to have a single evolutionary root and very low population structure, but our coalescent analyses find an estimated age of ~500 Kya, and no evidence of reduced nucleotide diversity, inconsistent with a very recent colonization event (Linder & Hardy, 2004). These overall geographic patterns are consistent with higher phylogenetic diversity, but older divergence F I G U R E 5 Popgraph of genetic connectivity among 51 sampled locales of Leucadendron salignum chloroplast DNA (see Figure 1 for geographic sampling, Table S1 for locale information). Color-coding reflects the five major phylogenetic tree clusters (Figure 2) of Cape clades in the west, compared to lower phylogenetic diversity and more recent species radiation in the east (Richardson et al., 2001;Schnitzler et al., 2011;Valente et al., 2010;Verboom et al., 2009).
While the few previous CFR plant phylogeographic studies were limited in sample size and scope, they revealed consistently high population structure (Lexer et al., 2013(Lexer et al., , 2014. Here, population genetic signatures are consistent with the hypothesis that southwestern CFR high environmental and evolutionary diversity is the result of Pleistocene climate stability, whereas eastern CFR diversity is the result of climate instability during the same period Huntley et al., 2016;Linder, 2008;Verboom et al., 2009)essentially reflecting "two CFRs" (Cowling et al., 2017). This theory predicts that biome stability in the southwestern CFR would foster a deep, evolutionarily stable L. salignum genealogical diversity, but with reduced gene flow and continuous lineage divergence, without extinction, over even short geographic distances due to a highly heterogeneous landscape. While biome diversity in the east is exceptionally rich , its instability has evolved a L. salignum genetic diversity with signs of constant bottlenecks (which erode genealogical diversity), and evidence of invasion from older southwestern lineages. This latter explanation comes from the observation in our phylogenetic analysis that eastern lineages have a single common ancestor (derived from southwestern lineages) with reduced phylogenetic diversity. The southwestern CFR may act as an evolutionary genetic "source," typified by genetic structure that is higher and older, that feeds into the eastern and northern CFR, which may appear as evolutionary genetic "sinks" typified by re-occurring lineage extinction and colonization from the southwest.

| Genetic connectivity in the CFR is shaped by biotic and abiotic factors
Landscape genetic results show that genetic connectivity among geographically distant locales is significantly associated with lower elevation and higher winter-rainfall. This overall pattern mirrors previous studies showing declines in species richness moving west to east that are consistently associated with lower winter-rainfall (Cowling et al., 2017). By marrying our patterns of gene flow with patterns of topographic heterogeneity (Figures 1 and 5; Table S4 and Figure S4), we can identify geographic areas with conservation implications. Our social network analyses find the greatest drivers of gene flow are geographically distant locales in a stretch of lowland elevation along the southwestern coast, for example, Cape Peninsula, across False Bay, and the Cape Agulhas Plain, into the eastern CFR. Interestingly, popgraph results show that this southern coastal route also maintains connectivity among clusters of southwestern locales separated along low elevation routes moving in two directions. For example (Figures 1 and S4), clusters are separated by the Cederberg, Boland, and Hottentot and Holland mountain ranges that run north to south, with smaller pockets separated by the Hex River and Riviersonderend ranges. Running west to east we see clusters separated by the Cape Fold Mountains into the Langeberg and Swartberg ranges. Although locales in these ranges are distantly separated, they maintain high genetic connectivity running west to east linked by lowland elevation, whereas moving from coastal to montane locales, we see low genetic connectivity as a result of the parallel mountain ranges acting as barriers. These high-resolution genetic patterns strongly coincide with species richness patterns predicted by large studies of biome stability and topographic heterogeneity at the global scale as well as within the CFR (Colville et al., 2020).
Certain CFR environmental features, such as elevation and precipitation, may be biogeographically correlated, along with vegetation type; thus, our statistical model was designed to account for these associations to identify factors that may play independent roles in shaping population structure. For example, the observation that fynbos thicket is associated with reduced gene flow in the east is also coincident with seasonal winter-rainfall showing a similar pattern; however, these variables are not always naturally associated (Becker et al., 2015;Potts et al., 2013). The subtropical thicket of southern Africa was hypothesized to have undergone great declines during colder glacial periods of the Pleistocene Vlok et al., 2003), which as noted previously, may have resulted in rapid species turnover and extreme bottlenecks in the east. This consistent upheaval may have had an effect on the evolutionary dynamics of population connectivity of flora associated with subtropical thicket (Potts et al., 2013). In fact, our exhaustive search identified only a few eastern pockets of L. salignum (Figure 1) nestled within dense thicket, and indeed, these eastern L. salignum lineages exhibit more recent, late Pleistocene coalescence. Our population genetic analyses suggest that the consistent evolutionary turnover of eastern CFR types, such as subtropical thicket, may have potentially both depressed invasion of and increased fragmentation of eastern fynbos lineages, such as L. salignum.

| The impact of human activity on CFR diversity
Our project was initiated by the need to document how anthropogenic activity impacts biodiversity hotspots, like the CFR, through the use of evolutionary genetics. Thus, the observed genetic signature of negative anthropogenic impact is of interest. Environmental factors, such as elevation and precipitation, appear to have long-impacted L. salignum genetic diversity, dating to at least the Pleistocene; however, human land-usage is expected to impact the CFR in only the last ~300 years (Gaertner et al., 2016;Rebelo et al., 2011). Recent studies show that anthropogenic effects on allele frequencies and gene flow may take a considerable number of generations, depending on effective population sizes and standing genetic diversity, unless these effects are severe (Johnson & Munshi-South, 2017;Leigh et al., 2019;Miles et al., 2019). Thus, our result with human landusage underlies the severe impact that anthropogenic activity may have on gene flow, especially in coastal and agricultural areas of the southwestern CFR, where development has exponentially expanded (Horn & Van Eeden, 2018). This development into coastal lowlands includes not only areas most vital to maintaining gene flow among the most geographically separated CFR locales, but also where we identify some of the most ancestral, oldest, and highest levels of genetic structure in the CFR overall.
This anthropogenic signature is of great concern given that 75% of South Africa's human population is expected to live in urban areas by 2030 (Cooperative Governance & Traditional Affairs, 2016). As anthropogenic activity fragments the landscape and reduces gene flow, we artificially shift the adaptive landscape (Wright, 1982) and select for very different life-history traits that inhabit narrow geographic areas. Urban evolutionary studies have shown that in selecting for potential urban "ecotypes" (Johnson & Munshi-South, 2017;Miles, Johnson, et al., 2018;Rost et al., 2009;Winchell et al., 2016), we stand to bottleneck vital CFR genotypic and phenotypic variation that has accumulated over millions of years of evolution. This study reminds us that while there is a need to focus on rapidly growing urban areas, we must be cognizant of other land-usage such as agriculture as it also erodes vital biodiversity.

| Future directions and applications
In identifying both evolutionary and environmental features associated with population connectivity in the CFR, several future directions emerge. First, L. salignum is an appropriate biogeographic and neutral model to generate hypotheses, but we need replication from other biogeographically diverse taxa for comparison, and not for plants alone. Lexer et al., (2013) noted that the heterogeneous and complex mosaic landscape is similar for both plants and animals alike, and examining diverse taxa can identify how these drivers are similar across phylogenetic time, mating systems, and genomes. The limited number of animal studies to date in the CFR shows far less phylogeographic diversity and evidence of range expansion consistent with the impact of Pleistocene climatic changes (Daniels et al., 2007;Portik et al., 2011;Tolley et al., 2006). Yet, this result is taken with caution given the limited sampling design, as well as the fact that some studies focus on animals of conservation interest but with current ranges and population sizes that deviate greatly from historical ones (Smit et al., 2007;Van Hooft et al., 2002). Thus, using larger and more sampling locales to disentangle the complex demographic histories of these taxa will be hard-pressed to reveal the underlying drivers of CFR diversity and species divergence. In this respect, studies of unthreatened species with low conservation priority, but which have broad distributions such as the current study, have great conservation implications as they can help build neutral evolutionary models with which to test patterns from species of management concern.
A related challenge is teasing apart the relationship between geographic range and population genetic diversity, especially for neutral traits (McDonald et al., 1996;Morueta-Holme et al., 2013;Slingsby & Verboom, 2006). As environmental history is an important factor in driving CFR evolutionary change, then sampling southwestern and eastern species with similar geographic ranges may show similarly high diversity/older lineages and low diversity/recent lineages, respectively. In addition, comparative phylogeographic analyses among Mediterranean biomes with similarly high diversity (Feliner, 2014) can reveal drivers unique to the CFR. That is, we may predict that signatures of genetic diversity and phylogeographic structure are similar across these biomes, but that landscape genetic analyses identify unique drivers and thus unique conservation priorities. Even within the CFR, we find that no one landscape feature alone can explain the phylogeographic diversity. Lower elevation and vegetation types alter connectivity in different ways, and thus, an intense sampling along specific mountain ranges, substrates, and watersheds across the CFR is needed to separate their influence on fine scales. One example here is from the association that "thicket" has with lower L. salignum gene flow in the eastern CFR. Potts (2017) conducted phylogeographic analyses of Albany subtropical thicket on fine scales and showed that watersheds have been long-term barriers to gene flow creating isolated catchments that reflect evolutionary conservation units. By fine-scale sampling within the CFR, we may reveal that these catchments isolate L. salignum populations as well. These similar patterns across species lend support to these local ecosystems needing consideration as evolutionary conservation units.
While our results shed further light on long-term environmental features to aid in conservation priorities (such as lowland and coastal fynbos), other CFR species may reveal other priorities. For example, our study supports the need of fine-scaled population network analyses that target contemporary factors of high threat such as human impact. We need a multi-tiered approach to understanding CFR biome evolution  that includes CFR population genetic studies to create a narrative of the consistent eco-evolutionary factors driving long-and short-term diversity. As the affordability of large-scale genetic data collection continues, our approach can be applied to any CFR species to help guide conservation plans.

ACK N OWLED G EM ENTS
We thank Jan Vlok for detailed information on locale and plant iden-

CO N FLI C T O F I NTE R E S T
The authors declare no competing interests with research described in this paper.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data for this study are available at the Dryad Digital Repository: https://doi.org/10.5061/dryad.612jm 642glink.