The interplay between hunting rate, hunting selectivity, and reproductive strategies shapes population dynamics of a large carnivore

Abstract Harvest, through its intensity and regulation, often results in selection on female reproductive traits. Changes in female traits can have demographic consequences, as they are fundamental in shaping population dynamics. It is thus imperative to understand and quantify the demographic consequences of changes in female reproductive traits to better understand and anticipate population trajectories under different harvest intensities and regulations. Here, using a dynamic, frequency‐dependent, population model of the intensively hunted brown bear (Ursus arctos) population in Sweden, we quantify and compare population responses to changes in four reproductive traits susceptible to harvest‐induced selection: litter size, weaning age, age at first reproduction, and annual probability to reproduce. We did so for different hunting quotas and under four possible hunting regulations: (i) no individuals are protected, (ii) mothers but not dependent offspring are protected, (iii) mothers and dependent offspring of the year (cubs) are protected, and (iv) entire family groups are protected (i.e., mothers and dependent offspring of any age). We found that population growth rate declines sharply with increasing hunting quotas. Increases in litter size and the probability to reproduce have the greatest potential to affect population growth rate. Population growth rate increases the most when mothers are protected. Adding protection on offspring (of any age), however, reduces the availability of bears for hunting, which feeds back to increase hunting pressure on the nonprotected categories of individuals, leading to reduced population growth. Finally, we found that changes in reproductive traits can dampen population declines at very high hunting quotas, but only when protecting mothers. Our results illustrate that changes in female reproductive traits may have context‐dependent consequences for demography. Thus, to predict population consequences of harvest‐induced selection in wild populations, it is critical to integrate both hunting intensity and regulation, especially if hunting selectivity targets female reproductive strategies.


