Reproductive biology and population structure of the endangered shrub Grevillea bedggoodiana (Proteaceae)

Narrowly endemic species are particularly vulnerable to catastrophic events. Compared to widespread species, they may also be less capable of adapting to shifts in environmental pressures as a result of specialisation on a narrow range of local condition and limited ability to disperse. However, life-history traits, such as preferential outcrossing and high fecundity can maintain genetic diversity and evolutionary potential, and boost species resilience. The endangered Grevillea bedggoodiana (Enfield Grevillea) is an understorey shrub restricted to an area of ca. 150 km2 in south-eastern Australia with a legacy of large-scale anthropogenic disturbance. Prior to this study little was known about its biology and population structure. Here, its breeding system was assessed through a controlled pollination experiment at one of its central populations, and eight populations were sampled for genetic analysis with microsatellite markers. The species was found to be preferentially outcrossing, with no evidence of pollination limitation. In most populations, allelic richness, observed heterozygosity and gene diversity were high (Ar: 3.8–6.3; Ho: 0.45–0.65, He: 0.60 − 0.75). However, the inbreeding coefficients were significant in at least four populations, ranging from Fi -0.061 to 0.259 despite high outcrossing rates. Estimated reproductive rates varied among sampled populations but were independent of gene diversity and inbreeding. Despite its small geographic range, the species’ populations showed moderate differentiation (AMOVA: FST = 0.123), which was largely attributable to isolation by distance. We interpret these results as suggesting that G. bedggoodiana is reproductively healthy and has maintained high levels of genetic diversity despite recent disturbance.


