Harvest impacts on caribou population dynamics in South West Greenland

We examined the effects of hunting on caribou populations in South West Greenland from year 1999 to 2007. In the Ameralik area a reported average annual harvest of 2950 caribou coincided with a population decline from 31 000 (90% CI: 22 000 44 000) animals in 1999 to 8900 (90% CI: 5800 13 000) in 2007. A survey estimate from 2006 indicates that a suggested target caribou density of 1.2 / km 2 was met. A Bayesian population model estimates the annual replacement for Ameralik at minus 170 individuals (90% CI: -550 460), which indicates that the target density may or may not be maintained even in the total absence of a hunt. For the Qeqertarsuatsiaat area an average annual harvest of 230 caribou appears to have left the density unaffected, remaining steady on target with an abundance of approximately 5000 individuals. The harvest in this area increased from 100 animals in 2000 to 440 in 2006. With an estimated 2007 replacement of 190 (90% CI: -190960) caribou per year the target density may not be maintained in the future unless hunting restrictions are implemented. The density of caribou in Qeqertarsuatsiaat may, however, be maintained over the short term if the emigration of animals from Ameralik into Qeqertarsuatsiaat continues.


Introduction
Caribou (Rangifer tarandus groenlandicus) in West and Northwest Greenland are segregated into approxi¬ mately ten populations that are separated by the Greenland Ice Cap, glaciers and wide fjords.The fjords often penetrate from the sea to the Ice Cap and generally run parallel to the seasonal migra¬ tion of caribou.Most caribou populations in West Greenland are small and relatively isolated with only a small degree of gene flow between them (Jepsen, 1999).Caribou have been abundant, however, in mid-West Greenland for at least the last decade, and the three largest populations are found in what has become known as the North, Central and South regions (Fig. 1).
There are no large predators in the terrestrial ecosystem in West Greenland, and herbivore diver¬ sity is low.In the South region the only permanent resident herbivores that share range with caribou are arctic hares (Lepus arcticus) and ptarmigan (Lagopus Fig. 1.Caribou regions in West Greenland.This paper describes the impact of hunt on the caribou population in the South region. Fig. 2. Details of the South region; area within solid black outline is approx.13 500 km 2 .The hatched line separates the northern Ameralik (8400 km 2 ) and the southern Qeqertarsuatsiaat (5100 km 2 ) areas.
mutus).The South region is located immediately south of Greenland's capital Nuuk, and it contains the third largest caribou population in West Greenland.
In the absence of predators, the limiting factors on the populations include intra-specific competition, range capacity, pathogens and human harvest, as well as stochastic weather events and climatic changes.Since 2001 quotas were increased resulting in greater harvests and may now have become one of the strongest limiting factors on some of the populations.In this paper we examine the impact of recent harvests on the caribou popu¬ lation in the South region, specifi¬ cally in the northern Ameralik and the southern Qeqertarsuatsiaat areas.

Background
The South region, known as caribou hunting regions 4 and 5, has a sea¬ sonally ice-free area of approximately 13 500 km 2 .Steep sided fjords and a rugged alpine terrain characterize this region (Fig. 2).Coastal lowlands are minimal and much of the region has an elevation greater than 200 metres.The northern border is the large Godthåbsfjord that cuts from the sea to the Greenland Ice Cap.The Frederikshåb Isblink, which is a large glacial tongue of the Greenland Ice Cap, forms the southern boundary.To the west is the Davis Strait, and to the east is the Greenland Ice Cap.The region is divided into two areas that provide two sub-populations of caribou; the northern Ameralik (8400 km 2 , hunting region 4) and the southern Qeqertarsuatsiaat (ca.5100 km 2 , hunting region 5).The region's largest human settlement is Nuuk, the capital of Greenland, with about 15 000 inhabitants.It is situated on the region's northern border at the mouth of Godthåbsfjord.There are two small hamlets; Kapisillit, located in the inner reaches of Godthåbsfjord, and Qeqertarsuatsiaat, situated on the seacoast about 100 km south of Nuuk.Development of the region is limited to a hydro power plant at the head of Buksefjord, with a transmission line to Nuuk.Roads are limited to the settlements, and do not link settlements nor penetrate the terrain.

