Range- wide population viability analyses reveal high sensitivity to wildflower harvesting in extreme environments

1. The ecological effects of harvesting from wild populations are often uncertain, es pecially since the sensitivity of populations to harvesting can vary across species’ geographical ranges. In the Cape Floristic Region (CFR, South Africa) biodiversity hotspot, wildflower harvesting is widespread and economically important, provid ing an income to many rural communities. However, with very few species studied to date, and without considering range-


| INTRODUC TI ON
Understanding the impacts of harvesting on wild populations is central to ecology, population biology and conservation management (Allendorf et al., 2008;Beissinger & Westphal, 1998;Ticktin, 2004).
Harvesting can affect the local population dynamics and viability of exploited species (Ghimire et al., 2008;Lamont et al., 2001;Peres et al., 2003;Ticktin, 2004), causing population extirpations in extreme cases (e.g. Nantel et al., 1996). The viability of natural populations depends foremost on the key demographic rates of survival and reproduction, which vary among populations and environmental conditions across species' geographical ranges (Treurnicht et al., 2016; see also Merow et al., 2014). Range-wide population viability analyses that account for local demographic variation and environmental conditions are therefore required to improve our understanding of the effects of harvesting across large spatial extents (Campbell et al., 2018;Crone et al., 2011;Reed et al., 2002). This will promote the development of spatially adapted guidelines to safeguard exploited plant species (Jansen et al., 2018;Ticktin, 2004), which is increasingly important in the face of progressing climate change threatening the future persistence of many species (Butchart et al., 2010;Pimm et al., 2014).
Spatially replicated and multigenerational experiments are ideally required to assess the effects of harvesting on species' population dynamics across their geographical ranges, but performing such extensive studies, particularly for long-lived species, is challenging (Campbell et al., 2018;Reed et al., 2002). Hence, conservationists often rely on limited demographic data from few populations to extrapolate population viability analyses over large spatial extents (Beissinger & Westphal, 1998;Cabral et al., 2011;Morris & Doak, 2002). For plant species, population viability analyses that integrate variation in individual demographic rates, density dependence and local environmental conditions to assess the sensitivity of species and their constituent populations to harvesting are particularly lacking (see Campbell et al., 2018 and references therein for non-plant studies). This limits our understanding of inter-and intraspecific variation in sensitivity to harvesting, and how the local environment affects the sensitivity of species to wildflower harvesting across their geographical ranges.
The harvesting of wildflowers is an important and rapidly expanding industry in the Cape Floristic Region (CFR, South Africa) biodiversity hotspot, providing employment and income for many rural communities. The economic value of exported CFR wildflowers is estimated to be at least US$16 million gross income per year (Privett et al., 2020). Most target species of the wildflower industry are Proteaceae shrubs with large, attractive inflorescences that are harvested for both domestic and international markets (Turpie et al., 2003). Many Proteaceae are serotinous and have a firedependent life cycle ( Figure 1): plants accumulate canopy-stored seed banks which release seeds after fire (Bond et al., 1984;Lamont et al., 2020). Serotinous Proteaceae may thus be particularly vulnerable to wildflower harvesting which reduces the size of canopy seed banks and may limit post-fire seedling recruitment (Figure 1; Cabral et al., 2011;Lamont et al., 2001;Maze & Bond, 1996;Mustart & Cowling, 1992). To safeguard the canopy seed banks of these species, regional conservation guidelines recommend that at least 50% of flowers (reproductive units, including inflorescences or cones) should remain in the canopy after harvesting (reviewed in Van Wilgen et al., 2016). However, local-scale experimental studies have shown that harvesting 50% of flowers can compromise post-fire recruitment (Maze & Bond, 1996). Such experiments are typically restricted in spatial and temporal extent, which leads to an insufficient understanding of the impacts of wildflower harvesting across species' geographical ranges. This hampers the development of guidelines to manage the CFR's wildflower resource sustainably, such as setting harvesting limits that vary by species and/or geographical location (Pressey et al., 2007;Turpie et al., 2003;Van Wilgen et al., 2016).
Here, we examined spatial variation in sensitivity to wildflower harvesting across the geographical ranges of 26 Proteaceae species in the CFR. We built on data and analyses of range-wide demographic variation across environmental gradients (Treurnicht et al., 2016) to parameterise models of environment-dependent population dynamics for the study species . In a population viability analysis, we then simulated the effects of 0% and 50% wildflower harvesting on population extinction risk as a function of local environmental variables across species' geographical ranges.
To evaluate the sensitivity of population viability to harvesting, we compared population-level probabilities of extinction with and Our findings caution against the application of general harvesting guidelines irrespective of species, geographical location or local environmental conditions. Our range-wide population viability analyses provide insights for developing speciesspecific, spatially nuanced guidelines for conservation management. Our approach also identifies species and areas to prioritise for monitoring to prevent the overexploitation of harvested species. without wildflower harvesting for the 26 study species. Specifically, we investigated: (a) how species differ in sensitivity to harvesting, (b) how sensitivity to harvesting varies within the geographical range of each species and (c) whether harvesting sensitivity of local populations is associated with environmental conditions. This serves to identify the species and populations most impacted by harvesting.
The resulting understanding of how range-wide demographic variation and local environmental conditions affect harvesting sensitivity for multiple plant species should inform spatially explicit conservation management to prevent the overexploitation of species.