Introduction
Effective conservation of threatened species requires a thorough understanding of multiple interacting aspects of their biology (Aitken and Whitlock 2013;Allendorf et al. 2013;Scheele et al. 2018). The 'endangered' status is afforded to species that are at high risk of extinction, with one of the criteria for listing being acute restriction of the species' geographic range. The rationale for this is that most threatening processes are spatially correlated, meaning that if all individuals of a species occur in a one small area, a single catastrophic event can impact its entire population (IUCN 2012). Moreover, compared to widespread congeners, geographically restricted species tend to have narrower niche breadth (Gaston et al. 1997;Slatyer et al. 2013) and be less genetically diverse (Gitzendanner and Soltis 2000;Cole 2003). When faced with novel challenges, such species may be unable to draw on standing genetic diversity to evolve and adapt (Aitken et al. 2008;Anderson et al. 2012). This is a precarious situation because narrow endemics tend to be poor dispersers (either due to lack of traits facilitating dispersal and/or because of being confined to a habitat 'island') (Gaston 1996), and so are unlikely to be able to track the conditions of their niche. However, the evaluation of extinction risk relating to the species' extent of occurrence is strongly context dependent. Different species with similar 1 3 ranges may be vulnerable to different threats depending on their life history traits such as lifespan (Thrall et al. 2014), breeding (Willi et al. 2005) and mating systems Nora et al. 2016;Ritchie et al. 2019), reproductive rates (Rymer et al. 2005) and seedbank dynamics (Kenny 2000;Pickup et al. 2003), as well as genetic diversity (Frankham 2015) and population size (Leimu et al. 2006;Caballero et al. 2017). Moreover, to add to the complexity of the task, all of the above should be considered in the context of the species' landscape, disturbance history and interspecific interactions, particularly, pollination ecology (Cranmer et al. 2012;Llorens et al. 2013).
In order to persist and evolve, outcrossing plant species need to maintain genetic diversity and produce enough viable offspring to replace the adults. At the level of local populations the plants' ability to do so depends on mate availability, but, in the case of animal-pollinated species, also on the local abundance and effectiveness of pollinators (Richardson et al. 2000;Devaux et al. 2014; Barrett and Harder 2017). Scarcity of effective pollinators may cause reduced levels of ovule fertilisation (pollination limitation) and reduced rates of outcrossing. This can lead to reduced seed set and/or increased inbreeding, potentially leading to low reproductive rates and low offspring fitness (inbreeding depression) (Copland and Whelan 1989;Wilcock and Neiland 2002;Coates et al. 2007;Bezemer et al. 2019). Moreover, for species without adaptations for long-distance seed dispersal, the maximum foraging range of the pollinators is expected to limit gene flow between geographically distinct populations (England et al. 2001). Low levels of gene flow can make small populations vulnerable to loss of genetic diversity through drift (Frankham 2015) and reduce the likelihood of the spread of adaptive traits (Lowe and Allendorf 2010;Aitken and Whitlock 2013;Ottewell et al. 2016). The knowledge of reproductive health and levels of inter-population gene flow is, therefore, crucial for evaluating resilience of a threatened species and deciding on appropriate conservation action. However, this knowledge cannot be gained solely by measuring seed set or by inference from population genetic parameters such as inbreeding coefficients and heterozygosity. This is particularly true for long-lived plants which may be reproducing sporadically but may employ reproductive strategies such as staggered germination of soil-stored seedbanks which may buffer local populations against loss of genetic diversity. Consequently, assessments of the reproductive health of plant populations should be based on estimates of actual reproductive rates followed by tests of whether they are correlated with measures of heterozygosity and inbreeding (Szulkin et al. 2010;Breed et al. 2015). This is relevant for two reasons. First, because the number of seeds produced by the adult plants in a population is a poor indicator of their reproductive potential as not all seeds will evade granivores, germinate and establish.
For example, in Proteaceae, most of the seeds produced are commonly lost and do not contribute to the next generation (Lamont and Groom 1998;Auld and Denham 1999). Second, while there is evidence that low fitness is associated with low population genetic diversity and high inbreeding, the reverse is not always the case (Reed and Frankham 2003;Frankham 2015).
Grevillea is a large and diverse genus of woody shrubs and trees in the Proteaceae (Makinson 2000). A wide array of reproductive strategies are known in grevilleas from full self-compatibility (Richardson et al. 2000;England et al. 2001;Gross et al. 2012) to obligate outcrossing (Hoebee and Young 2001;Caddy and Gross 2006;Holmes et al. 2008Holmes et al. , 2009, and at least two species, G. infecunda (James et al. 2017) and G. renwickiana (James and McDougall 2014), are thought to be sterile and reproduce exclusively by rootsuckering. The pollination strategies are also varied in this genus with native insect, bird (usually of the Meliphagidae), and mammal pollinators involved (Olde and Marriott 1994). The introduced honeybee, Apis mellifera, has become a common floral visitor to many grevilleas, leading to concerns over their potential to disrupt pollination and mating systems of native species (Vaughton 1996;Paton 2000;Smith and Gross 2002;Celebrezze and Paton 2004;Whelan et al. 2009). Pollination disruption is of particular concern for the primarily bird-pollinated forest understory grevilleas of south-eastern Australia, owing to habitat fragmentation and, in recent decades, the decline of forest birds (Ford et al. 2001;Ford 2011).
Since the early 2000s, several studies have assessed the genetic diversity and population structure of grevilleas (Hoebee and Young 2001; England et al. 2002;Llorens 2004;Nistelberger et al. 2015), including some species of the Grevillea aquifolium group (the 'holly-leaved grevilleas') from south-eastern Australia (Holmes et al. 2008(Holmes et al. , 2009James and McDougall 2014;James et al. 2017). Several species in this group appear to be reproductively constrained as a result of pollinator and/or mate limitation (Holmes et al. 2008), variation in ploidy (Holmes et al. 2009) and sterility (Kimpton et al. 2002;James et al. 2017). This study focused on a narrowly endemic holly-leaved grevillea, G. bedggoodiana J.H.Willis ex McGill. This species is restricted to a small area of Palaeozoic sedimentary soils in central Victoria surrounded by soils derived from Quaternary volcanic activity (Joyce 2004). There are no historical collections of this species outside of its current range (Australasian Virtual Herbarium, accessed 08 June 2022) suggesting that the newer igneous soils are unsuitable for its growth. Currently, all Victorian grevilleas are protected under the state's Flora and Fauna Guarantee Act (1988), and the largest populations of G. bedggoodiana are managed for conservation within Enfield State Park. However, in the recent past (ca. 160 years), much of its habitat has been impacted by mining and logging while the surrounding landscape of fertile basaltic plains have been largely cleared for agriculture. This relatively recent disturbance likely caused population bottlenecks and disrupted pollination, potentially leaving a legacy of reduced genetic diversity which could affect the species' reproduction and long-term evolutionary potential. Grevillea bedggoodiana is listed as 'endangered' (Flora and Fauna Guarantee Act, 1988) but prior to this study it had not been a focus for research. Consequently, we sought to establish some of the fundamental aspects of its reproductive biology and population genetics, specifically, (1) to determine the species' breeding system through a controlled pollination experiment; (2) to assess the genetic diversity and structure of its populations through microsatellite-based analyses, and (3) to test for bottlenecking and reduction of reproductive rates in populations with lower-than-average genetic diversity.

Site and study species
Grevillea bedggoodiana is endemic to an area of approximately 150 km 2 comprised of Enfield State Forest (ESF) and Enfield State Park (ESP) (here, collectively referred to as Enfield Forest) in Victoria, Australia (Fig. 1). Since the 1850s areas of Enfield Forest were disturbed by gold mining (until the 1930s) and logging (until the 1960s) (Parks Victoria 1998). During the last 70 years, fires have occurred in Enfield Forest every 15-20 years. In 1995, a large wildfire burnt through about 90% of the Forest (Parks Victoria 1998). Subsequently, fuel reduction burns were conducted in some parts of the forest, particularly near major roads and infrastructure (Irena Cassettari, Department of Environment, Land, Water and Planning, pers. comm.). The consequences of these multiple disturbances for the local flora and fauna assemblages are unknown. However, they likely contributed to the decline of the native vertebrate pollinators, (Ford et al.  (Young et al. 1996).
Grevillea bedggoodiana is an understory shrub of eucalypt-dominated open forest, growing on gravelly clay soils, predominantly on ridges and north-facing slopes in the southern and western parts of Enfield Forest. It tends to form discrete stands of several hundred individuals (Carter et al. 2006), with scattered individuals occurring occasionally particularly on roadsides (SW, pers. obs.). In the context of this study, the term 'population' is applied following a standard definition of a spatially discrete group of organisms where interbreeding is expected to occur primarily within the groups but not excluding the possibility of gene flow between groups. This is consistent with the National Recovery Plan for Grevillea bedggoodiana (Carter et al. 2006) which defined important populations of this species as those of over 1000 individuals (Carter et al. 2006). There are 11 such populations, containing an estimated total of 20,000 individuals (ca. 50% of all known plants), and all of them are within ESP (which has been established specifically for their protection) (Carter et al. 2006) (Fig. 1). The populations in ESF are considerably smaller and tend to be separated by greater distances (SW, pers. obs.). Consequently, this study focused primarily on the large populations in ESP (Table 1).
Grevillea bedggoodiana is a low to prostrate shrub, typically to 0.5 m tall with trailing lateral branches up to 1.9 m long (Fig. 2a), sometimes more upright. It flowers in spring from late October through November. Its conflorescences are secund (toothbrush-like) (Makinson 2000) and consist of 10-20 pairs of flowers ( Fig. 2b) which open sequentially over 5-7 days (SW, pers. obs.). The individual zygomorphic flowers are protandrous. During the male phase of a flower, pollen is transferred to, and presented at the tip of a 12-16.5 mm-long style. Removal of pollen by visiting animals exposes the stigma, which becomes viscous and receptive to pollen ca. 2-3 days after the beginning of anthesis (SW, pers. obs.).
Animal visitors to G. bedggoodiana conflorescences were surveyed through direct observation and camera trapping in a concurrent study (Aristidou and Hoebee, unpublished data). The vertebrate visitors included one bird: the yellow-faced honeyeater (Lichenostomus chrysops) and several mammal species: antechinus (Antechinus sp.), eastern pygmy possum (Cercartetus nanus) and common brushtail possums (Trichosurus vulpecula). The invertebrate visitors of this species include native bees, wasps, and hoverflies (Emily Noble, Field Naturalists Club of Ballarat, pers. comm.). However, the most common day-time floral visitor at the sites during the study was the introduced European honeybee (Apis  Fig. 2c) contain up to two seeds, which are approximately 6 mm long when mature (Fig. 2d). The species is non-serotinous, and its seeds are thought to be dispersed primarily by gravity (Makinson 2000). However, the presence of an elaiosome ( Fig. 2d) suggests that myrmecochory may contribute to secondary seed dispersal over short distances (Auld 1995;Edwards and Whelan 1995;Auld and Denham 1999). It is not known whether the seeds can germinate without fire/ smoke (Carter et al. 2006) and how long they can persist in the soil. Field observations suggest that wild plants grown from seed are able to flower and set seed within five years (SW, pers. obs.). Older plants form a lignotuber and may be able to resprout after fire, but do not appear to reproduce clonally by root suckering (SW, pers. obs.).

Breeding system
A controlled pollination experiment was conducted at MCS ( Fig. 1) between 25th of October and 4th of November 2018. Ten healthy adult plants with unopened flowers and growing at least 10 m apart were selected. Four pollination treatments were applied to each plant: autogamy (spontaneous selfing), geitonogamy (animal vector-mediated self-pollination), xenogamy (cross-pollination), and open pollination (control). The conflorescences in the autogamy treatment (two per plant) were bagged in bud using fine mesh organza bags to prevent visitation by pollinators and otherwise not manipulated. Those in the geitonogamy (one per plant) and xenogamy (two per plant) treatments were bagged in bud, had self-pollen removed shortly after anthesis and were hand-pollinated using either pollen from other flowers on the same plant (geitonogamy), or pollen from two different plants growing more than 10 m away (xenogamy). The open pollination conflorescences (three per plant) were tagged but otherwise not manipulated, then bagged several days after the last flowers had opened. The unbalanced design of this experiment was dictated by the amount of pollen that could be collected for geitonogamous pollination without affecting other treatments, while the use of three conflorescences in the open pollination treatment aimed to provide a more accurate estimate of natural levels of seed set. The developing confructescences remained bagged until collection on 16th of December 2018. At that time, the follicles were welldeveloped but not yet open. They were counted, stored at room temperature for four days and then opened manually to assess seed set. The numbers of filled seeds were compared between pollination treatments using paired t-tests in base R (R Core Team 2021).

Sampling, DNA extraction and genotyping
Eight representative populations were selected for the population genetic study: six in ESP and two in ESF (Table 1, Fig. 1). Fresh leaf tissue was collected from 30 adult plants at each of the eight sites (240 in total) -selecting plants growing at least 5 m apart and avoiding nearest neighbours to reduce the chance of sampling potential clones (Holmes et al. 2009) or closely related or individuals (Krauss et al. 2009). Sampled leaves were placed in individual paper envelopes and desiccated using silica gel beads (ChemSupply, Adelaide, Australia). The GPS coordinates of each plant sampled were recorded and representative voucher specimens for each population were collected and deposited at LTB with duplicates sent to MEL. For each sample, 20 mg of dried leaf tissue was disrupted using a Tissue Lyser II (Qiagen, Hilden, Germany). DNA was then extracted using a DNeasy Plant Mini Kit (Qiagen, Hilden, Germany) following the manufacturer's protocol. The concentration of the DNA isolates was measured using a NanoDrop Lite (Thermo Fisher Scientific, Waltham, MA, USA). All samples were stored at − 20 °C until required.
The plants were genotyped using the three primer method implementing four fluorescently labelled universal primers (Schuelke 2000;Blacket et al. 2012) and fourteen sets of primers designed to amplify microsatellite loci in two closely related species (Grevillea aquifolium and G. infecunda; James et al. 2012). DNA was amplified in two multiplex reactions, each consisting of seven primer pairs. Sequence tails of the universal primers were added to the forward primers, such that each multiplex PCR amplified up to two loci labelled with the same fluorescent primer avoiding size overlap of the expected amplicons (James et al. 2012) (Table S1). Multiple Primer Analyzer (Thermo Fisher Scientific, Waltham, MA, USA) was used to predict potential primer dimerization. A third multiplex reaction was designed to most efficiently re-amplify samples which initially failed to amplify or produced ambiguous amplicon peaks for some loci (Table S1). A Type-it Microsatellite Kit (Qiagen, Germany) was used to amplify the loci in 15 μL reactions containing 10 − 20 ng of template DNA, 0.1 μM of each of the tailed forward primers, 0.2 μM of each of the reverse primers and 0.2 μM of the FAM, VIC, NED and PET fluorescent labelled primers. Polymerase chain reactions (PCR) were performed using BioRad MyCycler™ Thermal Cycler (BioRad, Hercules) with a denaturation step at 95 °C for 5 min, followed by 30 cycles of 95 °C for 30 s, 57 °C for 90 s, and 72 °C for 30 s, followed by a final extension step at 60 °C for 30 min. Resultant amplicons were processed by the Australian Genome Research Facility (AGRF, 1 3 Melbourne) using LIZ500 as an internal standard. Geneious Prime v.2019.0 (Biomatters Ltd, Auckland, New Zealand) was used in-house for manual genotype determination.

Population genetics analyses
Linkage disequilibrium between loci and the frequencies of null alleles were estimated using Genepop (Rousset 2008; R Core Team 2013). Basic population genetic parameters, including number of alleles per locus (N a ), number of private alleles (A priv ), observed heterozygosity (H o ), gene diversity (expected heterozygosity, H e ) were calculated using GenAlEx 6.5 Smouse 2006, 2012). Allelic richness (A r , number of alleles per locus with rarefaction to the smallest sample size) were calculated using SPAGeDi (Hardy and Vekemans 2002). SPAGeDi was also used to estimate the individual inbreeding coefficients (F i , with a test of significance based on 999 permutations) and population outcrossing rates (1-s) based on standardized identity disequilibrium using the kinship coefficient (Loiselle et al. 1995). The 95% confidence intervals for A r , H o , H e and 1-s were approximated by multiplying the standard errors of the estimates by ± 1.96. Bottleneck 1.2.02 (Piry et al. 1999) was used to test for excess of heterozygotes (evidence of recent population bottlenecks) employing the Wilcoxon signed rank test (owing to the low number of loci) with sequential Bonferroni correction (Holm 1979) and the infinite allele model (IAM) of marker evolution. Other models (stepwise mutation and two-phase) could not be used because our data included irregular amplicons (presumably resulting from mutations in the microsatellite flanking region England et al. 2002;Holmes et al. 2009;Putman and Carbone 2014)).
Population differentiation was assessed through F STbased analysis of molecular variance (AMOVA) (Excoffier et al. 1992) as implemented in GenAlEx. Resulting values of F ST from this approach are equivalent to θ of Weir and Cockerham (1984) (Holsinger and Weir 2009 (Rousset 1997;Séré et al. 2017). The genetic structure of the sampled populations was also assessed using Structure 2.3.4 (Pritchard et al. 2000(Pritchard et al. , 2009 testing K values of 1-8 (1,000,000 MCMC steps, with 100,000 dememorization cycles, replicated ten times for each K value). Analyses were run under the admixture model of ancestry with sampling sites as priors and the correlated allele frequencies model. Structure Harvester (http:// taylo r0. biolo gy. ucla. edu/ struc tureH arves ter/) was used to select the model that best explained the data using the Evanno (Delta K) method (Earl and VonHoldt 2012). Discriminant analysis of principal components (DAPC, adegenet) was used as an alternative method to assess the clustering of the populations according to individual plants' genotypic similarity (Jombart, 2008;R Core Team, 2013).

Regression analysis: reproductive rates vs gene diversity and inbreeding
Reproductive rates were estimated for each of the eight populations by the number of juveniles (i.e., plants with flexible stems, fewer than about 40 leaves and not flowering) per 100 adults (N j ) in two 2 × 50 m belt transects. The transects were selected such that they passed through areas of typical plant density at each site (avoiding roadsides and population edges). The counts of adult plants in the transects were also used to estimate the total number of adult plants at the sampled populations (N est ). This was done by extrapolating the number of plants in the transects to the total area occupied by each population, where the area was estimated using QGIS v3.4.8 (QGIS Development Team 2017) as half that of the minimum convex polygons enclosing the plants sampled for the genetic analyses. Linear regression (R Core Team, 2013) was used to test for correlation between N j and gene diversity (H e ) and inbreeding coefficients (F i ) of the adult populations. Because reproductive rates are potentially related to time since fire, detailed fire history data was obtained for Enfield Forest area from DataShare (https:// datas hare. maps. vic. gov. au/ search; 'Fire history overlay of most recent fires'). The study sites were overlaid with the fire history data using QGIS v3.16.1, and linear regression used to test for correlation between population N j values and time since last fire. The R package 'gvlma' (Peña and Slate 2006) was used to validate the regression models.

Microsatellite data
All of the 14 primer pairs used amplified G. bedggoodiana DNA. However, five loci (AQ11, AQ14, AQ34, IN32, IN41) were excluded from analyses due to missing or ambiguous PCR products for more than 75% of individuals in any population. The remaining nine loci were successfully amplified in at least a third of the samples (10 individuals) in all populations. Evidence of linkage was detected in one out of 36 loci pairs (IN18-IN34) when testing across all populations (Chi-sq > 49.6, df = 16, p < 0.001, which remained significant after Bonferroni correction for 36 simultaneous tests, lowest adjusted α = 0.00139, all other p-values > 0.05). However, per population tests showed significant linkage between these two loci in only two of the eight populations (LGR and VR) and they were retained in the data set.
The frequencies of null alleles across all loci and populations were low: mean = 0.10 ± 0.02 SE (median 0.06). However, locus-by-locus analysis showed high frequencies of null alleles in two out of nine loci, IN12: mean = 0.21 ± 0.06 SE (median 0.18), and IN17: mean = 0.25 ± 0.07 SE (median 0.17). Excluding these loci from the data had negligible effect on estimates of allelic richness and expected heterozygosity (H e ) but biased upward the estimates of observed heterozygosity (H o ) in seven out of eight populations (mean difference: + 0.06) and considerably lowered the mean population fixation indices in all populations, except OT (F i , mean difference: − 0.0875). It also lowered the estimates of global F ST by 0.03 and increased by 10% the coefficient of determination (R 2 ) in the Mantel test of isolation by distance, but it had negligible effect on the topography of the inferred population structure (i.e., clustering of the populations). In the interest of retaining the information provided by these loci, both were retained in the data set, and two alternative estimates of F were reported (with and without IN12 and IN17).

Genetic diversity
Across the nine loci, a total of 107 alleles were identified (7 − 25 alleles per locus). The allelic diversity was distributed unevenly among the populations, with 1 − 6 private alleles per population and mean allelic richness (A r ) ranging from 3.8 to 6.3, mean 5.25 (after rarefaction to n = 10) ( Table 2). The population estimates of observed and expected heterozygosity (H o and H e ) also varied considerably (from 0.45 to 0.63 and from 0.60 to 0.75 respectively). The inbreeding coefficient (F i ) ranged from -0.061 to 0.259 (0.149-0.260, when estimated excluding loci IN12 and IN17) and were significantly greater than 0 for all sampled populations except VR, when estimated with all 9 loci, and for four out of eight populations when estimated excluding IN12 and IN17 (Table 3). The estimates of outcrossing rates (1-s) were high for all populations, ranging from 0.90 to 1.0 (Table 2). Bottleneck analysis with sequential Bonferroni correction indicated significant excess of heterozygotes only at IGR (two-tailed Wilcoxon's test, p = 0.002 compared to α = 0.006).

Population structure
AMOVA indicated that the populations were moderately differentiated with 12% of the genotypic variation due to among-population differences (F ST = 0.123, p = 0.001). Pairwise F ST values ranged from 0.012 to 0.286 (all p-values < 0.01) and were consistently highest for all pairs involving the small and peripheral VR population (Table 3). The Mantel test indicated moderate and statistically significant correlation between genetic differentiation and geographic distances between sampled populations (R 2 = 0.49, p = 0.01) (Fig. 4); a similar correlation was found for standardised F ST *.
Structure Harvester's delta K plot indicated that the models with three and five clusters explained the data nearly equally well (Fig. 5a). In the K = 3 model, the small peripheral VR population appeared genetically distinct while the remaining populations formed two clusters intermixing at the MCS and MCN in the centre of ESP. The K = 5 model further assigned distinct ancestry to the OT, LGR and IGR with the remaining WT, HR, MCS and MCN forming a broad fifth cluster variously intergrading with their geographic neighbours but virtually no gene flow between the populations that were ≥ 6 km apart (e.g., OT and VR) (Fig. 5b). The associated tree plots show relative genetic distances between the clusters (Fig. 5c).
The results of DAPC agreed with those of Structure and indicated genetic divergence of the two smaller populations from ESF (IGR and VR) from a main cluster consisting of the six larger populations from ESP (Fig. 6a). Removing VR and IGR from the analysis, allowed to resolve finer structuring of the remaining populations, whereby the two smaller and somewhat peripheral OT and LGR populations appeared to diverge from the four central populations (Fig. 6b).

Reproductive rates vs gene diversity and inbreeding
The estimates of population areas ranged from 2500 and 22,000 m 2 and population sizes (N est ) from 300 to 4800 adult plants (Table 4). Reproductive rates (N j ) ranged from 5 to 93 juveniles per 100 adults (Table 4). Because N j for IGR was a clear outlier, the data was log-transformed to achieve linearity. Regression analysis revealed no significant correlation between log 10 (N j ) and gene diversity (H e ) (adjusted R 2 = − 0.50, df = 6, 95% CI: − 0.89 -0.31) or between log 10 (N j ) and population mean inbreeding coefficients (F i ) estimated based on both full data set (9 loci) (adjusted R 2 = -0.04, df = 6, 95% CI: − 0.89 -0.31) and without IN12 and IN17 (adjuster R 2 = − 0.06, df = 6, 95% CI: − 0.73 -0.67). The fire history data for the area revealed that two the study sites have been control-burnt since the 1995 wildfire (OT and WT, 10-9 and 7-6 years ago, respectively). A third site (WT) was partially burnt 14 years ago (Table 4). The differences in time since last fire did not appear to influence N j values (adjusted R 2 = − 0.001, df = 6, p = 0.36). The models satisfied the assumptions of linear regression.

Discussion
This study aimed to determine the breeding system of the endangered Grevillea bedggoodiana, and to assess the genetic diversity and population structure across most of its range. It also aimed to quantify the reproductive output of the eight sampled populations, and to test whether it correlated with the estimates of gene diversity and inbreeding coefficients. Based on controlled pollinations of ten plants in a large central population, the species was found to be incapable of reproduction via autogamy. The small number of seeds produced through geitonogamy indicated the species is strongly and preferentially outcrossing but capable of self-fertilisation (at least at the MCS population). This result was supported by the genetic analyses which suggested high outcrossing rates, despite low-to-moderate levels of inbreeding in at least four of the populations studied. Open pollination and hand cross pollination resulted in equivalent levels of seed set indicating that seed production was not limited by pollen availability. Thus, the low number of seeds produced relative to total numbers of flowers per conflorescence (typical in Proteaceae; Hermanutz et al. 1998) is indicative of the conflorescence carrying capacity and/or or resource limitation. Allelic richness, heterozygosity and gene diversity of the populations were mostly high. There was a moderate degree of differentiation among the sampled populations, which could be largely attributed to isolation by distance. Analyses exploring population structure indicated that the large central populations in Enfield State Park were well connected but there was considerably less gene flow among the populations at the periphery of the species range and private alleles were detected in all populations. Numbers of seedlings per hundred adults were not Fig. 6 Discriminant analysis of the principal components (DAPC) scatter plots: a including all populations sampled; b excluding the two outlying populations from ESF for better resolution of the structure of the ESP populations. In both analyses retained were the first 30 principal components (dimensions of variation in the genotypic data) and all discriminant functions (synthetic variables optimizing the clustering of groups) correlated with gene diversity or inbreeding coefficients. This result suggests there was no reduction in reproductive rates in populations with lower-than-average gene diversity or elevated estimates of inbreeding as could be expected if they were declining due to inbreeding depression (Szulkin et al. 2010). Realised breeding systems are known to vary among populations in a number of angiosperms (Whitehead et al. 2018) including Grevillea (Hermanutz et al. 1998;Richardson et al. 2000). Likewise, the flowering intensity, abundance of pollinators and reproductive output are likely to vary both spatially and temporally (e.g., G. repens, a species closely related to G. bedggoodiana, Holmes et al. 2008; see also Copland and Whelan 1989). Consequently, the results of the controlled pollination experiment presented here need to be interpreted with caution as they are based on a small sample of plants from a single, central location and the experiment was conducted during a single flowering season. Moreover, the viability and fitness of the seeds produced through geitonogamy was not assessed. However, preferential outcrossing with a possibility of low rates of selfing (mixed mating) was supported by the high estimates of multilocus outcrossing rates in all populations. The proportion of initiated follicles that matured was high in G. bedggoodiana compared to, for example, G. barklyana (Vaughton 1996), in which consistently fewer than half of all initiated fruit developed to maturity. Furthermore, the presence of juveniles at all study sites indicates the species produces viable seeds capable of germinating at least sporadically under favourable conditions and without requiring fire/smoke cues (given that the most recent major wildfire in Enfield Forest occurred 23 years prior to this study and most of the populations have not been burnt since then). This is in contrast to for example G. barklyana (Vaughton 1998) which has very low recruitment rates in the absence of fire, and several other species in which germination is strongly cued by fire and smoke (Kenny 2000).
The floral characters of G. bedggoodiana suggest the species should be primarily bird-pollinated (Olde and Marriott 1994;Cronk and Ojeda 2008;Willmer 2011) whereby the flowers have no obvious scent, do not reflect UV light (Aristidou and Hoebee, unpublished data) and produce copious amounts of somewhat dilute nectar (°Bx: 25.8 ± 0.9 SE, n = 6 conflorescences; SW, unpublished data). Although, classical pollination syndromes may not always reliably predict actual primary pollinators (Waser et al. 1996;Rosas-Guerrero et al. 2014). A single bird species (the Yellow-faced honeyeater; Lichenostomus chrysops) was recorded visiting flowers during two weeks of remote camera trapping, there were also visits by small nocturnal mammals which may also be effective pollinators (Aristidou and Hoebee, unpublished data). However, the frequency of visitation by vertebrates was very low (only 16 records in total) and the most common floral visitor was the European honeybee (Apis mellifera, Aristidou and Hoebee, unpublished data), which was large enough to occasionally touch the pollen presenters/stigma while foraging for nectar (SW, pers. obs.). Thus, the species' potential vertebrate pollinators appear to have been effectively replaced by the honeybees. This is consistent with the well-documented decline of the woodland-and forestdependent birds in south-eastern Australia (Ford et al. 2001;MacNally et al. 2009;Ford 2011) and with the finding that selective exclusion of vertebrates from G. bedggoodiana conflorescences had no effect on seed set (Aristidou and Hoebee, unpublished data). However, recent studies have raised concerns regarding the consequences of pollination by honeybees for reproduction and gene flow of native plants adapted to vertebrate pollination. Due to their high mobility and limited grooming of pollen, avian pollinators appear to play an important role in maintenance of genetic diversity and long-distance pollen dispersal of species they pollinate (Hoebee 2002;Breed et al. 2015;Krauss et al. 2009b;Bezemer et al. 2019;Kestel et al. 2021). For example, evidence from the self-compatible Grevillea macleayana (another species likely adapted to bird-mediated pollination) suggests that honeybees are more likely than birds to distribute pollen among flowers within individual plants, potentially reducing outcrossing rates (England et al. 2001). While in self-compatible population of G. barklyana, selective exclusion of vertebrate pollinators (allowing honeybee visitation) resulted in 50% reduction in seed set (Vaughton 1996). Our results of high outcrossing rates and low geitonogamous seed set in G. bedggoodiana, suggest that its mating system at the population level should be comparatively resilient to disruption by honeybees. Table 4 Estimated population area, estimated population size (N est ), estimated number of juveniles per 100 adults (Nj) and time since last fire at the sampled populations of Grevillea bedggoodiana a These populations were control-burnt in two stages (using a road as a fire break) or only a portion of the population was burnt since the 1995 wildfire. In each case, the lower of the two values was used in the regression analysis 1 3 The genetic diversity of G. bedggoodiana was not uniformly distributed over its geographic range. The three largest populations (WT, MCS and MCN) had higher allelic richness than the two smallest populations (IGR and VR; Table 2). The observed heterozygosity and gene diversity also tended to be higher in the larger populations, however, the differences were not significant. OT appeared to break the trends in the data. Despite its relatively large size, its genetic diversity parameters, particularly allelic richness, were closer to those of the two small population from ESF (IGR and VR; Table 2), suggesting it may have gone through a recent bottleneck (Allendorf et al. 2013). This is consistent with OT and IGR bearing evidence of severe disturbance associated with historic gold mining (abandoned mine shafts and mounds of disturbed soil). Our Bottleneck analysis detected significant excess of heterozygotes at IGR, but not at OT. However, we note that because the short generation time and largely overlapping generations of G. bedggoodiana (the plants reach sexual maturity within five years), the bottleneck may have occurred too long ago (up to ca. 160 years prior to this study) to be reliably detected by this method (Allendorf et al. 2013). Nonetheless, the presence of G. bedggoodiana at these highly impacted sites, suggests that the species is either able to recover after severe and prolonged (several decades) disturbance, or can colonise disturbed areas relatively quickly. Additionally, the lack of correlation between the reproductive rates and gene diversity suggests that even the least genetically diverse populations were not affected by inbreeding depression. In fact, contrary to expectation, the population with the highest ratio of juveniles to adult plants (indicative of high fitness of the adult plants) was the visibly disturbed, small, and putatively bottlenecked IGR population. While it must be noted that the power of our test was sufficient only to detect large effect size with high probability (80%, https:// www. stats kingd om. com/ sample_ size_ regre ssion. html), this result suggests that G. bedggoodiana is capable of vigorous reproduction, possibly stimulated by soil disturbance breaking seed dormancy (Edwards and Whelan 1995;Pickup et al. 2003;Holmes et al. 2008). Taken together, our results suggest that the sampled populations of G. bedggoodiana are reproductively healthy, resilient to disturbance and have life history traits likely to promote genetic diversity, such as preferential outcrossing and overlapping generations (due to the staggered germination of soil-stored seeds). However, the species may still be at risk of reproductive failure should large wild-fires occur repeatedly in short succession (England et al. 2003). Frequent small fires such as the prescribed fuel reduction burns undertaken by the land managers in parts of the species range may also impact negatively on the species.
The genetic diversity and population structure of other species of Grevillea has been studied by several authors (Hoebee and Young 2001;England et al. 2002England et al. , 2003Hoebee et al. 2008;Holmes et al. 2009;Nistelberger et al. 2015;James et al. 2017;Kanjere and Hoebee, unpublished data). The methodological differences among these studies (molecular markers used, population genetic parameters reported and variation in sample size) make comparisons with our results difficult. Nonetheless, some comparisons with closely related congeners can be made. Grevillea bedggoodiana appeared to be less genetically diverse that either of the two grevilleas for which the microsatellite markers used in this study were originally developed: G. infecunda, a clonal narrow endemic, and G. aquifolium, a species with considerably wider distribution (H e = 0.78 and 0.79 respectively, compared to 0.68 for G. bedggoodiana) (James 2012). However, G. bedggoodiana was considerably more diverse than two other narrowly endemic holly-leaved grevilleas: G. obtecta (mean H e = 0.57 (Kanjere and Hoebee, unpublished data)) and G. repens (mean H e = 0.58 (Holmes et al. 2009)). Despite this, the population inbreeding coefficients (F i ) were on average higher in G. bedggoodiana than in the latter two species (mean F i = 0.030 and 0.024, respectively, compared to 0.164 for G. bedggoodiana). This is likely a consequence of the differences in the three species' breeding systems: G. obtecta and G. repens are self-incompatible (Kanjere and Hoebee, unpublished data; Holmes et al. 2008), while G. bedggoodiana, although preferentially outcrossing, appears to be capable of low levels of selfing.
Given the restricted distribution of G. bedggoodiana and considering that birds and potentially other vertebrates appear to contribute to its pollination, one could reasonably expect the species to be panmictic, or nearly so. However, our results indicate that the species' populations are moderately structured with at least three distinct clusters of ancestry detected by the Structure and DAPC methods. It must be noted that Structure may overestimate the number of genetic clusters when inbreeding or isolation by distance are significant (Pritchard et al. 2000(Pritchard et al. , 2009. In such cases the authors of the program note that the may not be a single 'correct' value of K and the Structure results need to be interpreted in the context of the species biology (Pritchard et al. 2000(Pritchard et al. , 2009. Nonetheless, DAPC analysis of the main cluster of central five populations suggested further structuring even among the populations separated by only a few kilometres. Compared to two other holly-leaved grevilleas for which microsatellite-based population genetic studies were undertaken, the global F ST for G. bedggoodiana was intermediate: considerably lower than that reported for G. repens -a species with strongly disjunct distribution (Holmes et al. 2009) -but higher than for G. obtecta, which has a continuous distribution across a similar geographic range and showed little genetic differentiation (Kanjere and Hoebee, unpublished data). Our results also suggest that considerable amount of genetic diversity of G. bedggoodiana is attributable to variation among its small peripheral populations, as indicated by AMOVA and presence of private alleles in all populations (Table 2). This suggests a degree of reproductive isolation. However, we acknowledge that our data does not allow us to infer inter-population migration rates (Holsinger and Weir 2009;Lowe and Allendorf 2010). Given the importance of nectarivorous birds for long-distance pollen dispersal in other species (Hoebee 2002;Llorens et al. 2012;Phillips et al. 2014) and their current decline (Ford et al. 2001), future studies should test whether avian and/or mammal pollinators are essential and currently sufficiently abundant to maintain gene flow and demographic connectivity across the range of G. bedggoodiana (e.g., Sork and Smouse 2006;Krauss et al. 2009;Lowe and Allendorf 2010;Bezemer et al. 2016).

Implications for conservation
Currently, conservation management of Grevillea bedggoodiana is focused on protecting its largest populations in Enfield State Park (ESP). The management plan for ESP (Parks Victoria 1998) outlines the key conservation actions which include implementation of a fire management program to mitigate the risk of large-scale wildfires and control of invasive weeds and pest animal species (Parks Victoria 1998). It also postulates maintenance of genetic diversity of the native plant communities. Our results indicate that the most genetically diverse populations of G. bedggoodiana are already protected within ESP. However, the smaller populations in ESF were found to be considerably distinct. Conservation management should include provision of greater protection to the ESF populations to ensure maintenance of their unique genetic diversity. Moreover, while our results suggest G. bedggoodiana can flower and reproduce within six years post-fire, populations may require more time to replenish their soil-stored seedbanks and to build up their genetic diversity. In relation to this, our result of significant isolation by distance leads us to recommend that managers allow at least five years between burning adjacent patches of the species habitat to allow unburnt older plants an opportunity to breed with the younger plants grown from seeds in the burnt areas. This would increase the degree of generational overlap and mitigate the risk of loss of genetic diversity through drift. This is particularly relevant for the small, geographically isolated populations which are less likely to receive pollen migrants from the large core populations. Finally, the apparent decline of the vertebrate pollinators, is a cause for concern. While the European honeybees appear to be effective pollinators of this species, pollination solely by honeybees may be insufficient for maintenance of high outcrossing rates and gene flow among populations. Consequently, an integrated conservation action should aim to protect and increase local abundance of native vertebrate pollinators, particularly nectarivorous birds.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.