Recent history
Following aerial surveys in the 1990s, the Greenland government's caribou managers thought that caribou populations through¬ out West Greenland were low in number, while the consensus from local hunters indi¬ cated substantially larger populations (Cuyler et al., 2002(Cuyler et al., , 2003(Cuyler et al., , 2005(Cuyler et al., , 2007;;Cuyler, 2007).New improved abundance surveys beginning in 2000 supported the hunters' knowledge.
The 2001 aerial survey resulted in a caribou estimate of 32 000 (CV: 18%, Table 1) for the Ameralik portion of the South region (Fig. 2).This was seven times greater than the estimate of 4500 from 1996 (Ydemann & Pedersen, 1999), and caribou density was almost four times greater than a proposed target caribou density of 1.2 / km 2 (Kingsley & Cuyler, 2002;Cuyler et al., 2003Cuyler et al., , 2005;;Cuyler, 2007).Following the large numbers of caribou, hunters observed habitat degradation; car¬ ibou-lichen heaths had become overgrazed where they had been deep and plentiful, and the once lush expanses of crowberries (Empetrum nigrum) were trampled (Cuyler et al., 2007).In the years following 2001, the government's caribou managers decided to attempt to reduce West Greenland caribou populations towards the target density.Har¬ vest levels were raised by greatly increased caribou quotas in 2000, 2001 and 2002, and by may favour the preservation of vegetation quantity, extending the hunting season.Free/unlimited harvest quality and availability.If true, the latter would over the extended season was permitted in 2003, and ultimately benefit the body condition, health and this raised harvest levels further by permitting sev-productivity of West Greenland caribou, and may eral thousand sport hunters to partake in the caribou provide the foundation for sustainable harvests.hunt, which until then had been monopolized by With hunting seasons lengthened and quotas commercial hunters.
increased or unlimited, the reported total harvest The target caribou density of 1.2 / km 2 was based increased (Fig. 3, Table 2) to a maximum of about on studies of carrying capacity in North America 7000 caribou in 2003, and then fell to between and Finland (Haber &Walters, 1980, Helle et al., 3000 and4000 annually (Greenland Self Rule, 1990) and correlations between observed densities unpubl.).Since harvest reporting was voluntary, the and changes in caribou productivity, dispersal or the actual harvests may have been far larger than the condition of the range (Kingsley & Cuyler, 2002) their caribou harvest.Meanwhile, the Greenland commercial hunter organization (KNAPK) stated that 80% of the commercial harvest, which was sold at public market in Nuuk, came from the Ameralik area (Cuyler et al., 2007).The trend of the total harvest over time in the Ameralik area shows an initial increase, coinciding with increasing quotas, followed by a peak harvest of 4700 caribou in 2003, which matched the first year of unlimited harvest.Although unlimited harvests continued in the years following, the caribou harvest from Ameralik declined.This may have reflected a decline either in total caribou abundance, or in the number of animals in the areas that were accessible to hunters.In contrast, the caribou harvest from Qeqertarsuatsiaat increased from 100 caribou in 2000 to 440 in 2006.
Prior to 2000, the reported harvest was comprised of 90% males (Loison et al., 2000).Following recommendations from the Greenland Institute of Natural Resources to harvest more females, however, more females were reported harvested than males in Ameralik (Fig. 4, Table 3) in all but one year.Mean¬ while sex ratios in the Qeqertarsuatsiaat harvest were almost even.
The 2006 aerial survey of the Ameralik area resulted in an estimated caribou population of approximately 9680 (90% CI: 6515 -13 147), which constituted a 70% decrease in abundance from 2001.This indicated that the target density had been reached.Local commercial hunters agreed with the decrease in abundance, and reported that there had been no increase in natural mortality based on their own observations (Cuyler et al., 2007).No significant change could be shown in Qeqertarsuatsiaat, which remained at just over 5000 caribou and a density of 1 / km 2 (Cuyler et al., 2007).

