Genetic and ecological consequences of recent habitat fragmentation in a narrow endemic plant species within an urban context

Understanding the timescales that shape spatial genetic structure is pivotal to ascertain the impact of habitat fragmentation on the genetic diversity and reproductive viability of long-lived plant populations. Combining genetic and ecological information with current and past fragmentation conditions allows the identification of the main drivers important in shaping population structure and declines in reproduction, which is crucial for informing conservation strategies. Using historic aerial photographs, we defined the past fragmentation conditions for the shrub Conospermum undulatum, a species now completely embedded in an urban area. We explored the impact of current and past conditions on its genetic layout and assessed the effects of genetic and environmental factors on its reproduction. The historically high structural connectivity was evident in the genetics of the species. Despite the current intense fragmentation, we found similar levels of genetic diversity across populations and a weak spatial genetic structure. Historical connectivity was negatively associated with genetic differentiation among populations and positively related to within-population genetic diversity. Variation partitioning of reproductive performance explained ~ 66% of the variance, showing significant influences for genetic (9%), environmental (15%), and combined (42%) fractions. Our study highlights the importance of considering the historical habitat dynamics when investigating fragmentation consequences in long-lived plants. A detailed characterization of fragmentation from 1953 has shown how low levels of genetic fixation are due to extensive gene flow through the non-fragmented landscape. Moreover, knowledge of the relationships between genetic and environmental variation and reproduction can help to implement effective conservation strategies, particularly in highly dynamic landscapes.