| Study region and study species
The CFR is characterised by a Mediterranean-type climate, high topographic and environmental variation with predominantly nutrient-poor soils (Allsopp et al., 2014). The CFR constitutes fire-prone shrublands where many plant species depend on periodic fires that occur on average every 10-21 years (Kraaij & Van Wilgen, 2014) and burn most above-ground biomass (Bond & Van Wilgen, 1996). Proteaceae are a characteristic family of mostly overstorey woody shrubs that contribute substantially to the region's plant diversity (Manning & Goldblatt, 2012). CFR Proteaceae are functionally diverse and regarded as model organisms for ecological research (Schurr et al., 2012).
We studied 26 Proteaceae species from the genera Leucadendron (dioecious) and Protea (hermaphrodite) that are endemic to the CFR and differ in their geographical ranges (see Table S1 in Supporting Information). All study species are serotinous and have a fire-driven life cycle (Figure 1), accumulating long-lived canopy-stored seed banks (not soil-stored seed banks) and retaining seeds in fire-protected woody cones until fire occurs (Bond et al., 1984;Lamont et al., 2020;Treurnicht et al., 2016Treurnicht et al., , 2020. Fire generally kills adult plants and triggers seed release, seed dispersal and subsequent seedling recruitment. Populations thus regenerate and establish depending on the availability of seeds (Bond & Van Wilgen, 1996;Treurnicht et al., 2016).
Seedling recruitment is confined to the first few years after fire and once recruited plants are more than 3 years old, they generally experience low mortality between fires (Lamont et al., 2020;Treurnicht et al., 2016). While fire is the predominant cause of mortality for adult plants (Bond & van Wilgen, 1996), the 26 study species fall into two distinct fire response types: adults of resprouting species (n = 7) frequently survive fire by regrowing from fire-protected buds, whereas adults of nonsprouting species (n = 19) are generally fire-killed (Clarke et al., 2013;Rebelo, 2001;Treurnicht et al., 2016).
Experiments with serotinous Proteaceae (including three of our study species) showed that removing reproductive units (inflorescences or cones) by wildflower harvesting reduces the size of the canopy seedbank (Mustart & Cowling, 1992;Witkowski et al., 1994; see also Figure 1), and does not trigger compensatory responses of plant individuals (such as increased inflorescence production or increased investment in the remaining inflorescences or seeds).
Wildflower harvesting may even reduce subsequent inflorescence and cone production, while reducing seed set and increasing seed predation in the remaining cones (Mustart & Cowling, 1992).
However, evidence for such reinforcement of negative harvesting impacts is rather weak and varies between species. We therefore assume that harvesting a fraction of a plant's inflorescences causes a proportional reduction of seed numbers in the canopy seed bank.

F I G U R E 1
The fire-driven life cycle of serotinous Proteaceae depends on key demographic rates of adult fecundity (size of the canopy seed bank), post-fire seedling recruitment and adult fire survival (blue-grey boxes). Wildflower harvesting impacts this life cycle by reducing the number of seeds in the canopy seedbank. Environmental effects on demographic rates and changes in population density (green text) affect the sensitivity of populations to wildflower harvesting

| Demographic data collection
The study species' fire-driven life cycle ( Figure 1) facilitates measurements of key demographic rates of total adult fecundity, post-fire seedling recruitment and adult fire survival that can be directly related to stand age, thereby describing how these demographic rates change with post-fire age (Treurnicht et al., 2016. Adult fecundity since the last fire (the size of the individual canopy seed bank) was determined as the product of total closed cones per plant and number of fertile seeds per closed cone. In recently burnt populations, post-fire recruitment per pre-fire adult was quantified as the ratio between post-fire recruits and all pre-fire adults (fire survivors or burnt skeletons). Additionally, adult fire survival was measured in recently burnt populations as the ratio between post-fire surviving adults and all pre-fire adults. To estimate density dependence of demographic rates, intraspecific population density was recorded as the density of either alive adult plants (fecundity), all pre-fire adults (survival) or alive post-fire adults (recruitment). Demographic sampling was designed to cover variation in population densities and range-wide variation in environmental conditions (climate factors, soil fertility and fire return interval), which are major drivers of the study species' population dynamics (for details see Treurnicht et al., 2016). In total, we assembled 3,617 population-level records of demographic rates (68-370 records per species with a median of 85; see Pagel et al., 2020).

| Stochastic simulation models of local population dynamics
We developed a stochastic simulation model that describes the local and population density (D t = N t /A). Furthermore, the simulation model accounts for stochasticity in each demographic process.
For each time step, the number of adults that survive the fire is drawn from a binomial distribution with a probability parameter π.surv(X, T, D t ): The generation of new recruits depends on total seed production and seed establishment. The total number of seeds in the population's canopy seed bank is drawn from a Poisson distribution, where the mean λ t is a product of the per-individual fecundity μ.fec (X, T, D t ) and the number of pre-fire parents (adjusted by the sex ratio p.fem in the case of dioecious Leucadendron species): Wildflower harvesting is implemented as a proportional reduction of the canopy seed bank by a harvesting rate f H : Remaining seeds can then establish as new recruits with a probability described by the establishment rate π.recr(X, SD t , AD t ): In contrast to adult survival and fecundity, seed establishment rates are not affected by pre-fire adult densities, but by the density of available seeds (SD t = Seeds.act t /A) and surviving adult densities

| Species-specific model parameterisation
Parameterisation of the stochastic simulation model was based on the 3,617 population-level records of demographic rates across species ranges Treurnicht et al., 2016).
With these data, we estimated species-specific functional responses of the key demographic rates (fire survival rate: π.surv; fecundity: μ.fec; per-seed establishment rate: π.recr) to population density, fire return interval (T) and four environmental covariates (X 1,…,4 ; see Table S2 for interspecific variation in estimated maximum demographic rates and Pagel et al. (2020) for details on parameter estimation and the functional form of demographic responses). Fire history information (time since last fire and length of the previous fire interval) for the demographic study sites was deduced from either field observations minimum temperatures in the month of July (T min , °C), (c) the mean of daily maximum temperatures in the month of January (T max , °C) and (iv) soil nutrient status (or 'soil fertility'), which was represented by an index (ranging from 0 to 10) that incorporates soil texture and base status (Schulze, 2007).