Bayesian population dynamics model
We used Bayesian statistical analysis to examine the recent population dynamics of caribou in the Ameralik and Qeqertarsuatsiaat areas.We fitted an age and sex structured population dynamic model to the abundance data, subtracting the annual esti¬ mates of the sex and age-specific harvest from the population.As the time period for which we have area specific harvest data is relatively short, from 2000 to 2006, we applied an exponential popula¬ tion dynamic model.Although we have only two abundance estimates, we allow ourselves to fit the model to the abundance data in order to maintain a best first estimate of the current production levels.
The exponential model implies that we have no con¬ trol over regulating factors and therefore we cannot make long-term predictions.But under all circum¬ stances, long-term predictions in caribou are often problematic because caribou may have fluctuating or cyclic dynamics.This implies that traditional den¬ sity regulated models cannot likely describe caribou population dynamics: relevant models need also to consider delayed density dependent factors that may have a strong influence on the dynamics.A clearer understanding of such processes in West Greenland caribou requires harvest and abundance data for a much longer period of time.
During our analysis we focus on the effects of the recent increase in harvest on caribou population dynamics for the South region.During the 2000¬ 2006 period, the reported caribou harvest from the South region increased from under 3000 to a maxi¬ mum of about 6500 individuals in 2003, and then declined to about 3,400 individuals in 2006.We attempt to determine whether observed abundance changes are the direct result of harvest.Further, we examine whether the target density was reached, and attempt to estimate what harvest levels are needed in the future to maintain or approach the target.

Data
From 2001 to 2006 total abundance was estimated twice by aerial strip surveys (Table 1).Both sur¬ veys used the same design and included identical transects.
For 2003, 2004 and 2005 the annual harvest of caribou (Table 2) from each herd in the South region was estimated using details from hunter harvest reports, which include information of among other things, the location, sex, age category and rump fat depth of each animal.Location was often missing from the hunter reports and the Greenland Self Rule government's total annual harvest data is only avail¬ able per municipality and not per caribou population.We estimated the harvest from the Ameralik and Qeqertarsuatsiaat herds by comparing each munici¬ pality's total annual harvest from all six caribou populations in West Greenland to the relative per¬ centage of harvest from each population contained in the hunter reports with location data.The individual estimates of caribou killed per population by each municipality were then summed to obtain a total annual harvest from each population.Detailed har¬ vest databases, however, were not available for 2000, 2001 and 2002.To obtain estimates of annual harvest for these years, a population's average percentage of the total annual harvest for the years 2003, 2004 and 2005 was applied against the total harvest for 2000, Table 1.Estimates of total caribou abundance (N) for the Southern population divided into Ameralik and Qeqertarsuatsiaat.N includes all age classes, and cv is given in %. (Cuyler et al., 2003(Cuyler et al., , 2007;;Cuyler, 2007).The age-structure for the male and female harvest for each area (Table 3) was estimated from hunter reports.These reports were used to separate the har¬ vest into calves, juveniles and adults.For each area the age-structure in the harvest was estimated from the average age structure over the whole period.The adult harvest included all animals aged over 3-years, under the assumption of no age-class selectivity, a stable age-structure, an annual adult survival rate of 0.91.

Population dynamic model
The applied population dynamics is exponential with constant survival and fecundity rates in an age-and sex-aggregated model.
Let the number of animals in age classes larger than zero be

Statistical methods
The population dynamic models were estimated from the abundance data by projecting the population given the historical harvests, with the initial abundance drawn from a prior distribution of the abundance in the first year of the iteration (assuming a stable age-structure given the fecundity and har¬ vest of that year).A Bayesian statistical method (e.g., Berger, 1985;Press, 1989) was used, and posterior estimates of the model parameters and other manage¬ ment related outputs were calculated.This implied an integration of the product between a prior distri¬ bution for each parameter and a likelihood function that links the probability of the data to the different parameterisations of the model.