Introduction
Flowering plants are a crucial component of most terrestrial ecosystems (Ollerton et al. 2006). They provide food and shelter for vertebrate and invertebrate populations, which in turn are important for plant pollination and seed dispersal. Yet, as a consequence of human modification of landscapes, the extinction risk of many species is increasing worldwide (Newbolt et al. 2015). Land use changes are predicted to drive the most significant effects on biodiversity throughout this century (Sala et al. 2000) and among human-built environments, urban areas are the fastest growing (Hobbs et al. 2006).
The structural changes imposed by anthropogenic habitat fragmentation can have a range of potentially deleterious genetic and demographic consequences (Haddad et al. 2015). For plants, such consequences may result from the loss of pollinators and other alterations to plant-pollinator interactions, as well as from the reduction of reproductive output due to drastic changes in the size and isolation of populations (Aguilar et al. 2019;Lienert 2004;Potts et al. 2010). Understanding the consequences of habitat fragmentation and the main drivers of declines in plant reproductive output is therefore critical in implementing effective restoration strategies.
Population genetics theory predicts that small and highly fragmentated populations will lose genetic variation due to genetic drift, leading to highly structured populations. Although some gene flow between remnant populations may be enough to maintain genetic diversity lost by genetic drift in small populations, increasing level of fragmentation is likely to eventually disrupt pollen and seed dispersal vectors from exchanging migrants between remnants. Under such conditions, extinction risk for plant species is expected to increase due to genetic erosion, expression of inbreeding depression (Charlesworth and Willis 2009), and/or a loss of mating types (e.g. self-incompatibility alleles; Pickup and Young 2008). Empirical studies and meta-analyses, however, have found mixed support for these theoretical predictions (e.g. Aguilar et al. 2008;Kramer et al. 2008;Lowe et al. 2015) with often idiosyncratic species response to habitat fragmentation, highlighting the complexity of such systems. Indeed, for many plant species, attempts to identify the effects of fragmentation on population genetics measures over ecological time scales are often confounded by life history and reproductive traits that buffer populations from expected loss of genetic variation and reproductive output (Vranckx et al. 2012;Broadhurst et al. 2017;Miles et al. 2019). These include overlapping generations, long generation times, seed banks, and diverse mating system and pollination dispersal mechanisms.
In particular, adult cohorts of long-lived plant species may not show the expected population genetics consequences of fragmentation simply because insufficient time has elapsed since land use change for such signals to be detectable (Vranckx et al. 2012;Lowe et al. 2015). Accounting for the interaction between the time span of the fragmentation process and the life cycle of the focal species, therefore, appears to be particularly important when studying anthropogenic fragmentation effects on long-lived plants (Kramer et al. 2008;Lowe et al. 2015). Combining study of genetic information with current and past environmental characteristics of fragmented populations of plant species provides an integrated approach to understanding how present-day and historical fragmentation conditions may correlate with population genetics measures (Bacles and Jump 2011). This approach, however, requires systems where quantitative data on fragmentation can be placed into historical context together with genetic diversity. Nonetheless, a detailed characterisation of ecological indices (i.e. fragment size, isolation) through time is not always feasible since fragmentation can be a complex process over 1 3 time, with insufficient historical records. As a consequence, relatively few studies have incorporated spatially explicit historical isolation measures and characteristics of the surrounding landscape as possible factors in explaining the genetic diversity of fragmented populations (e.g. Da Silva Carvalho et al. 2015;. Preventing the loss of genetic diversity is widely considered vital for maintaining the adaptive potential of populations in the face of environmental change (Booy et al. 2000;Reed and Frankham 2003). However, we also need to understand whether current reproductive output in fragmented landscapes correlates with components of genetic variation and attributes of the fragmented populations. It has become increasingly accepted that population genetics and reproductive ecology are complementary approaches to address central questions about the conservation of plant biodiversity (Lowe et al. 2009: Real 2017. Yet, interactions between the biotic and abiotic environment, genetics, and reproductive output of endangered flora remains poorly understood. It is known that genotype-environment interactions are crucial in explaining phenotypic variation, including reproductive traits (Via and Lande 1985). Therefore, looking at interactions between genetic and ecological variation and reproductive output is especially important in fragmented habitats where the ecological and environmental condition varies greatly and may be confounded.
Conospermum undulatum is a threatened insect-pollinated long-lived lignotuberous shrub with a narrow natural range embedded in the rapidly expanding metropolitan area of the city of Perth (Department of Environment and Conservation 2009). The area has been impacted by fragmentation because of relatively recent urban development (Wardell-Johnson et al. 2016). This species represents an ideal candidate to study the genetic consequences of recent fragmentation by urbanization on plants as it is feasible to quantify pre-fragmentation environmental characteristics across its entire distribution range. In this system, historic aerial photographs provide valuable information to distinguish historical from contemporary effects of environmental conditions on the present-day genetic diversity. Delnevo et al. (2019a, b) recently showed that habitat fragmentation negatively influenced the reproductive output of C. undulatum in terms of seed production and germination rate. Further investigations showed that this was due to both pollen quantity and, particularly in small populations, quality limitation (Delnevo et al. 2020a, b). However, it remains unclear to what extent the genetic diversity of this species has been affected by population size reductions and fragmentation or whether the genetic consequences of fragmentation are playing a role in changing the reproductive dynamics of the species. Moreover, the historical data enables a test of whether the genetic consequences of habitat fragmentation pre-date recent population declines or reflects deeper time-scales of population distributions.
To address these gaps, we performed a genetic characterisation of the remnant populations of C. undulatum and investigated how variation in genetic diversity is related to contemporary and historical environmental factors and to its current reproductive performance. The naturally restricted distribution of this species provides the opportunity to conduct a genetic survey at the whole species level to investigate its recent response to a large reduction in habitat. We quantified the extant genetic diversity of these populations embedded within a large urban area. This was used to address three questions: (1) Do recent decreases in population size and increased isolation explain current population genetic structure? (2) To what extent is the spatial distribution of genetic diversity associated with contemporary versus historical fragmentation and environmental factors? (3) Are genetic diversity and environmental factors correlated with reproductive output of this threatened species? 1 3

Study area
The study was conducted in southwest Western Australia within the Swan Coastal Plain bioregion (Fig. 1), a low-lying coastal plain that is part of the southwest Australia global biodiversity hotspot (Mittermeier et al. 2004). This region has experienced extensive land clearing for urbanisation, which is centred around the capital city of Perth. Approximately 35% of natural vegetation remains on the Swan Coastal Plain, although only the 10% occurs in protected areas (Wardell-Johnson et al. 2016). Our target species, C. undulatum, Fig. 1 a Species distribution model for Conospermum undulatum (SDM; darker grey) superimposed on a shaded relief map with 10 m contour lines (http: //services.arcgisonline.com/ArcGIS/rest/services/Elevation/World_Hillshade/MapServer/tile/{z}/{y}/{x}). Darker green represents historic (1953) extent of native bushland (native bushland growing over incompatible environment unlikely to have contained C. undulatum). Bright green represents current extent of native bushland containing C. undulatum. Coordinates cannot be provided for specific locations of Threatened Flora. Pie charts shows the assignment probability of each population to genetic clusters inferred at K = 4 in Bayesian clustering. b Genetic ancestry of 293 individuals sampled from 14 populations, estimated using STRU CTU RE analysis of microsatellite markers. Populations are ordered based on genetic assignment from southern to northern, eastern, and western clusters grows on rapidly drained soils on the eastern part of the Swan Coastal Plain adjacent to the Darling Scarp. The entire distribution range of the species is characterised by highly leached sandy soils of the Bassendean, the Pinjarra, and the Forrestfield systems ( Fig. 1).

Biology of study species
Although occupying a narrow range, Conospermum undulatum (Proteaceae) is fairly widespread in remnant patches of native bushland within this range where compatible sandy soils occur. It is a long-lived non-clonal species that can survive at least several decades as individual shrubs can resprout from semi-subterranean woody lignotubers after disturbance such as fire or mechanical damage from strong winds or trampling by animals (Bennet 1995). It is an entomophilous self-incompatible species, although it is unclear whether this is predominantly a genetic-or physiological-based incompatibility. The flowers possess an active pollination mechanism. Insect visitors apply pressure with their mouthparts at the base of the calyx which triggers the style to flick away from the fertile anthers and strike the insect. This characteristic pollination system is a physical barrier to autogamous selfing and is activated by a restricted group of native insect pollinators, with the small specialised native bee Leioproctus conospermi being the most effective pollinator, follwed by ants (Delnevo et al. 2020a, b). Highly mobile visitors, such as the introduced European honeybee, are not able to pollinate the small and specialised flowers of C. undulatum (Delnevo et al. 2020a, b). Other generalist pollinators, such as dipterans, have been recorded visiting Conospermum flowers, but they are mostly unable to trigger the mechanism, resulting in extremely small to no pollen load (Delnevo et al. 2020a, b).

Sample collection and genotyping
In 2018, we collected leaf samples from 24 individuals of C. undulatum from each of 14 of the 17 known populations (Fig. 1). Of the three unsampled populations, two were inaccessible, being located on private land, and one only contained three individuals, making the sampling potentially detrimental to the population and the estimation of genetic indices unreliable. Census population sizes ranged from eight to approximately 880 plants. In populations with less than 24 plants, we collected leaf samples from all individuals. Populations with > 24 plants were sampled using a stratified random approach, where sampling reflected the pattern of distribution of plants within a population and covered the full geographic extent of the population. We ensured each sampled plant was at least 10 m apart where possible.
A total of 293 samples were genotyped for the full set of 19 microsatellite loci developed for Conospermum undulatum; details of DNA extraction, polymerase chain reaction conditions and scoring of electropherogram profiles are provided in Delnevo et al. (2019a, b). The final genetic dataset was assembled after discarding, as a precaution, 5 loci (Cu15, Cu17, Cu29, Cu41 and Cu45) that showed the possible presence of null alleles in the majority of populations at average frequencies ranging from 0.167 to 0.402.

Genetic diversity and population structure
To detect the possible presence of null alleles and genotyping failures in each population, we used the Bayesian approach implemented in INEST v2.2 (Chybicki and Burczyk 2009).

3
Each model was run for 5 × 10 5 cycles discarding the first 5 × 10 4 iterations and keeping results every 100 updates. After model selection based on the BIC criterion, we discarded loci for which null alleles were detected at > 50% populations. Standard genetic diversity parameters (number of alleles, A; effective number of alleles, Ae; observed and expected heterozygosity, H O and H E , allelic richness, Ar; and private allelic richness, PAr) were then calculated using GenAlEx 6.5 (Peakall and Smouse 2012) and hp-rare (Kalinowski 2005). Deviations of genotypic distributions from Hardy-Weinberg expectations (HWE) and genotypic linkage disequilibrium (LD) between all pairs of loci were assessed using GENEPOP 4.0 (Rousset 2008). A measure of the contribution of each population to the total allelic diversity (A T ) and its within-(A S ) and between-population (D A ) components was estimated with metapop2 (López-Cortegano et al. 2019). The same program was used to estimate the percentage of individuals from each population contributing to a virtual pool of 1000 individuals created to maximise the total average number of alleles (maxK).
Parameters of genetic fixation (G ST ;Nei 1977) and differentiation (Jost's D; Jost 2008) among populations were estimated with GenAlEx 6.5, and their statistical significance was assessed with 999 permutations. The same software was used for calculating pairwise population genetic distance matrices used in both Principal Coordinate Analysis (PCoA) and Mantel test.
The Bayesian clustering algorithm implemented in STRU CTU RE v. 2.3 (Pritchard et al. 2000) was used to assess the putative number of differentiated genetic clusters. The mostlikely number of clusters (K) in which individuals can be divided was assessed with two approaches: using the empirical statistic ΔK calculated by STRU CTU RE HARVESTER (Evanno et al. 2005;Earl and von Holdt 2012), and by estimating the log likelihood for each K (Ln P(D) = L(K)). The default parameter settings of the admixture model with correlated allele frequencies and no LOCPRIOR were used, with K varying from one to 15. Each run consisted of 1 × 10 5 burn-in iterations and 5 × 10 5 data collection iterations, and ten runs for each K were performed. Different runs for the same K were averaged using the software CLUMPP (Jakobsson and Rosenberg 2007).

Species distribution model and spatial analysis of genetic differentiation
The likely pre-fragmentation distribution of C. undulatum was reconstructed based on three environmental variables. These variables encompassed topographic (elevation and slope) and soil type (Appendix S4). Given the limited extent of the study area, climatic variables such as temperature and precipitations were not included in the model as they were equal throughout the area and thus uninformative. Elevation and slope were obtained from 10 m contour lines using QGIS 3.18.2-Zürich (https:// qgis. org/ en/ site). Soil types were obtained from soil-landscape mapping covering Western Australia at a scale of 1:20,000 (https:// catal ogue. data. wa. gov. au/ datas et/ soil-lands cape-mappi ng-best-avail able) and rasterised with QGIS. All layers were imported in R software and adjusted to the same coordinate projection system (EPSG: 4326 or + proj = longlat + ellps = WGS84 + datum = WGS84 + no_defs) and cropped to an extent delimited by longitude: 115.876, 116.127, and latitude: -32.079, -31.873.
The statistical machine-learning model based on the principle of maximum entropy, Max-Ent (Phillips et al. 2006;Sánchez-Ramírez et al 2018), was used to model habitat suitability based on the environmental data as predictor variables, and presence and pseudo-absence coordinate data. Habitat suitability scores were inferred using the predict function of the package raster and a threshold value of 0.18 was estimated with the threshold function.
The estimated membership coefficient matrix generated by TESS3 (Caye et al 2015) was used to predict ancestry coefficients on the spatial grid generated by the species distribution model (SDM). This was achieved with an R implementation of Kriging interpolation in the package tess3r (Caye et al. 2015).
To further identify the role of landscape structure on the genetic variation of C. undulatum we performed the Moran's eigenvector map (MEM; Legendre et al. 2015) analysis based on genotypes and spatial coordinates implemented in the R package MEMGENE (Galpern et al. 2014). Results from MEMGENE were then superimposed on the SDM to identify spatial genetic neighbourhoods and possible barriers to pre-fragmentation gene flow (e.g. areas with environment incompatible with the presence of C. undulatum, hence falling outside our SDM; Fig. 5a).

Effective population size
The effective population size (N e ) of the genetic clusters detected was assessed based on the linkage disequilibrium method implemented in Neestimator v2 (Do et al. 2014), discarding alleles with frequencies < 0.02. The heterozygosity excess and M-ratio deficiency tests (Cornuet and Luikart 1996;Garza and Williamson 2001) as implemented in INEST v2.2 (Chybicki and Burczyk 2009) were used to infer possible occurrence of a recent bottleneck. The twophase-mutation model was used, with a proportion of multi-step mutations of 0.22. Statistical significance was assessed using the one-tailed Wilcoxon test after 1 × 10 6 permutations.

Fine-scale spatial genetic structure
Spatial autocorrelation analysis was used to assess the presence, intensity and extent of the spatial genetic structure possibly characterizing C. undulatum populations. It was performed by calculating the kinship coefficient F ij (Loiselle et al. 1995) for each pair of individuals. For each distance class, tests for the statistical significance of average F ij values were conducted by 5 × 10 3 permutations of individual geographical coordinates. Analyses were performed using 18 ad-hoc selected distance classes, from every 10 m at short spatial scale to every 1000 m at long spatial scale. The intensity of the spatial genetic structure (SGS) was measured by the Sp statistic (Vekemans and Hardy 2004), computed as Sp = b F /(F 1 -1), where b F is the regression slope of the kinship estimator F ij computed among all pairs of individuals against their geographical distances, and F 1 is the average kinship coefficient between individuals of the first distance class (0-10 m, for this calculation). An historical estimate of gene flow (gene dispersal distance parameter σ) was also obtained from the regression of pairwise kinship coefficients on the logarithm of the distance, computing b F within a restricted distance range (σ e to 20σ e , Hardy et al. 2006) and using a range of effective density (d e ) estimates from 10 ind/ha to 100 ind/ha. All analyses were performed using SPAGeDi 1.5 (Hardy and Vekemans 2002).

Association between ecology and genetics
Redundancy analysis (RDA) was performed to assess the influence of environmental characteristics on the genetic indices of populations. In this analysis, the five genetic indices (Ar 16 , H E , F IS , G ST and Jost's D) were the response variables, whereas the 1 3 following environmental variables were considered as predictors: number of plants, population area, fragment area, percentage of vegetated (both native and not native) land within a 500-m-radius area around the population centroid, current isolation index, historical connectivity index, and historical percentage of vegetated land ('historical' refers to the 1950s before large-scale clearing for urbanisation occurred within the species range; detailed methods are reported in Appendix S1 of supplementary materials). Prior to analysis, variables were scaled to zero-mean and unit-variance. Model selection was performed by assessing multicollinearity among and significance of predictors (respectively by means of the variance inflation factor, with a cut-off value of 3, and permutation tests).
Association between genetics and reproductive output was undertaken using extensive reproductive data collected in the 2017 flowering season (September-November 2017) by Delnevo et al. (2019a, b). Briefly, they obtained a total of 64,170 flowers and 3144 seeds from 198 selected plants (see Delnevo et al. 2019a, b for details on seed collection and germination methods).
To infer the relationship of genetic indices and current fragmentation conditions with plant reproductive variables, Principal Component Analysis (PCA) followed by variation partitioning (based on RDAs and partial-RDAs) were carried out. First, scaled genetic indices (i.e. Ar 16 , H E , F IS , G ST and Jost's D) and environmental variables (number of plants, population area, fragment area, percentage of vegetated land, floral display index, and current isolation index) were separately subjected to PCA in order to extract loading coefficients (i.e. scores of the populations in the multidimensional ordination space), providing a low-dimensional description of the datasets. Then, following the procedure suggested by Borcard et al. (2011), variation partitioning was performed on three variables related to plant reproduction (fruit production, seed production, and seed germination) constrained by the genetic and environmental loadings. Multivariate analyses were run in R 3.6.1 (R Development Core Team, 2018) with the "vegan" package (Oksanen et al. 2019).

Genetic diversity and structure
The levels of genetic diversity were mostly similar across populations, with the exception of populations H, M and N ( Fig. 2; Table 1). These populations showed the lowest values of allelic richness and heterozygosity and, generally, the highest differentiation from the other populations (Table 1; Fig. 3; pairwise population matrices of G ST and D in Appendix S2 of supplementary material). However, all populations carry at least one private allele, and private allelic richness varied between 0.10 in population H to 0.36 in population B. The generally even levels of genetic diversity and widespread presence of private genetic diversity reflects in the balanced percentage of individuals that should be collected from each population to maximise the allelic diversity (maxK; Table 1).
The most likely number of clusters assigned by both the ΔK and the log likelihood method was K = 4 ( Fig. A1 in Appendix S3). The results of TESS were congruent with the structure find by STRU CTU RE v. 2.3, as we also found four groups (Fig. A1 in Appendix S3). Analyses of genetic structure showed three populations (population H, and the small populations M and N) were markedly divergent from the rest of the populations that represented a genetically homogeneous group (Figs. 1, 3). Populations M and N are the only two populations where the genetic signature of a recent bottleneck was detected according to the test for deficiency in M-Ratio (Wilcoxon signed-rank test: population M, Z score = −1.72. P = 0.041; population N, Z score = −2.20. P = 0.012).

Spatial and landscape genetics
By interpolating the ancestry coefficients, four spatially defined clusters were identified (Fig. 4a). Cluster 1 stretched from the northern to the central populations (Fig. 4b). Cluster 2 was more heavily represented in the southern part of the distribution range, but also extended to north-central populations (Fig. 4c). Clusters 3 and 4 were confined to the west and east portions of the geographic range, respectively, and were less represented in the populations in the central part of the main distribution ( Fig. 4d and e). Within the main groups in the central (north-south) part of the geographic distribution we found an estimated Ne equal to 1589 [921-5082] with an isolation-by-distance (IBD) pattern where the pairwise genetic distance between populations is correlated with their geographical distance (r = 0.37; P < 0.01). This relationship becomes random when populations H, M, and N are included in the analysis (r = 0.15; P = 0.11). Consistently with the IBD result, only a small proportion of all genetic variation of C. undulatum was attributed to landscape structure (MEMGENE R 2 adj = 0.037). Despite the small proportion of genetic variation attributed to landscape patterns, the first MEMGENE variable, explaining 40.8% of the variation, found two main spatial genetic neighbourhood separated by a central gap of unsuitable environment identified by our SDM (Fig. 5a). This central gap, however, is not complete but is surrounded by suitable environment for C. undulatum. The second MEMGENE variable (explaining 29.3% of the total variation) showed an additional neighbourhood among individuals from the north to the south of the distribution separated from the westerly populations, and in particular from the western-most population H, which was separated from the main distribution by an almost complete gap of unsuitable environment (Fig. 5a). Noteworthy, the southern-most part of the historic natural bushland was already separated from the main bushland by a gap of ~ 1 km in 1953 (Fig. 1a), and the three southern-most populations (i.e. C, D, G) showed a weak sign of spatial genetic neighbourhood.
At a finer spatial scale, a significant positive spatial autocorrelation of pair-wise kinship coefficients reflecting an average, non-random spatial distribution of genotypes up to 1000 m was found (Fig. 5b). Interestingly, a clear change in the slope of the relationship between kinship coefficients and spatial distance is recognizable between the

Drivers of genetic diversity and reproductive performance
After excluding the variable number of plants from the model due to collinearity (VIF = 6.3), RDA showed that the variation in genetic indices of populations was significantly correlated to only one explanatory variable, namely the historical connectivity index (pseudo-F 1,12 = 7.46, P = 0.023). This parsimonious model, explaining 38.3% of the variation in genetic indices (based on the adjusted R 2 ), indicated that historical connectivity was negatively associated with among-population genetic differentiation and positively related to within-population genetic diversity (Fig. 6).
Model selection of RDAs performed on plant reproductive performance, constrained by the principal components (PCs) extracted by PCAs of genetic and environmental variables, retained the first genetic PC (pseudo-F 1,7 = 10.32, P = 0.009) and the first environmental PC (pseudo-F 1,6 = 12.66, P = 0.003). In particular, the first genetic PC had the highest correlation with genetic fixation index (G ST ), whereas the first environmental PC was correlated with natural vegetated land (fragment area) (Supplementary material Fig. S2). Variation partitioning of reproductive performances indicated significant influences for the whole set of sources of variation (i.e. genetic, environmental and combined fractions, explaining 65.9% of the variance in the reproductive output; F3,5 = 9.68, P = 0.001), as well as for genetic fraction plus the combined fraction (F1,7 = 10.34, P = 0.004), and the environmental fraction plus the combined fraction (F2,6 = 12.66, P < 0.001). Moreover, while the proportion of variation in reproductive performance explained by the environmental fraction alone was significant (F2,5 = 4.51, P = 0.020), the one attributed to the genetic fraction alone was only marginally significant (F1,5 = 3.21, P = 0.059).

Discussion
The relatively recent urban expansion of Perth allowed us to investigate the role of both current and historical environmental factors that have influenced the distribution of genetic diversity of C. undulatum. Genetic analysis of this threatened species revealed a relatively weak genetic structure and levels of genetic diversity and differentiation indices mainly influenced by historical levels of connectivity. We suggest that the low levels of genetic fixation and moderate differentiation in the core distribution area can be attributed to the combined effect of extensive gene flow through the pre-fragmented landscape and the overall lack of detectable reduction in population size of this long-lived resprouting plant. This indicates that recent fragmentation has yet to impact patterns of spatial genetic structure in the contemporary adult cohort. Our study also demonstrates that the reproductive performance of C. undulatum were correlated to a combination of current environmental conditions and level of genetic diversity, with important implications for its conservation.

Genetic diversity and structure
The patterns of genetic diversity and population structure found in our study were not consistent with theoretical expectations for highly fragmented populations. This result supports previous suggestions that the impact of fragmentation may depend on species-specific lifehistory traits and the time-scale of the fragmentation process (Gibbs 2001;Kramer et al. 2008). We found fixation index close to 0, moderate levels of genetic differentiation, and weak genetic structuring across the entire distribution range of C. undulatum. Although comparisons among studies involving different genetic markers have to be treated with caution, the only available study on the population genetics of C. undulatum, characterizing four populations using AFLP markers, showed similar weak genetic structuring (Close et al. 2006). The difference between fixation and differentiation indices (i.e. larger D than G ST ) might represent an early signal of failure of gene flow keeping populations genetically connected. The values of these indices depend on the accumulation of historic and recent migration and population sizes. Populations that are completely isolated from migration can still show low G ST due to the maintenance of within-population diversity if deme size is still large, whereas D could be more influenced by past variation in migration rates (Jost et al. 2018;Whitlock 2011). However, such indices should not be used to make inferences about current migration in populations which are likely not at equilibrium and, instead, should be interpreted at the empirical level to guide conservation decisions.
We also found a relatively even distribution of genetic diversity across the species' distribution. Ten out of 14 populations had similar levels of diversity and contributed a similar amount to the within-population component of total allelic diversity. Such a scenario, characterized by a rather genetically homogeneous group comprising the majority of populations, has been likely generated by extensive gene flow among populations over time. This was also supported by the small proportion of genetic variation that was attributed to landscape structure. In fact, all the currently fragmented populations of C. undulatum, with the exception of two populations (i.e. populations H and N), were part of a continuous bushland 70 years ago, only punctuated by areas of incompatible environment or cleared bushland. Peripheral populations that were historically isolated from the central distribution showed a marked genetic divergence from the main group of populations. This may be explained by the presence of a large patch of incompatible environment that spatially separated these populations from the nearest conspecific individuals. Such barriers may have kept the populations separated well before the recent fragmentation event that started in the late 1950s (Kelobonye et al. 2019), with levels of gene flow too low to avoid genetic drift. Indeed, populations H and N are separated by incompatible environment that have a minimum width of ~ 1500 m, well above the average distance that can be travelled by C. undulatum pollinators (based on the relation between body size and foraging distance in hymenopterans, Greenleaf et al. 2007).
The impact of recent fragmentation is also likely to be detected in the allelic diversity of populations that are now reduced to very small numbers of individuals (López- Cortegano   Fig. 4 Spatial interpolation of ancestry coefficients of Conospermum undulatum using Krig modelling in tesse3r package, superimposed on a shaded relief map with 10 m contours (http: //services.arcgisonline. com/ArcGIS/rest/services/Elevation/World_Hillshade/MapServer/tile/{z}/{y}/{x}). The first map a combines areas of maximal ancestry proportion for each of the four genetic clusters obtained with TESS3. b Northern, c southern, d western, e eastern clusters. In all maps, total area was trimmed using the species distribution model (SDM; darker grey area). Maps were generated with the R package raster v3.4 (https: // cran.r-project.org/web/packages/raster/index.html) 1 3 1 3 et al. 2019). In our study, population M, which was part of the main continuous bushland in 1953 but is now reduced to 10 individuals, is one of the only two populations where the genetic signature of a recent bottleneck was detected and has the lowest allelic richness among all analysed populations (Ar 16 = 3.810). Although geographically central, the genetic structure of population M differed from the main cluster of populations in the ordination space. This population also produced 20 fruits during the 2017 reproductive season, of which only two developed viable seeds.
Gene dispersal is important to counteract the possible negative effects of genetic drift, and our investigation has also provided useful information about the dispersal capabilities of C. undulatum. Pre-fragmentation dispersal was not impeded, at least within large populations, with historical estimates of gene dispersal distance of up to Correlograms from spatial autocorrelation analysis using the correlation coefficient F ij by Loiselle et al. (1995). The grey area represents the 95% CIs around the null hypothesis of absence of spatial genetic structuring, black lines around mean F ij values represent their 95% confidence intervals generated by jackknifing loci. Population genetics kinship coefficients are calculated with respects to the population average, and this explains negative values 1 3 100 m. Contemporary inter-population dispersal, on the other hand, is unlikely, given the current fragmentation and isolation of populations. The abrupt change in the slope of the kinship-distance relationship at ~ 20 m may reflect the impact of different dispersal curves between pollen and seed dispersal on fine scale spatial genetic structure. The limited, gravity-driven seed dispersal of C. undulatum (to date only inferred by field observations; Close et al. 2006) is consistent with high average kinship coefficients observed over short distances. In contrast, pollen dispersal by the native bee L. conospermi is likely responsible for genetic connectivity until the upper limit of where fine-scale structure was significant (~ 1000 m). The genetic data is consistent with the inferred foraging range of L. conospermi (~ 500 m; Houston 1989) and provides evidence for the spatial extent of pollen dispersal undertaken by this pollinator. The foraging range of L. conospermi is shorter compared with other common bird and insect pollinators encountered in similar habitats. Common pollinating birds such as honeyeaters have been tracked moving up to 12.5 km while foraging (Saunders and De Rebeira 1991). European honeybees have been shown to maintain extensive pollen flow in fragmented landscapes (Byrne et al. 2007;Dick et al. 2003); indeed, A. mellifera is able to extend its foraging range up to 9.5 km from the hive (Beekman and Ratnieks 2000). For the smaller native bee L. conospermi, built-up urban areas may generate more significant barriers to inter-population movement (Andersson et al. 2017), limiting pollen-mediated gene flow only to within continuous populations of C. undulatum. These differences in dispersal highlight the great importance of considering species-specific pollinator groups when inferring the potential impact of habitat fragmentation on contemporary gene dispersal.

Historical connectivity
In many cases, failure to detect a genetic response to habitat fragmentation may be primarily due to insufficient time having elapsed since fragmentation. This may particularly be true for long-lived plant species, since they can carry a 'loss of genetic diversity' debt for many generations after a fragmentation event (Aguilar et al. 2008;Honnay and Jacquemyn 2007;Vranckx et al. 2012). Our results for C. undulatum showed a geographical clustering of populations and a spatial distribution of genetic diversity not compatible with impeded gene flow. The genetic structure of these populations was related only to the connectivity among populations in 1953 inferred from historical records. Such effects of historical landscape structure on gene flow have also been reported for the shrub species Grevillea caleyi, another member of the Proteaceae (Llorens et al. 2004). Moreover, similar to our findings,  found that fragmented populations of the forb Globularia bisnagarica, historically belonging to the same grassland patch, still retained a spatial genetic structure mostly determined by pre-fragmentation processes. Our study indicated that the current spatial genetic structure of remnant C. undulatum populations reflects population distributions over deeper time scales. Population disturbance since 1953, has yet to impact patterns of spatial genetic structure in the contemporary adult cohort. Taken together, these results support hypotheses that plant longevity relative to time since disturbance and mean dispersal distances are key factors buffering populations from the effects of habitat fragmentation.

Reproduction
It is widely known that sexual reproduction plays a key role in migration, adaptation, and population persistence (e.g. Fenner and Thompson 2005). Recent work in C. undulatum (Delnevo et al. 2019a, b) found that reproductive output was significantly correlated with ecological consequences of habitat fragmentation, including reduction in population size, floral display and spatial connectivity among fragments. In line with these findings, in the current study we found a strong influence of fragmentation factors on the reproductive effort of C. undulatum. The addition of the genetic information of populations to the analysis highlighted that most of the variation in key stages of sexual reproduction (i.e. pollination, fertilization and germination) was explained by the combined fraction of genetic and environmental characteristics of C. undulatum populations. Although the genetic fraction alone was only marginally significant, the environmental and the combined fractions explained 56.5% of the variation, indicating a joint influence of genetic diversity and the current environmental condition.
Changes in the spatial arrangement of remnant populations following anthropogenic fragmentation can lead to altered dynamics of gene flow, with small and/or isolated populations being less visited by pollinators (Dauber et al. 2010). In addition, reduced dispersal between populations may reduce the ability of small populations to recover lost genetic variation from larger remnants. In the short term, a loss of heterozygosity following population size reduction can reduce fitness due to inbreeding depression (Charlesworth and Willis 2009). When this occurs in predominantly outcrossing plants, such as C. undulatum, the effects of inbreeding depression can impact early fitness traits such as the development of the embryo and seed germination (Aguilar et al. 2019;Vranckx et al. 2012). Our results highlight the current area of populations as an important variable influencing the reproduction of the species. A lack of contemporary inter-population pollen flow due to limited pollen dispersal may over time reduce the reproductive output in small, isolated populations. The weak but detectable effect of genetic factors on reproductive components may represent an early signal of the impact of increased inbreeding in small remnants. Given that inbreeding depression is often environmentally dependent (Yun and Agrawal 2014), it is therefore not surprising that the interaction between genetics and environmental variables best explained reproductive components.

Conservation implications
This study highlights the importance of combining population genetics, reproductive data and ecological characteristics, together with the biogeographic history of populations, to obtain a more complete picture of the impact of habitat fragmentation. Many plants may tolerate habitat fragmentation due to naturally high levels of inbreeding and low genetic diversity due to self-compatible mating systems (e.g. Llorens et al. 2018). Other species may cope with habitat fragmentation due to highly mobile pollinators able to extend their foraging range to maintain pollen dispersal between remnant populations (e.g. Byrne et al. 2007). However, for plants such as C. undulatum which are self-incompatible species and rely on small native bees for pollen flow (Delnevo et al. 2020a, b) these population changes may have a stronger impact on reproductive output. Soil stored seed banks may provide some temporal buffer to these impacts, although seed viability in the related species Conospermum triplinervium have been shown to progressively decline in viability over 15 months (Tieu et al. 2007). For C. undulatum, the combination of life-history traits, mating system and pollinators in this system may therefore provide further challenges for its ability to regenerate after disturbance in fragmented landscapes.
Therefore, avoiding further decline in the size of medium and large populations is crucial to mitigating the effects of genetic drift and maintaining genetic diversity of C. undulatum. Nonetheless, populations of C. undulatum, like those of many other highly fragmented species, are more likely to face extinction from demographic factors (González-Varo et al. 2012) long before the lack of genetic variation impedes adaptation to changing environments. Reproductive performance of C. undulatum are already affected by current fragment size and the combination of environmental and genetic features of remnant populations. Our study did not identify any clear source or sink populations to be prioritised for conservation based on their contribution to within-and between-populations allelic diversity (López-Cortegano et al. 2019). This is particularly imporant for planning conservation actions such as translocations and reintroductions in order to maintain, or increase, adequate population size to sustain the genetic diversity within and among populations. Conservation efforts aimed at reconnecting remnant populations are predicted to be especially fruitful in species with limited pollen dispersal and should be combined with ongoing efforts to monitor reproductive output.

Code availabilityd
not applicable.