| INTRODUC TI ON
Harvesting of wild animal populations serves economic, cultural, and management purposes, but when exerted at a high rate, it can threaten population persistence (Jackson et al., 2001) and induce trait changes in life history, morphology, and behavior (Palkovacs et al., 2018). Human harvest constitutes a unique form of "predation" that fundamentally differs from "natural predation," because harvest mortality is often higher than natural mortality and not always directed toward individuals that are most vulnerable to natural mortality (Allendorf & Hard, 2009;Darimont et al., 2015;Festa-Bianchet, 2003). Because of this, human harvest has emerged as an important driver of trait change in the wild (Darimont et al., 2009;Palumbi, 2001), inducing selective pressures that vary both in strength and in direction, depending on harvest levels and practices, as well as on the phenotypes being targeted (Darimont et al., 2015). Harvestinduced selection on life-history, morphological, and behavioral traits has been documented in both fishery and hunting systems (Allendorf & Hard, 2009;Boyce, 1981;Leclerc et al., 2017;Palkovacs et al., 2018;Van de Walle et al., 2018).
In addition to its direct effect on population growth rate, harvest can affect population structure and induce changes in phenotypic traits and behavior, which thus indirectly influence population growth rate (Gosselin et al., 2015;Milner et al., 2007;Pelletier et al., 2007). Harvest-induced selection on traits linked to female reproductive performance is likely to have the greatest impact on the dynamics and persistence of populations. For instance, overfishing and harvest-induced selection on body mass and size, strong drivers of individual performance in fish, are expected to lead to earlier sexual maturation at smaller sizes and reduction in population biomass (Jørgensen et al., 2007). Regardless of the mechanism generating them, phenotypic changes within fish populations have also been shown to have larger scale impacts (Fenberg & Roy, 2008). For example, the reduction in body size and egg production in Pacific salmon (Oncorhynchus spp.) over the past 60 years in Alaska resulted in reduced marinederived nutrient transport inland, with consequences for ecosystem functioning (Oke et al., 2020).
In long-lived mammals, hunting-induced selection commonly affects male secondary sexual traits, such as antlers, horns (Coltman et al., 2003;Jachmann et al., 1995;Pigeon et al., 2016), and body mass (Tenhumberg et al., 2004). However, hunting-induced selection on these traits is likely to have limited consequences for population dynamics of such species due to the weak correlation between body mass (or correlated secondary sexual traits value), and reproductive performance in mammals (Kuparinen & Festa-Bianchet, 2017), compared with fishes. Hunting-induced selection acting directly on female reproductive traits has a much greater potential to affect population dynamics (Rughetti & Festa-Bianchet, 2014;Servanty et al., 2011), although it has rarely been investigated (but see Proaktor et al., 2007;Rughetti & Festa-Bianchet, 2014). Therefore, a step forward in our understanding of the large-scale impacts of harvest would be facilitated by an evaluation of the demographic effects of changes in female reproductive traits in general, but especially in long-lived species.
Female reproductive traits are often the target of harvestinduced selection due to harvest intensity, regulations, and harvester preferences. Selection on female reproductive traits is generated by nonexclusive mechanisms that may act simultaneously. On the one hand, high rates of mortality should select for faster life histories and favor individuals that invest earlier and more into reproduction (Olsen et al., 2004;Stearns, 1992). This may explain why wild animal populations of the same species experiencing different levels of mortality often show contrasting life-history strategies (Servanty et al., 2011;Zedrosser et al., 2011). Modeling and empirical studies have revealed that an increase in extrinsic mortality can select for earlier age at maturation, higher probability to reproduce, and increased litter/clutch size (Olsen et al., 2004;Proaktor et al., 2007). On the other hand, the nonrandom and systematic removal of specific phenotypes from a population due to hunting can also generate selection toward "shielding" traits (i.e., traits that afford a certain level of protection to individuals).
Hunting regulations often aim at directing the harvest toward (or away from) individuals with specific traits within a population to achieve a management goal, for example, to manipulate the population growth rate. As such, hunting regulations can create harvest biases and, intentionally or not, induce selectivity (Bunnefeld et al., 2009;Festa-Bianchet, 2003;Hengeveld & Festa-Bianchet, 2011;Leclerc et al., 2016;Mysterud, 2011). In the case of reproductive traits, such hunting selectivity can affect the fitness pay-off of different female reproductive tactics (Rughetti & Festa-Bianchet, 2014; Van de Walle et al., 2018).
A common practice in the management of large mammal populations is to protect the female segment of the population to ensure population viability, because the survival and reproduction of prime-aged females have the greatest potential to affect population growth, size, and fluctuations therein (Gaillard et al., 1998;Pelletier et al., 2011). In species where it is difficult to differentiate between females and males from a distance, protection of females is often achieved through the protection of family groups (Miller, 1990), as males generally do not provide parental care in mammals (Clutton-Brock, 1991). In addition, the killing of mothers and dependent young may cause ethical concerns that often motivate the extension of the legal protection to dependent young. Family groups shield ANR, Grant/Award Number: ANR-16-EBI3-0003; BearConnect, Grant/Award Number: 96/2016; Naturvårdsverket K E Y W O R D S brown bear, harvest regulation, harvest-induced selection, multistate population model, population dynamics, reproductive strategies their members against hunting under such regulations, and selection on traits increasing the duration and frequency of the formation of family groups (e.g., age at first reproduction, reproductive rate, weaning age) can be expected. Even in the absence of regulations to protect offspring, hunters sometimes voluntarily refrain from killing members of a family group (Nilsen & Solberg, 2006;Rughetti & Festa-Bianchet, 2011). Moreover, when hunting of dependent offspring is allowed, producing offspring has also been suggested to shield mothers against hunters, as hunters will shot offspring first (Ericsson, 2001). Despite the widespread application of protective measures for the reproductive segment of hunted populations, there is still little empirical and theoretical evidence of their consequences for population dynamics.
Changes in environmental conditions or perturbations, such as hunting, can influence populations via feedback mechanisms (Lachish et al., 2020). Feedback loops in populations occur when F I G U R E 1 (a) Diagram showing the multi-level processes implied in the population dynamics of the brown bear population in southcentral Sweden. At the hunting level, hunting quotas directly affect hunting pressure on bears. Hunting regulations, by dictating which categories of individuals can and cannot be hunted, affect the ratio between the number of bears unavailable and the number of bears available for hunting (feedback level). This ratio further affects hunting pressure on bears available for hunting. Hunting pressure on available bears affects their survival rate, which, combined with reproductive rates (individual level), affects both population growth rate and population structure (population level). The population is further regulated by a feedback loop (feedback level) as changes in population structure can affect the ratio of unavailable to available bears and feedback on hunting pressure on available bears. (b) Annual life cycle of the brown bear. Transitions between age and state stages occur over one year, from den emergence at time t until den emergence at time t +1. Transitions from a given stage to another are represented by solid arrows, and indirect contributions of a given stage to the cub stage through fertilities are represented by dotted arrows. Definitions: 0 = female cubs, 1 = female yearlings, 2i = independent two-yearold females, 2d = dependent two-year-old females, 3 = three-year-old females, Al =solitary adult females, A0 = adult females with cubs, A1 = adult females with dependent yearlings, A2 = adult females with dependent two-year-olds. (c) Periodic (seasonal) life cycle graph of the brown bear with three seasons: (1) spring (April-July; mating period), (2) fall (August-October; hunting period), and (3) winter (November-March; denning period). Individuals transit between the seasonal states conditional on survival, but also on reproductive rates through the probability to be weaned early or to wean early yearlings in the spring, the probability to give birth and emerge from the den the following spring with cubs as three-year-olds, and the probability to emerge from the den in the spring as adults with cubs. The indirect contribution of females to the cub (0) stage the following year is represented in gray. Yearlings (1) are all dependent (d) upon their mother at den emergence in the spring but can become independent (i) if they are weaned in the spring or if their mother dies in the following fall or winter. Dependent two-year-olds (2d) at den emergence are weaned in the spring and are thus independent (2i) in the fall along with other two-year-olds that were already weaned as yearlings in the previous year. Adult females (A) can become solitary (Al) the following season if they have weaned their cubs (in the spring only) or have lost their litter (all seasons) of cubs or yearlings. Hunting pressure is calculated as hunting mortality rate (h) based on the proportion of individuals unavailable for hunting at the beginning of the fall (hunting season), which is conditional on the state of individuals (dependent vs independent for subadults, and solitary vs with offspring for adults) at this time demographic rates depend on current population properties (e.g., population size or composition). As the latter change in time, due to for example environmental changes, so does demographic rates and population dynamics (Kokko & López-Sepulcre, 2007). Examples of feedbacks in demography include density dependence (Coulson et al., 2008;Rughetti & Festa-Bianchet, 2014), frequency dependence (Jenouvrier et al., 2010), ecological feedbacks (Ransom et al., 2014), and eco-evolutionary feedbacks (Govaert et al., 2019). A sustainable management practice is to set hunting quotas based on population censuses and as a proportion of the population that can be harvested annually (e.g., Andrén et al., 2020). This implies that changes in the phenotypic or genetic composition of the population (e.g., shift in reproductive trait values) or in population management (e.g., shift in target individuals) can affect the proportion of individuals legally protected from hunting and ultimately redirect and exacerbate the hunt toward the remaining, unprotected ones.
In such systems, population composition can feedback on population dynamics through frequency-dependent nonlinearity (Caswell, 2008) between the frequency of protected (unavailable) individuals and the survival rates of available individuals ( Figure 1a).
Here, we took advantage of the long-term and individual-based monitoring of brown bears (Ursus arctos) in Scandinavia to identify the main processes by which hunting selectivity may impact population growth rate. To address this, we developed a state-of-theart multistate dynamic, frequency-dependent matrix population model. Other main features of the model are cause-specific mortality rates, intertwined fates between mother and offspring, and stage transitions dependent on reproductive traits. Using a demographic approach, we analyzed the model with five specific objectives. First, we characterized the brown bear population dynamics under observed conditions, accounting for frequency dependence.
Second, we evaluated the impact of changes in hunting quotas on the population growth rate. Third, assuming potential selection for productivity and shielding traits, we predicted the consequences for population growth rate of changes in four reproductive traits: (1) offspring age at weaning (the probability to wean offspring early, that is, after 1.5 years of maternal care in brown bears), (2) age at first reproduction (the probability to mate at the age of three years and produce cubs at the age of four years), (3) reproductive rate of adult females (an adult female's probability to produce cubs), and (4) number of offspring produced (litter size). We expected a different contribution to population growth rate for each reproductive trait, with adult female reproductive rate having the greatest contribution, as is typically the case for large mammals (Gaillard et al., 1998). Fourth, we evaluated the impact of four hunting regulation scenarios along a gradient of legal protections afforded to members of family groups: (i) no individual is protected, (ii) mothers but not offspring are protected, (iii) mothers and offspring of the year (cubs) are protected, and (iv) entire family groups are protected (i.e., mothers and dependent offspring of any age) (e.g., Swenson et al., 2017). We predicted that an increase in hunting quotas should reduce population growth, but that this reduction would be mitigated by changes in female reproductive traits because such traits have the potential to "shield" females against harvest. Finally, we assessed the interactive effects of simultaneous changes in reproductive traits and hunting quotas under different hunting regulations. We expected that the mitigating effect of changes in reproductive rates on population growth would be exacerbated under higher levels of protection afforded to females as hunting pressure on the remaining available individuals would increase along with the benefits or being in a protected category of individuals.

| Study species, population, and system
The brown bear is a large, nonsocial carnivore widely distributed across Europe, North America, and Asia (Schwartz et al., 2003).
Brown bears typically have a slow life-history strategy (Steyaert et al., 2012); however, reproductive rates vary greatly among populations, contingent mainly on resource availability (Nawaz et al., 2008). Brown bears have been and still are hunted in several populations (Zedrosser et al., 2001), sometimes at very high intensities, such as in Scandinavia .

Our study population of brown bears was located in Dalarna and
Gävleborg counties in south-central Sweden (~61 o N, 14 o E). The mean bear density in our study area was estimated at ~30 bears/1000 km 2 in 2002 . The brown bears' mating season in our study area is in May-July, with a peak in the first week of June (Dahle & Swenson, 2003). During hibernation, which spans from the end of October until the end of April for adult females (Friebe et al., 2001), pregnant females give birth to 1-4 cubs (median = 2 cubs) in January and start lactation while in the den. Survival of cubs-of-the-year (hereafter referred to as cubs) is relatively high, except in the spring, when the risk of sexually selected infanticide (SSI, the killing of unrelated offspring by males to gain access to reproduction with females; Hrdy, 1979) is high Swenson et al., 1997). After litter loss due to SSI in the spring, victimized females soon resume estrus (Steyaert et al., 2014) and can give birth to a new litter the following year . In the absence of complete litter loss, females provide maternal care for at least 1.5 years, and yearling offspring can be weaned after den emergence in their second spring. Alternatively, females can continue maternal care for an additional year, for a total duration of 2.5 years, and two-year-olds are then weaned after den emergence in their third spring (Van de Walle et al., 2021). All offspring from the same family group separate simultaneously or within the same week (Dahle & Swenson, 2003).
Whether or not a female weans her offspring as yearlings or two- year-olds depends on yearling mass in northern Sweden, but such a relationship has not been found in southern Sweden ( Van de Walle et al., 2021).
Bear hunting is allowed throughout the species' range in Sweden and anyone possessing hunting rights in an area and a weapon legal for big game hunting can shoot a bear (Bischof et al., 2008). Annual hunting quotas are set on a county basis and are determined by the national and county wildlife authorities . Successful hunters are required by regulation to report their kill to the authorities and provide the location of the kill, body measurements, sex, and a tooth for age determination (Bischof et al., 2008). There is no limit to the number of bears an individual can shoot, as long as the county-level quota has not been reached. Because there is little incentive for hunters to pass on an opportunity to kill a bear, bear hunting in Sweden is mostly considered as nonselective with regard to age, sex, and size (Bischof et al., 2009), although recent estimates show that hunting may now be slightly biased toward older males and larger individuals (Bischof et al., 2018;Leclerc et al., 2016). However, since 1986, all members of a family group of bears, that is, a female accompanied by dependent offspring of any age, have been afforded legal protection from hunting . By providing a survival advantage to members of family groups, this regulation artificially selects for longer periods of mother-offspring associations (Van de Walle et al., 2018), and any other trait allowing individuals to form and remain in a family group is also expected to be under selection.

| Bear monitoring
The brown bear population in southern Sweden has been monitored using radio-telemetry since 1985. The objective of the monitoring program is to follow female bears, ideally from birth until death . Bears were captured by darting (Dan-Inject, Børkop, Denmark) with an immobilizing drug in the spring, soon after den emergence (Arnemo et al., 2011). Most bears were first captured as yearlings with their mother during the annual spring capture season. For ethical reasons, cubs were not captured. For bears of unknown age when captured, age was determined by analyzing annuli cementum widths of an extracted vestigial premolar (Matson et al., 1993). Captured individuals were measured (e.g., weight, head circumference), identified with a uniquely coded tattoo on the lower lip and a microchip transponder, and their sex was documented. All females were equipped with a VHF transmitter (Telonics, model IMP/400/L HC) implanted in the peritoneal cavity.
From 2003 onward, female bears were also equipped with a GPS collar (GPS Plus; Vectronic Aerospace, Germany), except yearling females due to their rapid growth. After release, radio-marked females were located from the ground or a helicopter a minimum of three times during their active period to assess their reproductive status (solitary, with cubs, or with yearlings) and the number of offspring was counted. Handling of study animals in the monitoring program was approved by the appropriate authorities and ethical committees: Swedish Board of Agriculture (no. 35-846/03, 31-7885/07, 31 11102/12), Uppsala Ethical Committee on Animal Experiments (no. C40/3, C47/9, C7/12), and Swedish Environmental Protection Agency (no. 412-7327-09 Nv). The monitoring program provided information on female reproductive traits (i.e., litter size, age at first reproduction, probability of adult females to produce, and offspring age at weaning). In Sweden, all bears killed legally (e.g., legal hunting, management kills, defense of life and property) must be reported to the management authorities. Death due to other reasons (e.g., natural deaths, vehicle and train collisions, illegal hunting) has also to be reported, although an unknown proportion of mortalities remains undetected (Bischof et al., 2008(Bischof et al., , 2009).

| Model including frequency-dependent hunting mortality
We built a nonlinear matrix population model structured by age and state stages. The model projects the number of individuals within each stage, n, from year t to t + 1 based on the projection matrix A and the vector of parameters θ. The vector of parameters θ is function of the current population vector n (Caswell, 2001;Jenouvrier et al., 2010) because of the dependency of hunting mortality on the frequency (or proportion) of individuals unavailable for hunting in the population. Therefore, the population is projected as follows: If we define hunting quota (q) as the proportion of the total population size (N T,t ) to be harvested each year, the annual probability that an available (unprotected) individual will die from hunting (h t ) is: where N T,t is the sum of available (N a,t ) and unavailable (N u,t ) individuals at time t. It follows that, with a growing proportion of the population protected from hunting, hunting pressure (and thus the risk of being killed) increases for the remaining unprotected individuals. Availability to hunting depends on the state of an individual at the beginning of the hunting season (fall), which can differ from the one at the beginning of the spring. For example, a dependent yearling weaned in the spring becomes available for hunting as an independent yearling in the fall. To account for these seasonal transitions in our estimation of N u,t and N a,t , we built a periodic model with three seasons.

| Periodic model
We start by presenting the annual life cycle graph and then its season decomposition. Annual transitions are estimated assuming postbreeding censuses, that is, between den emergence at year t and den emergence at year t+1 to match with the sampling protocol of the monitoring program. The annual life cycle is based on females assuming a sex ratio of 1:1 within litters and includes 9 stages VAN de WALLe et AL.
( Figure 1b): 5 juvenile and 4 adult stages, based on the age/state in which females are at den emergence at year t: 1. Cubs (0): female cubs born in January of a given year that are dependent on their mother during the entire year.
2. Dependent yearlings (1): after their second hibernation with their mother, female yearlings emerging from the den are still dependent upon their mother.
3. Dependent two-year-olds (2d): 2-year-old females that have hibernated with their mother for a third winter and are still dependent on their mother in early spring of their third year.

5.
Three-year-olds (3): independent 3-year-olds. Note here that survival of cubs is assumed conditional on their own survival and that of their mother. If weaned or if the mother died, yearlings become independent in the fall. Adult females can lose their litter due to weaning or the death of all offspring in the litter. All A2 females will wean their two-year-olds in the spring. All adult females that have lost their litter become solitary in the fall.
Matrix M 2 projects the population from the 8 stages at the beginning of the fall (f) to the 8 stages at the beginning of the winter (Figure 1c), conditional on probability to lose a litter of cubs (l ) or yearlings (y) and stage-specific survival (S i ) rates: In the fall, yearlings that have lost their mother become inde- and thus depends on population structure, following Equation 2 above. N u,t and N a,t are calculated as the number of females within the stages that are protected (i.e., "unavailable") and not protected (i.e., "available") from hunting at the beginning of the fall. Hunting regulations determine the stages which are afforded protection. For instance, in the case of protection of family groups, stages 0, 1d, A0, and A1 in the fall are unavailable, whereas stages 1i, 2, 3, and Al are available. Because the model is females-based, it assumes that an increase in the frequency of unavailable females will spread hunting mortalities over the remaining available females only. Because the hunting quota is fixed based on total population size, and hunters cannot distinguish between males and females at a distance, hunting mortalities should be spread over males as well. Therefore, we included males in our calculations of N u,t and N a,t based on the female population vector, assuming a 1:1 sex ratio in the population. Availability of males was determined from their presumed stage as subadults, but since males do not contribute to parental care in brown bears, all adult males are available for hunting.
Matrix M 3 projects the population from the 8 stages at the beginning of the winter (w) to the 9 stages at the beginning of the spring at year t+1, conditional on the probability of losing a litter of cubs (l ) or yearlings (y), state-specific fecundity (f i ), and survival (S i ) rates: Here again, yearlings that have lost their mother become independent and adult females that have lost their litter become solitary, all other females remain in their previous stage conditional on their survival. Females aged three years and solitary adult females will emerge from their den with cubs with a probability r 3 and r Al and can contribute to the cubs stage the following spring through f 3 and f Al (See Appendix S1: Tables S1 and S2, for detailed mathematical descriptions of the transitions).

| Model parameterization
We used 29 parameters (see Table 1

Population dynamics under observed conditions
We parameterized the frequency-dependent population model with observed parameters (Table 1) to estimate population growth rate under the current hunting regulation protecting members of family groups. In years when population size was estimated, hunting quotas (q, in %) ranged from 3 to 11%, with an average of 5.5% (Figure 2a).
The parameter q was thus set at 5.5% in our simulations. Frequencydependent models eventually converge to an equilibrium population structure (p) and population growth rate (̂ ; Caswell, 2001).
We thus inspected the temporal dynamics of 10 simulations using random initial population vectors over 100 years to identify time at convergence. Once convergence was reached, we derived ̂ as the Then, we assessed how the parameters influence population growth rate using perturbation analysis, adapted for nonlinear models (Caswell, 2008). To account for nonlinearity in our model, we conducted elasticity (i.e., the proportional effect on population growth rate of a proportional change in a parameter) analyses at equilibrium (Caswell, 2008;Jenouvrier et al., 2010 equations 16-17).

Impact of hunting quota
We evaluated the effect of hunting quotas on population growth by re-estimating ̂ over simulated hunting quotas ranging from 0 to 25% by increments of 0.5%. Simulations were performed assuming legal protection of family groups. For each simulated hunting quota, we bootstrapped the procedure by randomly drawing the remaining parameters 1000 times within the parameter normal distribution to obtain a 95% confidence interval (CI) around the mean prediction.

F I G U R E 2 (a)
Temporal changes in Swedish brown bear population size, hunting quotas (number of individuals allowed to be hunted annually), and hunting quotas, q (proportion of total population allowed to be hunted annually). Data were taken from Swenson et al. (2017) for the period 1985-2015, except for q, which was calculated here as hunting quotas (numbers) divided by population size during years when estimates of population sizes were available. The parameter q was calculated only for the years when both population size and hunting quotas (individuals) were available. (b) Population growth rate at equilibrium (̂ ) based on observed vital rates, assuming q = 5.5%. (c) Changes in ̂ under simulated changes in hunting quota, q. Note that here a q value of 0.05 corresponds to a 5% quota. Bootstrapped distributions in (b) and (c) were obtained by randomly drawing parameters within their distribution of estimation 1000 times and under the current hunting regulation of protecting family groups. Only the 95% confidence interval is shown in (c). The dotted lines perpendicular to the x-axis in (b) and the y-axis in (c) indicate the threshold between population decline (̂ < 1) and population growth (̂ > 1). In (c), the dotted lines perpendicular to the y-axis indicate the minimum and maximum hunting quotas leading to stable population (̂ = 1) Year 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999

Effects of changes in reproductive rates
We evaluated the effect of changes in four reproductive traits, that is, a (early weaning probability), r 3 (probability to mate at the age of three years), r Al (probability to produce cubs), and n cub (litter size), on̂ . In our models, we allowed one reproductive parameter to vary at a time, while keeping the other parameters constant. We simulated probabilities (a, r 3 , r Al ) over the range 0-1 by increments of 0.01, whereas counts (n cub ) were simulated between 1 and 4 by increments of 0.01. Simulations were performed with q = 5.5% and assuming legal protection of family groups. As above, we bootstrapped the procedure by randomly drawing the remaining parameters 1000 times within the parameter normal distribution to obtain a 95% confidence interval around the mean prediction.

Effects of hunting regulation
We measured the impact of hunting regulation on population dynamics by recalculating p and ̂ under four hunting regulation scenarios: (i) no individual is protected, (ii) mothers but not offspring are protected, (iii) mothers and cubs are protected, and (iv) entire family groups are protected (i.e., mothers and dependent offspring of any age).
Depending on the scenario considered, availability to hunting in the fall was redefined based on which category of individuals was afforded legal protection. The resulting hunting probability, h, was applied to the available categories of individuals only (Appendix S1: Table S3). To highlight contrasts in hunting regulations, we simulated two hunting quotas: (1) 5%, that is, a hunting quota close to the average observed in our study population, and (2) 20%, that is, an extreme hunting quota.

Interplay between quotas, reproduction, and regulation
Finally, we investigated the interactive effects of hunting quotas and reproductive traits, which would indicate potential for huntingdependent selective gradients on reproductive traits. We started by estimating the elasticities of reproductive parameters over hunting quotas ranging from 0 to 25% to identify linear relationships that would suggest that the potential for reproductive traits to influence population growth rate is hunting-dependent.

| Population dynamics under observed conditions
Over the study period 1985-2015, using mean parameter values, the predicted ̂ was 1.029 (95% CI = [1.011, 1.045]). Of the 1,000 bootstrap iterations, only 0.3% of the simulations showed a population F I G U R E 3 Response of the equilibrium growth rate (̂ ) of the brown bear population in south-central Sweden to changes in (a) the probability to wean cubs at the yearling stage (i.e., early weaning probability), (b) the probability for females aged three years to mate at this age and emerge from their winter den with cubs the following spring (i.e., probability to mate at three years old), (c) the probability for adult solitary females in the winter to emerge from their winter den with cubs the following spring (i.e., probability to produce cubs), and (d) litter size produced at each reproduction event. The lighter area in each panel represents the 95% confidence interval of the observed values for each reproductive trait ( decline (̂ < 1; Figure 2b). The parameters to which ̂ showed the greatest elasticity values were survival rates of cubs at each season (elasticity for S 0,s : 0.070; S 0,f : 0.085; S 0,w : 0.086), litter size (elasticity for n cub : 0.074), probability for adult females to produce cubs (elasticity for r Al : 0.055), and hunting quota (elasticity for q: −0.054; Table 1).

| Impact of hunting quotâ
declined sharply and linearly with increasing hunting quotas ( Figure 2c). Within the hunting quotas interval [0.065, 0.100], the population was stable (̂ = 1.0). Hunting quotas below this and above this range would result in population growth or decline, respectively.

| Effects of changes in reproductive rates
With increasing simulated values of the four female reproductive traits investigated, ̂ increased ( Figure 3). Population growth rate increased linearly with early weaning probability ( Figure 3a) and with the probability to mate at three years of age (Figure 3b), but nonlinearly with the probability of adult females to produce cubs ( Figure 3c) and with litter size (Figure 3d). The population increased (̂ > 1) regardless of early weaning probability and probability to mate at three years old. However, the population was predicted to increase only when probability to produce cubs as adult was ≥0.

| Effects of hunting regulation
Hunting regulation affected both the population structure p and growth rate ̂ (Figure 4). At low (5%) hunting quotas, hunting regulation affected the mean and distribution of ̂ only moderately, but when hunting quotas were high (20%), protecting certain categories of females within the population had a positive effect on̂ . However, ̂ did not increase linearly with the number of categories afforded protection. Indeed, the highest ̂ values were attained when only mothers, and not their offspring, were afforded legal protection ( Figure 4a). Moreover, the effect of hunting regulation on the population stage structure was more apparent under high hunting quotas ( Figure 4b). As protection from hunting was extended to more categories of individuals, there was an increase in the proportion of females forming extended family groups (stages 2d, A1, and A2). The effect was more pronounced under high hunting quotas and with the legal protection of family groups where the proportion of A2 and Al (solitary adult females) were the greatest and lowest, respectively.

| Interplay between quotas, reproduction, and regulation
The elasticities of ̂ to reproductive traits changed with increasing hunting quotas ( Figure 5). With increasing hunting quotas, the elasticity of population growth rate decreased for litter size and early weaning probability, but it increased for the probability to mate at three years old and the probability of producing cubs for adult females. Interestingly, the elasticity of population growth rate to early weaning probability became negative at very high hunting quotas.
This suggests that whereas an increase in early weaning probability would increase population growth rate under low hunting quotas, a similar increase would reduce population growth rate at high hunting quotas.
Hunting quotas, reproductive trait values, and hunting regulations interacted in shaping equilibrium population growth rate, ̂ ( Figure 6). Higher probabilities to produce cubs as adult (r Al ), to mate at three years old (r 3 ), and larger litter sizes would all increase population growth rate at any given hunting quota. When no categories of individuals were protected, ̂ declined sharply and tended to converge to very low values, regardless of changes in reproductive rates. However, when mothers and family groups were protected, increases in r 3 , r Al , and litter size led to a slower decline in̂ . This effect was also detected under the legal protection of mothers and cubs, but only for changes in r Al and r 3 . Across hunting regulations, a high early weaning probability (a) was associated with greater ̂ at low hunting quotas. However, when hunting quotas increased, a higher early weaning probability led to similar ̂ when only mothers were protected and even greater ̂ at high hunting quotas when family groups were protected. See Appendix S1 (Fig. S1) for interaction plots using the full range of reproductive trait values.

| D ISCUSS I ON
Through the removal of individuals, human harvest directly impacts population size and, potentially, the persistence of wild populations (Jackson et al., 2001). Through selection on phenotypic traits, harvest can also indirectly impact population structure and dynamics (Frank et al., 2017;Law, 2000;Milner et al., 2007). Here, we quantified and compared the impact of hunting intensity, hunting regulation, and expected hunting-induced selection on female reproductive traits on the population dynamics of brown bears. We found that the population should grow if hunting quotas are below 10% of the total population size annually under the current regulation F I G U R E 6 Interactive effects of hunting quotas, q, and changes in female brown bear reproductive traits on equilibrium population growth rate (̂ ) under four scenarios of hunting regulation: (1) no individual is protected (column 1; "None"), (2) mothers are protected (column 2; "Mothers"), (3) mothers and their cubs are protected (column 3; "Mothers and cubs"), and (4) mothers and their offspring of any age are protected (column 4; "Family groups"). Hunting quotas were simulated over the range 0-0.25 (note that here a q value of 0.05 corresponds to a 5% quota), changes in the probability to produce cubs as adult (r Al ; first row) probability to mate at three years old (r 3 ; second row), and early weaning probability a; fourth row) were simulated over the range 0-100%, and litter size (third row) was simulated over the range 1-4 λ Hunting quota, q protecting family groups. Among the reproductive traits considered, harvest-induced selection acting on litter size and the probability for adult female brown bears to produce cubs would have the greatest impact on the population growth rate. Considering that sensitivity of population growth rate can be interpreted as a selection gradient (Caswell, 2008), our results show increasing selectivity for producing a litter (conditional on having the possibility to do so) with increasing hunting quotas. We also found that hunting regulations aimed at protecting the female segment of the population are effective, but providing legal protection to dependent offspring may dampen this effect. When family groups are protected, females producing litters are shielded against hunting, but this increases pressure on the other demographic groups that remain available for legal hunting.
Our model predicted an increase in the Swedish brown bear population, with an average 2.9% annual growth rate, which is similar to previous studies on the same population (Gosselin et al., 2015;Van de Walle et al., 2018), but relatively high compared with brown bear populations in North America (Garshelis et al., 2005). In Sweden, the current management objective is to reduce the size of the brown bear population . Our results seem to suggest that this management objective was not met. However, hunting quotas used in our model (average of 5.5%) were estimated based on data from 1985 to 2013 , and considering that hunting quotas and management kills have increased dramatically in recent years , this is probably an underestimation. Therefore, our predicted 2.9% annual growth should be inter- Perhaps the closest form of association among mammals is that between a mother and her dependent young. Their mutual influence on each others' vital rates is apparent; a mother will typically have a positive effect on the survival of her dependent young (Clutton-Brock, 1991;Klug et al., 2012), whereas dependent young may have a negative influence on their mother's fecundity, because in many species females do not enter estrus until their young are weaned (Borries et al., 2014). Additionally, young may also affect their mother's survival (if being part of a family unit can make her more or less vulnerable) and, in the case of extended associations, a mother may influence the fecundity of her offspring (e.g., reproductive suppression; Abbott, 1987). Demographic models should thus account for their intertwined fates (e.g., Hunter et al., 2010). In exploited populations, the association between a mother and her young can be even more important, as hunting regulations often afford protection for mothers, dependent young, or both, and can even reduce the cost of reproduction (e.g., Ericsson, 2001;Krofel et al., 2012;Solberg et al., 2000;Van de Walle et al., 2018). Such regulations can be motivated by population-dynamic considerations, or to limit the effect of hunting on wild animal populations, as well as ethical concerns. In some cases, hunters avoid killing lactating female even though hunting regulation allows it (Rughetti & Festa-Bianchet, 2014 Despite the widespread potential for females and their dependent offspring to benefit from a legal, or ethical protection, seldom are the demographic consequences of those management actions quantitatively assessed. This may be due to the lack of comprehensive population models and the detailed data required to parameterize them. By incorporating cause-and season-specific mortality rates, as well as interdependencies of mother-offspring vital rates, we were able to efficiently predict the outcome of various management decisions (e.g., increase or decrease in hunting quotas and changes in the regulations). Indeed, our results show the effectiveness of legal protection of reproductive females, especially at high hunting intensities, but also the potential for hunting-induced selection on reproductive traits to dampen the effects of an increase in hunting quotas when mothers are protected. Extending the legal protection to all members of family groups has the additional potential to reverse selection gradient on age at weaning.
In ungulates, hunting rate is the most important factor driving the dynamics of exploited populations (Mysterud, 2011;Rughetti & Festa-Bianchet, 2014;Rughetti et al., 2017). That is because, in long-lived species, survival is typically the demographic rate with the highest elasticity and thus with the greatest potential to affect population growth (Gaillard et al., 1998). Here, we show that hunting rate in a large carnivore is also a key driver of population dynamics. Relative to other causes of mortality, population growth rate showed the greatest elasticity for hunting quotas, indicating that hunting has the largest impact on the dynamics of this population, compared with drivers of natural mortality. Our results also suggest that despite the protective effect of female reproductive traits, hunting-induced selection for higher productivity in the population may not suffice to avoid population decline, even under strictly enforced legal protection of family groups. Reaching management goals for long-lived species, such as the brown bear, is thus primarily dependent upon decisions on the level of hunting quotas to be issued, especially considering feedbacks from frequency-dependent hunting mortality.
High levels of extrinsic mortality are expected to induce selection for increased productivity as a compensatory mechanism (Darimont et al., 2009;Law, 2000;Stearns, 1992). In support of this hypothesis, we found that brown bear females producing larger litters at an earlier age and more frequently would lead to increases in population growth rate. Increased productivity through high cub production relative to female size has likely contributed to the persistence of European brown bear populations after centuries of persecution (Zedrosser et al., 2011 We found that hunting quotas have a greater influence on hunted animal populations compared with hunting regulations in our study population. This is also the case for alpine chamois (Rupicapra rupicapra), for which hunting rates affected the population dynamics more strongly than hunting selectivity for nonlactating females (Rughetti & Festa-Bianchet, 2014). The demographic consequences of hunting selectivity in the alpine chamois were only apparent under high simulated levels of hunting and selectivity. However, hunting regulations had measurable (albeit slight) demographic consequences even at low hunting quotas in our brown bear population, probably because the protection of lactating females in our system is strictly enforced and only a handful of mistakenly shot females with cubs have been reported over the last two decades in Sweden (Van de Walle et al., 2018). Our model also revealed that hunting selectivity on female reproductive traits has the potential to dampen the effect of harvest rates, which is in line with recent findings that selection can mitigate the effects of environmental changes on wild animal populations (Urban et al., 2020).
Like all modeling approaches, our predictions are based on a set of assumptions to reduce complexity, or because of the difficulty to estimate certain processes (Caswell, 2001). First, we did not account for potential age differences in vital rates within our adult female stages and used instead parameter estimates averaged for females aged between 6 and 10 years old (Bischof et al., 2018). Young adult females may still divert energy to growth and consequently show reduced reproductive output and survival probabilities due to lifehistory trade-offs (Stearns, 1992). Similarly, older females, through senescence, may show lower reproduction and survival (Kirkwood & Rose, 1991), however, considering that the onset of senescence in brown bears is ~27 years old (Schwartz et al., 2003) and that the average age at death is 4.8 years in Sweden (Bischof et al., 2008), this assumption appears reasonable. Second, we assumed that orphaned cubs would die. The assumption of death of orphaned offspring strongly affects the predicted demographic response of alpine chamois to selective hunting (Rughetti & Festa-Bianchet, 2014). Because we did not capture and equip cubs with radio collars in our study, the fate of orphaned cubs is typically unknown. Some orphaned cubs have been documented to survive , but considering that orphaned offspring can show reduced growth and future survival prospects (Festa-Bianchet et al., 1994), their contribution to the population as adults might be limited (Zedrosser et al., 2013). Third, we assumed males and females as similarly vulnerable to hunting in our calculations of hunting mortality. This assumption is sensible in our study population considering that male and female brown bears are not discernable at a distance for hunters (Bischof et al., 2009), but can be challenged in other populations. Our results would potentially be reinforced in a male-biased harvest system.
Lastly, we assumed density independence as a previous study on the same population did not find relationships between bear density and demographic rates (Bischof et al., 2018). This may mean that the population has not reached stationarity, a phase regulated by density dependence where large populations decrease and small populations increase as resources (food, minerals, space, etc.) availability fluctuates (Coulson, 2020;Coulson et al., 2008). Nonstationarity is a reasonable assumption in the case of heavily hunted populations, but it may not hold for all populations, especially if removal rates are low. Despite those assumptions, our model predicts population growth rates within the range of published estimates for this population of brown bears and others (Garshelis et al., 2005;Gosselin et al., 2015;Kindberg et al., 2011;Van de Walle et al., 2018). As such, it represents a valuable tool for the management of this and other brown bear populations.
As the human impact on wild species is increasing, it is critical to predict the consequence of those changes on population dynamics and evaluate the effects of novel selective pressures on the fate of populations (Lasky et al., 2020). Models that include both human-driven extrinsic mortality, as well as phenotypic trait changes, are therefore useful to explore the effect of multi-trait lifehistory changes in natural populations. In fisheries, harvest-induced evolution is expected to increase the population growth rate in such a way that populations adapted to harvest would support higher levels of harvest (Dunlop et al., 2015;Enberg et al., 2009;Heino et al., 2013). In hunted populations of long-lived species, there is less evidence of harvest-induced evolutionary changes (Pigeon et al., 2016), potentially because they may require a longer period to detect. A critical step would therefore be to expand population models to integrate and quantify the possible feedback of evolutionary changes in life-history traits on population processes (Govaert et al., 2019;Smallegange & Coulson, 2013) especially since evolution can offer a mean to promote population resilience to exploitation (Dunlop et al., 2015). Notwithstanding, what matters the most in terms of demography is to understand the consequences of phenotypic changes, regardless of their cause (Hendry, 2016). As such, phenotypic changes ought to be considered in management and conservation of wild populations, especially as they may offer mechanisms for mitigating the negative consequences of intense harvest.

ACK N OWLED G M ENTS
We are thankful to all the people who have contributed to the data collection. We thank M. Festa-Bianchet for useful comments on

DATA AVA I L A B I L I T Y S TAT E M E N T
The MATLAB codes used for the model construction and analyses are available in the Appendix S1 of this article. The data that support the findings of this study are available upon request in figshare at