| Simulations to quantify variation in sensitivity to harvesting
We investigated inter-and intraspecific variation in sensitivity to wildflower harvesting by analysing population viability responses to different harvesting regimes across the geographical ranges of the 26 study species. For each species, we identified all grid cells (1′ × 1′) across the study region in which the species is known to occur (Protea Atlas Project, Rebelo, 2001) and parameterised the stochastic simulation models with grid-specific values of environmental covariates. Simulation models used sequences of fire return intervals (T) drawn from a Weibull probability distribution per grid cell, as predicted by a data-driven model of geographical variation in fire return intervals for the study region (Wilson et al., 2015). For each population, we considered scenarios of no harvesting (f H = 0) and 50% harvesting (f H = 0.5), which corresponds to a regional conservation guideline that 50% of reproductive units (inflorescences Population viability responses to harvesting were quantified as the probability (P 100 ) of a population that is initially in the quasistationary (or established) phase (Grimm & Wissel, 2004) to become extinct within a 100-year period. To obtain P 100 from our stochastic simulations, we used the ln(1 − P 0 )-method which is based on the theory of Markov processes (Grimm & Wissel, 2004).
For each combination of study species, grid cell and harvesting scenario, we performed 10,000 replicated simulations. Each simulation started with an initial population size of N 0 = 1,000 individuals in an area of A = 10,000 m 2 and ran until an extinction occurred (or up to a maximum of 100,000 years). From the simulated extinction times, we then calculated P 0 (t), the cumulative probabilities of extinction until time step t (with t ranging from 1 up to 100,000 years). The inverse slope of a regression of − ln(1 − P 0 (t)) against t is the population's intrinsic mean time to extinction, T m (Grimm & Wissel, 2004). T m thus measures the mean time to extinction of a population in the quasi-stationary phase.
As the fundamental currency of population viability, T m is independent of the initial population size used in simulations (N 0 ; Grimm & Wissel, 2004). We then used T m to calculate the extinction probability over 100 years (P 100 ; see above) as: P 100 = 1 − exp(−100/T m ) (Grimm & Wissel, 2004). P 100 thus represents extinction probability over a timeframe relevant for conservation management and red list categorisation (IUCN, 2019; criterion E. Quantitative Analysis). Figure 2 provides an example of a simulation model and shows how T m and P 100 were evaluated to derive the responses of population viability to harvesting.

| Inter-and intraspecific variation in sensitivity to harvesting
To assess interspecific variation in sensitivity to harvesting, we summarised the proportion of populations of each species that correspond to different 'extinction risk categories' under the two harvesting scenarios. We categorised population-level extinction risk as VERY LOW (P 100 < 0.1), LOW (P 100 ≥ 0.1), INTERMEDIATE (P 100 ≥ 0.2) or HIGH (P 100 ≥ 0.5), which broadly corresponds to IUCN threat categories 'NT', 'VU', 'EN' and 'CR', respectively (criterion E. Quantitative Analysis; IUCN, 2019). For each species, sensitivity to wildflower harvesting was then summarised as the proportion of populations that shift to a higher extinction risk category when subjected to 50% wildflower harvesting.
To assess intraspecific and spatial variation in sensitivity to harvesting, we calculated the change in extinction probability (∆P) as the difference between the extinction probabilities of the two harvesting scenarios: ∆P = P 100 (50%) − P 100 (0%). This follows the definition of Morris and Doak (2002) to evaluate sensitivity across multiple populations from population viability and extinction risk analyses. To identify geographical locations across the range of each study species where populations are more sensitive to harvesting, we then highlighted the 1′ × 1′ grid cells where ∆P ≥ 0.1.
To investigate the relationship between sensitivity to harvesting and environmental variation, we fitted linear regressions that described the population-level responses of ∆P (log-transformed) to three climate variables (AI, T min and T max ), soil fertility and mean fire return interval. All environmental variables were scaled and centred to ensure comparability across variables and species. We included both linear and quadratic terms in the maximal models to allow for the possibility of monotonic or unimodal response curves to environmental gradients. We used automated model selection (r package MuMIn; Barton, 2015) among all combinations of explanatory variables. The best model for each species was determined by the lowest sample size corrected Akaike information criterion (AICc; Burnham & Anderson, 2002). We evaluated these best models to show the response of ∆P to each environmental gradient across the study region and summarised these relationships for the 26 study species. The shape of an environmental relationship (positive, negative, unimodal (negative quadratic) or U-shaped (positive quadratic) effects) was determined from whether the coefficient matrix of each best model contained a quadratic effect of an environmental variable and from the sign of the coefficients.

| Intra-and interspecific variation in sensitivity to harvesting
We detected varying degrees of intraspecific variation in sensitivity to harvesting among our 26 study species (Figure 3; Figure S1).
While for some species the majority of populations were largely indifferent to 50% harvesting (i.e. remained in the 'VERY LOW' and 'LOW' risk categories; Figure 3a), other species had a considerable proportion of populations that appeared sensitive to wildflower harvesting (shifted to 'INTERMEDIATE' and 'HIGH' risk categories; Figure 3b). In addition, some species had a small proportion of their populations at high risk of extinction even in the absence of harvesting (e.g. top-right quadrats of Figure 3c).
There was substantial interspecific variation in sensitivity to harvesting, including nonsprouters (n = 19) and resprouters (n = 7; Figure 4). The mean proportion of populations that shifted to a F I G U R E 2 Example of replicated stochastic simulations showing the dynamics of one population of Protea repens in response to 0% (blue lines) and 50% (orange-red lines) wildflower harvesting respectively. The time of a local extinction event (X; years) was recorded per simulation replicate. The population-level intrinsic mean time to extinction (T m , following Grimm & Wissel, 2004) and the probability of extinction over 100 years (P 100, representing a timeframe relevant for conservation management) were calculated for each harvesting scenario. These two measures were evaluated to represent population viability and to evaluate the population-level sensitivity to wildflower harvesting. Note that only a subset of 10 replicated simulations per harvesting scenario are depicted in this figure. The full simulation study used a time horizon up to 100,000 years to improve the estimation of the intrinsic mean time to extinction (T m ) for populations with low per-generation extinction probability and comprised 10,000 replicate simulations per harvesting scenario for each population of the 26 Proteaceae study species (see Section 2 for details) F I G U R E 3 Intraspecific variation in sensitivity to wildflower harvesting for three exemplary Proteaceae study species. P 100 0% (x-axis) and P 100 50% (y-axis) are estimated population-level extinction probabilities over 100 years in response to 0% and 50% harvesting, respectively. In each plot, dots represent populations in different environments across the geographical ranges of each study species. The diagonal line (grey) represents no change in extinction probabilities across populations due to harvesting, whereas a greater deviation of points from this line indicates a higher sensitivity (∆P) to harvesting for the respective populations. Dashed blue lines indicate thresholds between the different extinction risk categories (VERY LOW (P 100 < 0.1), LOW (P 100 ≥ 0.1), INTERMEDIATE (P 100 ≥ 0.2) or HIGH (P 100 ≥ 0.5)) and numbers give the percentages of populations that either remain in a category or move to a higher category due to harvesting. See Figure S1 for all study species higher extinction risk category in response to 50% harvesting was 12% across all species. Four nonsprouter species, however, had a high proportion of populations for which extinction risk increased substantially under 50% harvesting (>35% in Figure 4). Nine nonsprouter and two resprouter species had >10% of their populations shifted to a higher extinction risk category due to harvesting (47% F I G U R E 4 Interspecific variation in sensitivity to wildflower harvesting for 26 Proteaceae species. Individual bars (blue: nonsprouters (n = 19); grey: resprouters (n = 7)) show the proportion of populations per species that shift to a higher risk category due to 50% wildflower harvesting (see Figure 3). where the change in population-level extinction probability ∆P > 0.1. ∆P was calculated as the difference between extinction probabilities under 0% and 50% harvesting (see Figure 2). The white and black areas depict species-specific occurrence records and the geographical distribution of CFR Proteaceae, respectively (Rebelo, 2001). See Figure S3 for all 26 study species and 29% of the total number of nonsprouters and resprouters, respectively). Similar results were obtained when comparing species in terms of the proportion of populations for which harvesting increases extinction probability over 100 years by more than 10% ( Figure S2).

| Geographical and environmental variation in sensitivity to harvesting
Intraspecific variation in sensitivity to harvesting showed distinct geographical patterns across the ranges of the 26 study species ( Figure 5; Figure S3). Sensitivity to harvesting was often high at the edges of species' geographical ranges. These areas of high sensitivity tended to cluster in the arid, north-western part of the study region, especially for species with large geographical ranges (Figure 5a,d; see also Table S1 for species-specific range size). In contrast, other species were sensitive to harvesting in either the central or eastern parts of their geographical range, while a few species with restricted ranges showed a consistently uniform response across their geographical range (Figure 5b; Figure S3).
Range-wide intraspecific variation in sensitivity to harvesting was driven by environmental variation across the study region. The predominantly U-shaped responses along most environmental gradients show that sensitivity increased towards the environmental extremes of species' geographical ranges ( Figure 5). Populations that are associated with more extreme environmental conditions are thus more sensitive to harvesting than populations in more moderate environments. These responses were particularly common for two climate variables, AI and T max (Figure 6; see also Table S1).

| D ISCUSS I ON
We developed range-wide population viability analyses to quantify spatial variation in sensitivity to wildflower harvesting across the global geographical ranges of 26 Proteaceae species in the CFR biodiversity hotspot. Our novel approach demonstrates the importance of combining spatial demographic data, density-dependent population dynamics and local environmental variables when assessing the sensitivity of species to harvesting. We detected considerable inter-and intraspecific variation in sensitivity to harvesting for our study species. Intraspecific variation in harvesting sensitivity associated with environmental variation across species' ranges. Here, we discuss the mechanisms that potentially shape population-level variation in sensitivity to harvesting, provide insights for conservation management and identify priorities for the future monitoring of harvested species.
The substantial inter-and intraspecific variation in sensitivity to harvesting that we detected provides evidence for range-wide variation in seed-limitation for the 26 study species. High sensitivity to harvesting should notably arise when seedling recruitment is limited by an insufficient seed crop (Maze & Bond, 1996;Nathan & Muller-Landau, 2000). In contrast, low sensitivity to harvesting implies that population viability and persistence are less limited by seed availability.
Species that were relatively indifferent to 50% wildflower harvesting, notably a few nonsprouters with relatively large geographical ranges, show consistent seed saturation across their geographical ranges. Seed saturation occurs if successful germination and seedling recruitment are more limited by post-fire habitat conditions than by the populationlevel availability of seeds (Lamont et al., 1993).
The high intraspecific variation in sensitivity to harvesting across species' geographical ranges ultimately results from spatial variation in key demographic rates (Treurnicht et al., 2016;see also Merow et al., 2014). Interestingly, range-wide intraspecific variation in sensitivity may arise from 'demographic compensation'. Demographic compensation occurs when different demographic rates (such as fecundity and per-seed recruitment) respond to environmental gradients in opposing ways. This reduces variation in population growth rates (Doak & Morris, 2010;Villellas et al., 2015). However, demographic compensation does not necessarily reduce variation in seed limitation and sensitivity to harvesting. In fact, demographic compensation will increase variation in harvesting sensitivity when opposing environmental responses of fecundity (or seed production) F I G U R E 6 Shape of relationships between the sensitivity to wildflower harvesting (∆P, calculated as the difference between extinction probabilities under 0% and 50% harvesting) and five environmental variables for the 26 Proteaceae study species (see also Table S1). Barplots show the number of species for which relationships were estimated as positive (+), negative (−), unimodal (∩), U-shaped (∪) or not significant (ns) in a linear regression of population-level sensitivity against local environmental variables: an aridity index (AI = P/(T + 10°C), mean daily minimum temperatures in July (T min , °C), mean daily maximum temperatures in January (T max , °C), soil fertility (Schulze, 2007) and mean fire return interval and per-seed recruitment probability change the limiting roles of these demographic processes across different environments. High sensitivity to harvesting therefore arises when fecundity is low and recruitment is essentially density independent (seed limitation), whereas populations with high fecundity and strong negative density dependence of recruitment show low sensitivity to harvesting (seed saturation).
The distinct geographical patterns of high sensitivity to harvesting at the range edges for species with large geographical ranges likely occur because range limits often coincide with climatic extremes, constraining individual growth and reproduction (Holt, 2003;Sexton et al., 2009). In contrast, species with smaller, more restricted geographical ranges showed high sensitivity throughout their range.
These species may experience climatic and environmental limits throughout their relatively small geographical ranges. In turn, a few species had nearly uniform sensitivity across their range in the absence of harvesting: high-altitude, montane species had >10% of their populations at high risk of extinction (e.g. P. punctata; Figure 5c; Figure S1). Species occurring in these specialised environments are likely constrained across most of their geographical range by slow growth rates associated with short growing seasons, harsh climatic conditions and specific requirements for seedling recruitment (Mustart et al., 2012;Treurnicht et al., 2016).
Our range-wide population viability analyses for multiple species provide insights for developing species-and area-specific harvesting guidelines. In particular, a few range-restricted species had a large proportion of populations at high risk of extinction even without harvesting and were severely vulnerable to 50% harvesting (e.g. L. modestum; P. compacta; Figures S1-S3). This corresponds to these species' current endangerment status (Raimondo et al., 2009)  Future investigations of harvesting impacts should ideally focus on long-term monitoring and multigenerational experiments across species' geographical ranges (Ticktin, 2004). We identified areas and associated growing conditions where populations are likely more vulnerable to harvesting, providing guidance for establishing long-term experiments and the monitoring of populations across the study region. Furthermore, our results highlight that endangered and range-restricted species have a substantial proportion of populations at high risk of local extinction even without harvesting (e.g. P. compacta; Figure S3). These species are good candidates for monitoring the impacts of ongoing global change in addition to harvesting, especially since a species like P. compacta is an economically important resource for the local wildflower industry (Treurnicht, 2010;Van Wilgen et al., 2016). Wildflower harvesting is also spatially and temporally heterogeneous, and the consistent removal of 50% flowers in every generation is likely to be rare.
Whether spatial and temporal heterogeneity in harvesting may alleviate impacts on local population viability, however, requires further investigation .
The detected high sensitivity to harvesting of certain species concurs with previous studies (Cabral et al., 2011;Maze & Bond, 1996;Mustart & Cowling, 1992). Further experiments are needed to explore species-specific responses to harvesting and whether some species show reinforcement of negative harvesting impacts (Mustart & Cowling, 1992;Witkowski et al., 1994). The potentially resulting information on non-proportional responses of fecundity to harvesting could then easily be included in our modelling framework. Future studies should also consider that colonisation of distant habitat patches is inevitably more seed-limited than local recruitment (Cabral et al., 2011;Nathan & Muller-Landau, 2000). Consequently, wildflower harvesting may have negative consequences for the long-term persistence of species at the landscape-scale that are not evident in local populations (Cabral et al., 2011). Integration of our local population models with spatially explicit models of large-scale dynamics (Cabral et al., 2011) would help to further understand such effects of harvesting across large spatial extents.
Ongoing climate change is a major driver of biodiversity loss and species range shifts (Pecl et al., 2017;Urban, 2015), yet local population extinctions and the exploitation of species are often overlooked in forecasts of global change impacts. Exploring the relationships between range-wide population viability, range limits and ecological traits, and conducting trait-based assessments of species' vulnerability to ongoing environmental change is particularly important (Estrada et al., 2016;Treurnicht et al., 2020).
Notably, the patterns of sensitivity to harvesting that we detected in the north-western CFR ( Figure 5) coincide with regional assessments of climate change vulnerability owing to changing precipitation and fire regimes in the study region (Midgley et al., 2003;Wilson et al., 2015;Yates et al., 2010). The consistent high sensitivity to harvesting at the hot and arid extremes of species' geographical ranges ( Figure 6) stresses the urgency to jointly consider the potential impacts of harvesting and global change drivers in future investigations.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data and code available from the Dryad Digital Repository https:// doi.org/10.5061/dryad.547d7 wm7v .
Additionally, details of the demographic data can be found in the cited articles Treurnicht et al. (2016) and Pagel et al. (2020) and supporting information of these articles. These data were partly used under license agreements from provincial and national conservation organisations in South Africa (CapeNature and SANParks).