Prior distributions
The values and prior ranges of the different param¬ eters for all the assessment models are listed in Table 4. Annual survival rates of 0.90 and 0.92 have previously been applied to adult caribou in West Greenland (Bergerud, 1980;Cuyler & Østergaard, 2005).We applied a Beta (a : 98.5; b : 9.74) prior with mean 0.91 and variance 0.00075 to adult survival, with the choice of variance being rather arbitrary to capture uncertainty in survival beyond the point estimate of 0.91.Fifteen late winter herd-structure counts in West Greenland observed calf percentages per female between 0.16 and 0.77, with an average value of 0.47 calves per female and a variance of 0.036.We applied a Beta (a :2.72; b :3.12) prior with mean 0.47 and vari¬ ance 0.036 to annual reproduction for adult females, which implies that the majority of first year mortality is incorporated into our estimate of reproduction.
As our estimate of annual reproduction is based on the late winter calf percentage, the majority of first year mortality is incorporated into our reproduction estimate.Our estimate of first year survival should thus be correspondingly small.However, having no estimate of calf survival for the remaining time of the first year we applied the adult survival prior also to caribou in age-class zero.
Several female caribou in West Greenland in 1996¬ 97 had their first calf in their second summer, being less than two years old (Cuyler & Østergaard, 2005).Normally, however, female caribou have their first calf in their third summer (Dauphine, 1976, Adams & Dale, 1998, Russell & McNeil, 2005).We applied a uniform prior for age of first reproduction from one to three years capturing this range for reproductive maturity.Apart from the distributions given in Table 4, for each randomly selected parameter set, the upper bound on the juvenile survival rate was always set to be less than or equal to the randomly selected value for the adult survival rate.

Bayesian integration
The Bayesian integration was obtained by the sampling-importance-resampling routine (Berger 1985;Rubin 1988), where n1 random parameterisations # (1 < i < n1) are sampled from an importance function h( 9).This function is a probability distribution func¬ tion from which a large number, n1, of independent draws of 9 can be taken.h(9) shall generally be as close as possible to the posterior, however, the tails of h( 9) must be no thinner (less dense) than the tails of the posterior (Oh & Berger, 1992).For each drawn parameter set 9 i the population was projected from the first year with a harvest estimate to the present.For each draw an importance weight, or ratio, was then calculated where L(9') is the likelihood given the data, and h(9i ) and p(9 i ) are the importance and prior functions evaluated at 9 i .In the present study the importance function is set to the joint prior, so that the impor¬ tance weight is given simply by the likelihood.The n1 parameter sets were then re-sampled n2 times with replacement, with the sampling probability of the ith parameter set being This generates a random sample of the posterior distribution of size n2.The resample of the posterior distribution was set to n2 = 5000, and the n1 sample from the joint prior being 1 000 000.
The method of de la Mare (1986) was used to calcu¬ late the likelihood L under the assumption that observation errors were log-normally distributed (Buckland 1992)

L =n M-^^V
where Nt is the projected and Nt i the point estimate of the observed total abundance at time t, and cvt is the coefficient of variation of the abundance estimate at time .
If the importance function is adequately specified, the mean of the importance sample for each param¬ eter should approach the mean from the true poste¬ rior distribution, given a sufficiently large sample.To illustrate whether the sampled posterior quanti¬ ties can be assumed to be representative of the true posterior distribution, convergence diagnostics were calculated.One such diagnostic is the maximum importance weight of a parameter set relative to the total summed importance weight over all n1 draws, another is the total number of unique parameter sets in the resample of n2 parameter sets, and a third is the maximum number of occurrences of a unique parameter set in the resample.

Posterior distributions
The maximum importance weight of a parameter set relative to the total sum of importance weights for all drawn parameter sets was essentially zero for all assessments.The number of unique parameter sets in a resample of 5000 parameter sets was greater than 4805 for all models, while the maximum occurrences of a unique parameter set in the resample across all models was 4. The model specific statistics are given in Table 5.The posterior estimates and their 90% credibility intervals are given in Table 6.
Rangifer, Special Issue No. 19, 2011 Table 5. Sampling statistics for the Bayesian runs of the different assessments models.Sample is the number of draws from the importance function; Weight is the maximum importance weight of a draw relative to the total importance weight of all draws (given in percent); Unique is the number of unique parameter sets in the resa¬ mple of 5000 parameter sets; and Max is the maximum occurrence of a unique parameter set in the resample.The realised prior and posterior distributions of the population dynamic parameters for the Southern population, and Ameralik and Qeqertarsuatsiaat sub-stocks showed that in all cases, updating of the prior to the posterior indicated a smaller population dynamic growth rate than assumed by the joint prior.For all cases most of the updating was towards lower reproduction, while for the survival rates, the poste¬ rior distributions remain closer to their prior.These differences may reflect the constraints of the model, more than they reflect the true values for the parameters in West Greenland caribou.
While the abundance data and the estimated pro¬ jections showed a marked decline in abundance from 2000 to 2006 for Ameralik, the abundance remained relatively stable for Qeqertarsuatsiaat (Figs. 5,6).

Southern population
The Southern population is estimated to have declined from 37 000 (90% CI: 27 000 -51 000) individuals in 1999 to 13 000 (90% CI: 10 000 -17 000) individuals in 2007.The latter abundance corresponds to a density of 1.0 / km 2 , which is below the target density.For the hypothetical case of no hunting, the abundance should have remained relatively stable at about 37 000 individuals, given an estimated growth rate of 0% (90% CI: -7% -6%) per year.

Ameralik
For Ameralik the abundance is estimated to have declined from 31 000 (90% CI: 22 000 -44 000) individuals in 1999 to 8900 (90% CI: 5800 -13 000) individuals in 2007.The latter abundance corresponds to a density of 1.1 / km 2 , which is close to the target caribou density of 1.2 / km 2 .Fig. 5.The point estimates of abundance and the projected median and 95% CI for Ameralik.
Fig. 6.The point estimates of abundance and the projected median and 95% CI for Qeqertarsuatsiaat.

Qeqertarsuatsiaat
For Qeqertarsuatsiaat the abundance is estimated to have increased slightly from 5200 (90% CI: 2900 -9300) individuals in 1999 to 5300 (90% CI: 3200 -8700) individuals in 2007, with an estimated growth rate of 4% (90% CI: -5% -14%) per year in the absence of harvest.The 2007 abundance corresponds to a density of 1.0 / km 2 , which is just below the target of 1.2 / km 2 .

Discussion
We conclude that the strong decline in the Southern population is the result of a hunt with an average harvest of almost 3200 caribou per year since 2000.Should the annual harvest remain at the 2006 level of 1900 caribou, we expect from the model that this population will continue to decline to an abundance of 2790 (90% CI: 425 -9164) individuals by 2012.This corresponds to a density of only 0.2 / km 2 , which is considerably below the target of 1.2 / km 2 .The current negative replacement yield of -210 (90% CI: -760 -620) caribou per year, suggests that the target density may or may not be maintained in the com¬ plete absence of a hunt.
It is the steep decline of the caribou population in Ameralik that is caus¬ ing the decline for the overall popula¬ tion.For the hypothetical case where no hunting is allowed, the abundance in Ameralik would have remained rel¬ atively stable with an estimated expo¬ nential growth rate just below zero [0% (90% CI: -8% -8%) per year].Therefore we conclude that the steep decline of caribou in the Ameralik area is the result of hunting with an average harvest rate of 2950 individu¬ als per year in this area since 2000.
Should the annual harvest remain at the 2006 level of 1600 caribou we expect that Ameralik will continue to decline to an abundance of only 860 (90% CI: 0 -6123) individuals by 2012.This gives a density of only 0.10 / km 2 , which is far below the target.The negative 2007 replacement yield of -170 (90% CI: -550 -460) caribou per year, suggests that this sub-stock may or may not be able to maintain the target density in the complete absence of a hunt.
In contrast, replacement was posi¬ tive for the Qeqertarsuatsiaat compo¬ nent of the Southern population.We can conclude that the average harvest of 230 caribou per year in this area since 2000 has been close to the yearly recruitment level, thus maintaining the caribou den¬ sity close to the target of 1.2 / km 2 .
The 2007 replacement for Qeqertarsuatsiaat is estimated at 190 (90% CI: -190 -960) caribou per year.However, the harvest of caribou from this area has increased over the 2000-2006 period.The 2006 harvest of 440 individuals exceeds the point estimate of current replacement yield and, thus, given the scenario of stable harvest at present levels, we would predict the abundance to decline to 3900 (90% CI: 800 -12 000) individuals by 2012, which would correspond to a density of 0.8 caribou per square kilometre.Some harvest restrictions in this area may be required to align future harvests with the replace¬ ment yield of 190 caribou per year.
In brief, the results indicate that 1) the abundance of caribou in Ameralik declined from 2001 to 2006, 2) that the target density (1.2 / km 2 ) was reached, and 3) that the population dynamic growth rate between 2001 and 2006 was estimated at approximately zero if no hunting had occurred.It appears that the increased harvest, with annual harvest rates between 3000 to 4500 individuals, caused the strong decline in abundance over this five-year period.Should harvest levels remain constant, then a further population decline may occur.
Following the late winter aerial population survey in 2006, harvest restrictions were implemented for the autumn 2006 hunting season.The season was reduced from 14 weeks to five, but unlimited har¬ vests remained.
Although harvest is the likely cause of this great reduction in the caribou population of Ameralik from 2000 to 2006, questions still remain.Will the Ameralik caribou population continue to decline, and will the Qeqertarsuatsiaat caribou population remain stable?Although Qeqertarsuatsiaat caribou abun¬ dance appeared stable despite increasing harvest from 2000 to 2006, this may have been primarily due to an increased emigration from Ameralik during the same period (Cuyler et al., 2007).The observed movement of Ameralik caribou southward, expanding into Qeqertarsuatsiaat, is supported by local hunters.
The 2007 and 2008 harvest data are not yet available.Despite the model projections presented here, during the autumn 2008 local hunters (pers.comm.)subjectively observed that Ameralik caribou were once again abundant but skinny, and that there were many cow-calf pairs, while Qeqertarsuatsiaat caribou were fat.The next aerial survey of the South region will be March 2011.Given our incomplete under¬ standing, there is room for uncertainty and debate regarding recruitment and future abundance of these two caribou populations.

Table 4 .
Prior distributions for the different assessment models.The list of parameters: sd is adult survival, sjuv juvenile survival, b the birth rate, the age of reproductive maturity, $ the fraction of females at birth, and N0 the abundance in the first year of the iteration (given in thousands).The type of probability distribution is given by superscripts; u = uniform, i = discrete uniform, b = beta, and p a parameter with fixed value.The first number of an entry in the table is the minimum value if pd = u or pd= i, the mode if pd = b, and a fixed parameter value if pd = p.The second number is the maximum value if pd = u or pd = i, and the mean if pd = b.

Table 2 .
The estimated total annual harvest of male and female caribou in Ameralik and Qeqertarsuatsiaat.

Table 3 .
The average age and sex-structure of the caribou harvested from the Southern population.
is the harvest of males/females of age a during year t, with the relative age distribution of the harvests being sex specific, the same in all years, and given by Table3, provided that the harvest will not exceed the abundance in any age class.If instead the harvest exceeds the abundance in an age class, harvest in that age class is set to the abundance of that class.This distribution of harvests is continued until it is possi¬ ble to distribute the remaining harvests in accordance with the age-structure in Table3.Let the annual survival rate sa of animals of age a be Mt ais the number of mature females in age class a at the start of year t, defined as -1 N Rangifer, Special Issue No. 19, 2011 where sa is age specific annual survival, m/f Nl is the number of males/females of age a at the start of year t, x is the lumped age-class (x = 5),where sum is the survival rate of 'juveniles' given the survival of their mother, sad is the survival rate for adults, and ad = 1 is the greatest age at which the 'juvenile' survival rate applies.The number of births at the start of year t, Bt, is a=a m where am is the age of reproductive maturity, and Bt,a, the number of births in age class a, is bM[ a where b is the fecundity rate for mature females, and

Table 6 .
Parameter estimates for the different assessment models.The estimates are given by the median and the 90% credibility intervals of the posterior distributions.r is the exponential growth rate, N the abundance in 1999, NT the abundance in 2007, and ry the replacement yield in 2007.