2018 status muskoxen, Maniitsoq & Sisimiut, West Greenland. Technical Report No. 119

1 2018 Status Muskoxen, Maniitsoq & Sisimiut, West Greenland Technical Report No. 119, 2022 Pinngortitaleriffik – Greenland Institute of Natural Resources 2 Title: 2018 status muskoxen, Maniitsoq & Sisimiut, West Greenland Authors: Christine Cuyler1, Tiago A. Marques2, Iúri J.F. Correia3, Aslak Jensen4, Peter Hegelund1 and Jukka Wagnholt5 1 Pinngortitaleriffik – Greenland Institute of Natural Resources, P.O. Box 570, 3900 Nuuk, Greenland 2 CREEM University of St Andrews, School of Mathematics and Statistics, Scotland 3 University of Lisbon, Faculty of Sciences, Portugal 4 Solviaq 15, 3900 Nuuk, Greenland 5 Tusass, P.O. Box 1002, 3900 Nuuk, Greenland Series: Technical Report No. 119, 2022 Date of publication: 16 February 2022 Publisher: Pinngortitaleriffik – Greenland Institute of Natural Resources Financial support: Government of Greenland and Pinngortitaleriffik – Greenland Institute of Natural Resources Cover photo: Christine Cuyler: Group of 23 muskoxen in the Sisimiut sub-area, between town of Kangerlussuaq and the Isortoq river. ISBN: 9788797297735 ISSN: 1397-3657 EAN: 9788797297728 Cited as: Cuyler, C., Marques, T.A., Correia, I.J.F., Jensen, A., Hegelund, P. & Wagnholt, J. 2022. 2018 status muskoxen, Maniitsoq & Sisimiut, West Greenland. Pinngortitaleriffik – Greenland Institute of Natural Resources. Technical Report No. 119. 113 pp. Contact address: The report is only available in electronic format. PDF-file copies can be downloaded at this homepage: https://natur.gl/forskning/rapporter/ Pinngortitaleriffik – Greenland Institute of Natural Resources P.O. Box 570 3900 Nuuk Greenland Phone: +299 36 12 00 E-mail: [email protected] www.natur.gl 3 2018 status muskoxen, Maniitsoq & Sisimiut, West Greenland By Christine Cuyler1, Tiago A. Marques2, Iúri J.F. Correia3, Aslak Jensen4, Peter Hegelund1 and Jukka Wagnholt5 1 Pinngortitaleriffik – Greenland Institute of Natural Resources, P.O. Box 570, 3900 Nuuk, Greenland 2CREEM University of St Andrews, School of Mathematics and Statistics, Scotland 3CEAUL, University of Lisbon, Faculty of Sciences, Portugal 4Solviaq 15, 3900 Nuuk, Greenland 5 Tusass, P.O. Box 1002, 3900 Nuuk, Greenland Technical Report No. 119, 2022 Pinngortitaleriffik – Greenland Institute of Natural Resources 4 [Empty page] 5 Table of Contents Summary (English) .................................................................... 8 Eqikkaaneq (kalaallisut) ....................................................... 90 Resumé (dansk) ........................................................................ 13 Introduction ........................................................................... 136 Methods .................................................................................... 19 Results ..................................................................................... 266 Discussion .............................................................................. 439 Acknowledgements ............................................................... 566 Literature cited ...................................................................... 566 Figures 1. North region surveyed by helicopter in 2018, … Page 16 2. Reported harvest, commercial, recreational, and combined, for the … Page 17 3. Area covered by 2018 caribou survey of the North region (23,303 km2) … Page 22 4. The 19 line transects used in the 2018 caribou survey of the North region… Page 23 5. Flow diagram of the modelling process during the analysis, from DS… Page 25 6. Muskox survey 2018: exploratory analysis plots… Page 30 7. Muskox survey 2018: exploratory analysis: group size distribution… Page 31 8. Muskox survey 2018: exploratory analysis of non-truncated… Page 32 9. Observer effect: histograms: detected distances for the two observers Page 33 10. Histogram of the different binning options for the muskox distance data… Page 33 11. The detected distances with the estimated detection function overlaid… Page 35 12. Estimated probabilities of detection, given non-truncated and truncated … Page 36 13. Muskox encounter rate (groups per km) per line transect, … Page 36 14. Relative distribution of muskox numbers along the line transects… Page 37 15. Distance sampling design-based muskox density estimates with CIs… Page 39 16. Fitted DSM plots corresponding to the fitted smooth of elevation… Page 40 17. Predicted muskox density across silhouette map of entire North region. Page 42 18. Fitted DSM variability silhouette map of entire North region. Uncertainty … Page 42 19. Virtually snow-free terrain of the inner Uniiviit Fjord, Angujaartorfiup … Page 44 20. Virtually snow-free lowland (foreground <200 m) of Angujaartorfiup … Page 45 21. Large expanse of virtually snow-free terrain in Angujaartorfiup sub-area … Page 45 22. Index of Maniitsoq muskox abundance in the Angujaartorfiup sub-area … Page 53 6 23. Plot sampling grid example of total area A divided into smaller plots… Page 63 24. Example of a patch of tundra with the transect in the middle… Page 66 25. Half-normal (top row) and hazard-rate (bottom row) detection functions… Page 68 26. Possible shapes for the detection function when cosine adjustments are…. Page 69 27. A good model for the detection function should have a shoulder… Page 72 28. Group of 30 muskoxen, which may have included 2-3 calves… Page 79 29. Group of 38 muskoxen, which may have included 3 calves… Page 79 30. Group of 23 muskoxen, which may have included 1-2 calves… Page 80 31. Group of 38 muskoxen, which may have included 7-8 calves… Page 80 32. Map illustrating elevation in the surveyed area. Legend colour codes… Page 81 33. Map illustrating aspect in the surveyed area. Legend colour codes… Page 81 34. Map illustrating slope in the surveyed area. Legend colour codes… Page 82 35. View from surface of lake ice during a pause between survey lines... Page 83 36. View as helicopter flew over surface of lake ice, illustrating the typical… Page 83 37. View while flying line transect 17, illustrating sparse snow-cover… Page 84 38 Overview of the north side of Tasersuaq Lake (large flat surface on left… Page 84 39. Overviews illustrating flat-light (no shadows) and the windswept terrain… Page 85 40. Minimal March snow cover typical of lowland valleys in Angujaartorfiup… Page 86 41. Observation platform for aerial survey: AS350 helicopter… Page 86 42. Map for the 2010-2019 period illustrating the entire Maniitsoq muskox … Page 88 43. Map illustrating the entire Maniitsoq muskox harvest management area … Page 89 Tables 1. Current naming of region, municipality, muskox population, harvest … Page 19 2. Summary of unprocessed results: Survey of the North region, … Page 26 3. Behavioural reaction of muskoxen to helicopter fly-by, as per … Page 27 4. Approximate elevations for muskox groups observed … Page 29 5. Summary of the coefficient characteristics of the GLM … Page 34 6. Model comparison across three Conventional Distance Sampling models… Page 38 7 Encounter rate estimates per sub-area for muskox groups … Page 38 8. March 2018, the Distance Sampling design-based muskox abundance … Page 38 9. March 2018, the Distance Sampling design-based muskox density … Page 38 10. GAM model summary table relative to smooth terms for covariates…Tweedie Page 39 11. GAM model summary table relative to smooth terms for covariates…NegBin Page 40 7 12. Comparison of DS design-based and GAM/DSM analyses, for entire North … Page 41 13. Maniitsoq and Sisimiut muskoxen 2000-2014: winter minimum counts, ... Page 61 14. Rough densities for Maniitsoq muskoxen in the Angujaartorfiup sub-area. Page 62 15. Commonly used key functions & series expansions for the detection function… Page 67 16. Muskox hunting seasons and quotas for Sisimiut and Maniitsoq … Page 87 17. Reported number of muskoxen harvested by commercial and recreational … Page 90 Appendices 1. Minimum counts 2000-2014 for Maniitsoq and Sisimiut muskoxen Page 51 2. Statistical methods behind Distance Sampling Page 63 3. Statistical methods behind GAM and DSM modelling Page 75 4. Distance Sampling Assumptions – short summary Page 78 5. Sisimiut sub-area inland, within 40 km of the Ice Cap, March 2018: … Page 79 6. Maps illustrating elevation, aspect, and slope in the North region Page 81 7. Angujaartorfiup sub-area, March 2018: Photographs survey conditions … Page 83 8. Hunting seasons, quotas & reported harvests, primarily for Maniitsoq … Page 87 9. 2020 Supplementary materials to biological advice … Page 91 10. Recommendations for improving future surveys for muskoxen Page 112 Raw data may be accessed by contacting the Pinngortitaleriffik – Greenland Institute of Natural Resources, Department of Mammals and Birds, Nuuk, Greenland. 8 Summary This report presents results from the first systematic aerial survey of muskoxen in southwest Greenland. The survey involved the North region (66°-68°N), which was divided into three sub-areas (Angujaartorfiup, Sisimiut, Sisimiut South). The survey occurred early March 2018 and provides the first estimate for muskox abundance in two muskox harvest management areas, Maniitsoq (66°-67°N) and Sisimiut (67°-68°N). The behavior of most muskox groups was unaffected by the helicopter fly-by at 40 m altitude. Regardless of group size, 77% of groups simply stood still. Detecting stationary groups is clearly essential for accurate population estimates. Sisimiut muskoxen were most often observed at elevations of ca. 200 m, which is typical for this species as they prefer to forage in lowland elevations even in winter. In sharp contrast, Maniitsoq muskoxen used 700-800 m elevations, despite their documented year-round preference for lowlands <400 m. Further, at the time of the survey, Maniitsoq muskoxen were clumped into two ‘hotspots’ almost inaccessible by motor vehicle and relatively far from human habitation. Since 84% of all muskox harvest (commercial and recreational), as well as most trophy hunting and qiviut (muskox inner wool) production in Greenland are taken from the Maniitsoq muskox population, these activities may have a role in the disruption of normal lowland distribution of the Maniitsoq muskoxen in winter. From a subset of the Maniitsoq data, the calf (age <1-year) percentage was ascertained ca. 18% for Maniitsoq muskoxen. Considering the absence of large predators, the current value is considered low. Factors involved may include density-dependent issues associated with the Maniitsoq muskoxen now foraging on high elevation suboptimal habitat in winter. Since this coincides with late gestation for muskoxen, calf production could be negatively affected. Whether current calf percentage is sufficient to support population size stability or growth is debatable. Population size & density estimates Conventional Distance Sampling (DS) design-based methods and analyses, as well as Generalized Additive Models/Density Surface Modelling (GAM/DSM) based analyses were applied to the dataset to obtain estimates of muskox population size and density. Most muskoxen occurred in the Maniitsoq muskox harvest management area (surveyed Angujaartorfiup sub-area). Far fewer muskoxen inhabited the Sisimiut muskox harvest management area (surveyed Sisimiut sub-area), and zero muskoxen were observed in the surveyed Sisimiut South sub-area. Densities from the GAM/DSM model-based analysis supported that in early March, muskoxen strongly preferred 9 Southwest facing slopes and specifically for Maniitsoq muskoxen the highest densities coincided with high elevations. For the entire North region, the DS design-based March 2018 muskox population abundance was estimated at 21,746 muskoxen (95% CI: 11,061–42,751; CV = 28.5%; SE: 6,194), with a density of 1.1 muskoxen/km2 (95% CI: 0.559–2.160; CV = 28.5%; SE = 0.313). Alternately for the entire North region, the GAM/DSM model-based March 2018 muskox population abundance was estimated 23,256 muskoxen (95% CI: 18,102– 29,877; CV = 11.36%), with mean density of 3.69 muskoxen/km2 (95% CI: 2.87-4.74). DS estimates as per specific sub-area: Maniitsoq muskoxen (Angujaartorfiup sub-area) : The DS design-based estimate was ca. 18,906 muskoxen (95% CI: 8,726–40,960; CV = 0.315; SE = 5,948), with a density of ca. 2.6 muskoxen/km2 (95% CI: 1.223–5.742; CV = 0.315; SE = 0.834). Sisimiut muskoxen (Sisimiut sub-area) : The DS design-based estimate was ca. 2,840 muskoxen (95% CI: 662–12,178; CV = 0.568; SE = 1,613), with a density of ca. 0.22 muskoxen/km 2 (95% CI: 0.052–0.962; CV = 0.568; SE = 0.127). The population size estimates from the two approaches were similar since the 95%CIs overlap. Despite the good survey coverage (10.6%), the high variability within the dataset was responsible for substantial uncertainty in the DS estimates and less so in the GAM/DSM. Regardless, calculating the CV for probability of detection of muskoxen permitted comparison of the two approaches. The DS probability of muskox detection had a CV of 5.98%, which was better than the CV of 11.36% for the GAM/DSM. Thus, for the 2018 survey for muskoxen, we recommend using the DS design-based abundance and density estimates when making management decisions. Alone, the 2018 estimate cannot indicate population trend. That requires at least two additional aerial survey estimate points using similar methods. Meanwhile, past counts, densities, harvests, and calf percentages do not suggest recent population growth but possibly a decline prior to the 2018 survey. Regardless, specifically the 2018 population size estimate for Maniitsoq muskoxen is larger than previous estimates for populations anywhere in Greenland. Also, Maniitsoq muskox density is much higher than elsewhere in the Arctic. This could increase exposure of individuals to infectious pathogens. Given possible density-dependent influences acting at current population size, population growth may not be advisable. Whether the 2018 population size and density for the Maniitsoq muskox harvest management area are within the current herbivore carrying-capacity of the pasture/range remains to be seen. 16 Introduction Muskoxen (Ovibos moschatus Zimmermann) are native only to the north and northeast part of Greenland. Nevertheless, there are currently several muskox populations in west and northwest Greenland. These resulted from translocations, excepting one (Inglefield Land), which is a combination of native and translocated muskoxen. Translocations began in the 1960s when 27 muskoxen live-captured in East Greenland were transported and introduced to Kangerlussuaq (Søndre Strømfjord) in West Greenland (Fig. 1). These animals became firmly established and until 2015 Pinngortitaleriffik – Greenland Institute of Natural Resources (GINR) publications referred to these as the Kangerlussuaq muskox population. Today, the Government of Greenland designates these as the Maniitsoq population. So named because they inhabit what was once the Maniitsoq municipality. From 1986 to 1991, the new Maniitsoq population became the source for several further muskox translocations along the west and northwest coast. Further, by the early 2000’s some Maniitsoq animals expanded northward into what was then the Sisimiut municipality. Although not truly another population, harvest was managed separately in accordance with the then two separate municipal jurisdictions, and it became common to regard them as two populations, Maniitsoq and Sisimiut. In 2009, those two municipalities merged into one, Qeqqata Kommunia, however, harvest continues to be managed separately. Figure 1. North region surveyed by helicopter in 2018, illustrating Maniitsoq and Sisimiut muskox management areas and the 1960’s release location for the 27 translocated muskoxen originating from East Greenland. 17 Government regulated harvesting of the Maniitsoq muskox population began in 1988 with commercial harvest, and since 1993 has supported both recreational and commercial harvesting. Until 1998, annual harvests were typically under 500 muskoxen (Fig. 2). Starting in 2002, annual harvests increased sharply to more than 2,500 muskoxen by 2008 and generally declined thereafter with the 2017-2018 harvests almost half the peak value. Decreasing harvests contrast with steadily increasing hunter effort. For example, prior to the winter harvest season 2000 only hunters from Sisimiut and Maniitsoq participated. The Sisimiut hunters used primarily dogsleds, while the Maniitsoq hunters used ca. ten snowmobiles and in total there were ca. < 25 hunters. Since 2000, the number of motorized vehicles used for transportation to and from the hunting-areas has increased steadily (Hans S. Mølgaard & Nuka M. Lund pers. comm.). Recently, in the 2020 winter-hunt there were 60 hunters with motorized vehicles and in 2021 there were 70 hunters with motorized vehicles (Nuka M. Lund pers. comm.). The decrease in muskoxen harvested despite increased hunter effort suggests fewer muskoxen available. Causes would include decline in muskox population size or that muskoxen are increasingly adept at avoiding hunters. Figure 2. Reported harvest, commercial, recreational, and combined, for the Maniitsoq muskox population from 1962 to 2018. The 1988-1992 period was solely commercial harvest. Trophy harvest not included. Data from Piniarneq records. Further to commercial and recreational hunting (Fig. 2), in the early 2000s trophy hunting began for Maniitsoq muskoxen. By the 2011/2012 season 108 muskoxen were killed as trophies. In 2012/2013 there were 120 trophy muskoxen killed (Cuyler & Raundrup 2014). Trophy hunting is now well established with specific hunting seasons and area concessions (i.e., area allocated for use by one trophy agent). As an industry, trophy hunting continues to grow (Naalakkersuisut 2018). Foreign trophy hunters pay up to Danish kroner 50.000,00 (ca. 7,900.00 US Dollars / 6,700 EUR) per trophy. In 2016, 0 500 1000 1500 2000 2500 3000 1950 1960 1970 1980 1990 2000 2010 2020 2030 Reported muskox harvest YEAR Total harvest Commercial Recreational 18 236 muskoxen were killed as trophies and provided the Greenland government with revenues totally Danish kroner 312.000,00 (ca. 49,200.00 US Dollars / 42,000 EUR). In 2017, trophy hunting rose by 46% with 354 trophy muskoxen killed, which provided Danish kroner 548.000,00 (ca. 86,500.00 US Dollars / 73,700 EUR). Additionally, the Greenland qiviut (muskox wool) industry is founded primarily on the winter harvest of Maniitsoq muskoxen, and owing to increasing demand worldwide, over the past decade the qiviut industry expanded. All the above attests to a substantial economic contribution to local communities from their use of today’s Maniitsoq muskox population, and to a lesser extent also the Sisimiut muskoxen. Given declining annual harvests since 2008, ascertaining population abundance and trend are essential for appropriate management decisions for sustainable harvesting in the future. In the past, infrequent winter ground surveys by snowmobile, over limited areas, provided minimum counts of the number of muskoxen observed (Appendix 1). Within the Maniitsoq management area (Fig. 1, Table 1), minimum counts in the period 20002006 ranged from 4186 to 5092 observed muskoxen. Applying Bayesian analysis to the 2000-2004 minimum counts resulted in a 2004 population estimate of ca. 7,312 (90%CI: 5538-10202) muskoxen in the Maniitsoq management area that was covered by the winter ground surveys (Cuyler & Witting 2004). Ground surveys of the Sisimiut management area were rarely completed. When these occurred, few muskoxen were observed, and ground effort was tiny relative to the management area. This report focuses on the muskox population size estimates for Maniitsoq and Sisimiut attained in conjunction with the 2018 aerial helicopter survey of the KangerlussuaqSisimiut (KS) caribou population in the North region. The name, North region, indicates a relatively northern geographical position within the context of West Greenland and delineates KS caribou distribution. The Government of Greenland’s muskox harvest management areas, Maniitsoq and Sisimiut, are contained within the North region’s boundaries (Fig. 1). Present survey This is the first Conventional Distance Sampling (DS) aerial survey of muskoxen in the Maniitsoq and Sisimiut management areas. This report investigates the DS data set for muskox observations obtained for those areas during GINR’s March 2018 caribou survey, which is described in Cuyler et al. (2021). Initially, we use DS analyses to present the first ever pre-calving population estimates of muskox density and abundance for the Maniitsoq and Sisimiut muskox harvest management areas. Then, 19 we create a Density Surface Model (DSM) for the muskoxen, where density can be spatially represented as a function of additional covariates collected during surveying. The DSM produces alternative estimates of density and abundance for Maniitsoq and Sisimiut muskox harvest management areas. The report provides the first ever population size estimates for muskoxen in the North region and then presents separate estimates for the Maniitsoq and Sisimiut muskox management areas. It also presents information on immediate muskox reaction (movement or lack thereof) to the helicopter fly-by of muskox groups detected, and an approximate calf percentage for the Maniitsoq muskox management area. Note that an earlier analysis for 2018 muskox abundance and density in the Maniitsoq and Sisimiut management areas was run on the same data (Marques 2018), however, the then known areas (km2) were incorrect. Marques’ paper from 2018 is an internal CREEM (Centre for Research into Ecological and Environmental Modelling (St. Andrews, Scotland)) report, which is available upon request. Table 1. Current naming of region, municipality, muskox population, harvest managementand surveyed areas for 2018 in West Greenland that are specific to this report. Region Municipality Greenland Government Muskox harvest management area Muskox Population Surveyed sub-area 2018 North Qeqqata kommunia Sisimiut management area Sisimiut Sisimiut North Qeqqata kommunia Maniitsoq management area Maniitsoq1 Angujaartorfiup 1 Previously referred to as the Kangerlussuaq muskox population, in GINR documents. Methods Study area The North region is within Qeqqata Kommunia in West Greenland. Although Qeqqata Kommunia has ca. 9,400 inhabitants (in 2020), not all live within the boundaries of the North region. The only large settlement within the region is the coastal city of Sisimiut, with ca. 5,600 inhabitants, followed by the ca. 500 residing in the town of Kangerlussuaq. The latter is located on the eastern inland side of region near the Greenland Ice Cap and is also the site of Greenland’s primary international airport. Together, the tiny coastal villages of Itilleq and Sarfannguit contain a further 200-300 people. The North region is seasonally ice-free and covers an area of 23,303 km2, (excluding lakes, rivers, sand, glaciers, and islands). Previous surveys reported a less precise land area of ca. 26,000 km2 (Cuyler et al. 2002, 2005, 2011). Located between 66-68° N Lat, the 20 Arctic Circle (66.5° N Lat) passes through the North region. The northern border is provided by the Nassuttooq Fjord (Nordre Strømfjord). The southern border is formed by two ice caps (i.e., Kangaamiut Sermiat (Sukkertoppen Ice Cap) and Tasersiap Sermia) and the western portion of the Kangerlussuaq Fjord. Elevations reach ca. 1700 on the Kangaamiut Sermiat and ca. 1800 m on the Tasersiap Sermia. The outer half of the Kangerlussuaq Fjord is ice-free year-round and dominated by cliffs of ca. 1000 m. The western border of the region is the permanently ice-free seacoast of the Davis Strait, and eastern border is the Greenland Ice Cap. The North region’s coastal topography is mountainous with peaks whose elevation can be 1000-1800 m and glaciers are common. The Kangerlussuaq Fjord penetrates the region, from SW to NE, ending just before the town of Kangerlussuaq, i.e., close to the Greenland Ice Cap. The Kangerlussuaq Fjord is the primary boundary separating the Maniitsoq and Sisimiut muskox management areas. In the Sisimiut area, moving eastward, the coastal mountains gradually give way to rugged terrain generally ranging 10-900 m elevation. In the Maniitsoq area, immediately north of the two ice caps, Kangaamiut Sermiat and Tasersiap Sermia, the terrain is generally barren highlands with elevations > 1000 m. Continuing northward towards the town of Kangerlussuaq, elevations decrease, and the terrain includes lowland valleys < 400 m elevation and highlands of generally < 1000 m elevation. Common to West Greenland, the North region exhibits a climate gradient on a westeast axis. The western seacoast is wet maritime; however, the climate becomes dry continental as one moves east towards the Greenland Ice Cap. Climate and weather in the west are under the maritime influences of the ice-free Davis Strait and the lowpressure oceanic storm systems that sweep in from the southwest. The climate in the inland of the North region, specifically the Maniitsoq area, is influenced by two ice caps at its southern boundary, the Kangaamiut Sermiat and Tasersiap Sermia. These have elevations of 1,700-1800 m respectively, and act as a barrier to the oceanic storm systems mentioned above, creating a precipitation shadow on the northern side, which in combination with the dominating high pressure over the Greenland Ice Cap creates the inland’s xeric continental climate. Loess/sandstorms are common near the town of Kangerlussuaq (Cuyler et al. 2005). These are caused by katabatic winds, föhn winds descending off the Greenland Ice Cap, which are dry and can have speeds of 30-60 m/s (Putnins 1970, Rasmussen 1989, Tamstorf 2004). Between Kangerlussuaq, and the Greenland Ice Cap, there are two heavily braided rivers, the Akuliarusiarsuup Kuua and Qinnguata Kuussua. The associated valleys exemplify the above conditions, and their Danish names, Sandflugtdalen and Ørkendalen, translate loosely into ‘Blown Sand Valley’ and ‘Desert Valley’, respectively. Further, föhn winds can cause sharp 21 increases in ambient temperature, which in winter or spring can result in extensive snowmelt (Hansen 1999). At Kangerlussuaq winter föhn winds produce large snowfree expanses (Fredskild 1996). In general, the vegetation of the North region may be described as open or alpine tundra. Vegetation is dominated by low arctic species of mainly dwarf shrub heath, which changes to predominantly steppe and grassland when moving east towards the Greenland Ice Cap, where lichen heaths are rare (Tamstorf et al. 2005). Specifically, the Angujaartorfiup sub-area is a grass steppe landscape (Nellemann 1997). Higher elevations are often fell field, abrasion plateaus and bare ground (Tamstorf et al. 2005). Aside from the now firmly established translocated population of muskoxen, native wild mammals present in the North region are caribou (Rangifer tarandus groenlandicus), arctic hare (Lepus arcticus Rhoads) and arctic fox (Vulpes lagopus Linnaeus). Large mammalian predators are absent. Field methods This study was possible owing to the 2018 aerial survey for caribou. The aerial survey occurred 01-15 March 2018 using a helicopter AS350 as the platform for observation. Period and platform were chosen as per criteria for caribou (details in Cuyler et al. 2021). A constant altitude above ground level was maintained while flying low (40 m) and slow (ca. 65 km/hour). Participants included three observers, all with previous survey experience: GINR’s senior scientist Christine Cuyler, GINR’s project coordinator Peter Hegelund, and professional hunter Aslak Jensen (Greenland Association of Professional Hunters (KNAPK)) from Nuuk. Jensen and Hegelund were seated in the rear of the helicopter and observed animals for all distances from the side they were sitting, which alternated each time the helicopter was refueled, which was usually once daily and sometimes twice. Cuyler always sat in front, observed the track line, including distances to either side up to 100 m, and was the data recorder. Verbal contact among the observers permitted the digital audio recording of all observations. Two audio devices (SONY IC recorder, ICD-SX712) were used to record separately the observations specific to the left and right side of the line transect. Audio recording devices were on continual recording for each line transect. At the end of each survey day, audio data was downloaded to computer for storage and back-up. Observations were later paired with Global Positioning System (GPS) coordinates of the helicopter at the time of observation. The audio recording included distance to (see below), and size of, each muskox group 22 observed and name of the observer. Often, behavioral reaction/flight by muskox groups and environmental conditions were recorded. Figure 3. Area covered by the 2018 caribou survey of the North region (23,303 km2). Three different colours illustrate the three sub-areas, designated as Sisimiut (blue), Sisimiut South (orange) and Angujaartorfiup (purple). The term ’Byer bygder’ refers to city/town/village identified in Fig. 1. Survey design The surveyed North region area, 23,303 km2, was divided into three sub-areas (strata), arbitrarily named Sisimiut (12,658 km2), Sisimiut-South (3,512 km2) and Angujaartorfiup (7,133 km2) (Fig. 3). The sampling design for the 2018 survey considered 19 systematic parallel line transects of variable length separated by 15 km and placed over the three sub-areas (Fig. 4). Those transects provide the maximum area coverage possible given the financial resources available. An initial line transect was computer generated at random, and the others followed 15 km apart. Aligning line transects perpendicular to known gradients within the surveyed area can maximize precision of the resulting estimate by lowering the encounter rate variance (Buckland et al. 2001). Thus, the transect axis direction was chosen as perpendicular to previously known animal distribution gradients in March. Lines 1 to 13 followed a west-east axis, which also reflects the climate gradient from wet maritime to dry continental. Line transects 14 to 19 followed a north-south axis, reflecting animal, climate, and topological gradients that exist between the combined ice caps (Kangaamiut Sermiat and Tasersiap Sermia) and the town of Kangerlussuaq. 23 Figure 4. The 19 line transects used in the 2018 survey of the North region, employing the same three colours as applied to the three sub-areas in above figure 3: Sisimiut (blue), Sisimiut South (orange) and Angujaartorfiup (purple). Distance collected was the perpendicular distance from the helicopter’s flown track line to a muskox group (object-of-interest), i.e., one or more animals. Tightly cohesive behavior identified groups of multiple individuals. Distance recorded was the observer’s instantaneous and subjective distance to the approximate center of the muskox group from the track line before any movement by the group occurred. Exact distance measurements were not possible primarily because of practical considerations for this survey (details in Cuyler et al. 2021). Instead, distance measurement used the following distance bins: 0, 50-, 100-, 200-, 300-, 400-, 500-, 750and 1500-meters perpendicular to the line transect. These values correspond to the upper limit for a specific bin that the muskox observation was included in. For analysis, these were recoded to the mid distance for a specific bin. Note, binning accuracy relies heavily on observer ability to correctly estimate distance to the observed animals. Thus, before starting the survey the helicopter hovered at the 40m altitude used during line transects, while each observer used a Leica laser range finder 1600 to gauge distances across the terrain. Then they marked their window with masking tape delineating the approximate distances for each bin. When possible while flying line transects, the laser range finders were used to double-check reported bin distances to detected groups. The recorded distances are used to estimate a detection function, then estimate the detection probability and finally to estimate the density of the muskoxen within the 24 surveyed area (Buckland et al. 2001). The detection function, 𝑔(𝑦), describes the probability of detecting an object-of-interest given that is at a distance 𝑦, from the centerline (0-line), thus being a non-increasing function of 𝑦 (Buckland et al. 2015). For line transects, 𝑦 is the perpendicular distance from the 0-line to the detected object. Within DS methods, the probability of detection is explained recurring to these observed distances (Buckland et al. 2001). Distance Sampling The muskox group was the selected sample unit for the DS analysis of the 2018 survey. Neither the individual muskoxen within a group, nor individual line transects were considered as the sample unit. The recorded distances to the observed muskox groups were used to estimate a detection function. With this, both the muskox detection probability and density within the surveyed area could be estimated (Buckland et al. 2001). The detection function describes the probability of detecting an object-of-interest (muskox group) that is at a distance 𝑦, from the centerline (track line), thus being a non-increasing function of 𝑦 (Buckland et al. 2015). For line transects, 𝑦 is the perpendicular distance from the track line to the detected object. Within DS methods, the probability of detection is explained recurring to these observed distances (Buckland et al. 2001). Prior to DS analysis, the raw data was first processed for inconsistencies, e.g., species or observer names written slightly differently had to be standardized before analyses to avoid being assigned a different category, and one observation lacked distance (replaced by the average observed distance). Then extensive exploratory data analysis was completed, including evaluation of observed distances, before proceeding to determining the detection function through model fitting and selection (Buckland et al. 2001; Marques et al. 2011; Thomas et al. 2010). To determine the detection function, several models were considered (Thomas et al. 2010). The model presenting the lowest AIC value was chosen. The subsequent analysis was based on Marques (2018). Details regarding DS theory, methods and analysis are available in Buckland et al. (1993, 2001, 2015), and a briefer summary provided in Appendix 2. For analysis, we used R Statistical Software (https://www.r-project.org/). Generalized Additive Model (GAM) Generalized Additive Models (GAM) are an extension to Linear Models (LM) and Generalized Linear Models (GLM) where non-linear responses with smoothing functions can be fitted to the data. Details regarding Generalized Additive Models, 25 methods and analysis are available in Wood (2017), and a briefer summary is provided in Appendix 3. Density Surface Model (DSM) Conventional DS methods provide average estimates of abundance over a region but no information about the distribution of the objects of interest within the survey region. An efficient option is to build a spatial model that incorporates spatially referenced environmental covariates. Density surface modelling uses the GAM framework (Wood, 2017) to build models of abundance/density as a function of environmental covariates, typically as part of a two-stage method (Fig. 5). In the first stage, the detectability via DS is modelled and in the second stage the counts, corrected for detectability, are modelled over space. Details regarding DSM , methods and analysis are available in Katsanevakis (2007) and Miller et al. (2013). Briefer summary is in Appendix 4. Figure 5. Flow diagram of the modelling process during the analysis, from the Distance Sampling analysis through to the Density Surface Modelling. 32 Regarding groups of muskoxen, the number of detections per unit transect length is the encounter rate. Considering the entire North region altogether, non-truncated data had a mean encounter rate of 0.016 muskox groups per kilometer. Specifically, the Sisimiut sub-area had a mean encounter rate of 0.003 muskox groups per kilometer, while the Angujaartorfiup sub-area had a mean of 0.026. No muskoxen were encountered in the Sisimiut South sub-area. Within the sub-area Sisimiut, encounter rates were zero for three line transects and sparse for the remaining five. For the Angujaartorfiup sub-area, the encounter rates were high in two possible hot spots, the east and west sides of the sub-area, and were least in the middle (Fig. 8). Histograms examining observer effects (Fig. 9) were somewhat dissimilar and therefore this covariate will be investigated as to whether it influenced results for detectability. Other potential covariates, like sun glare or snow covering, were available but there were too many missing observations and/or inconsistency when referring to the categories for these to be used in the analysis. The preliminary analysis provides the expectation of less precision in further analyses of detections within the Sisimiut sub-area, because of data variability. Precision is expected to be greater for Angujaartorfiup, as all line transects contained muskoxen. Nevertheless, the information agreed well with anticipated a priori, e.g., Angujaartorfiup sub-area would have more muskoxen than the other two sub-areas. Figure 8. Muskox survey 2018: exploratory analysis of non-truncated data for muskox encounter rate (groups per km) per line transect using and illustrating sub-area: Sisimiut (blue, line transects 1-8), Sisimiut South (orange, line transects 9-13) and Angujaartorfiup (purple, line transects 14-19). Line transect number 33 Distance Sampling analysis Before conducting any modelling, an analysis of the observed distances was made to evaluate whether any major assumption violation occurred or other data-related issue, Figure 9. Observer effect: histograms illustrating detected distances for the two observers (a covariate with two levels). Density, y-axis, refers to the density of observations. Figure 10. Histogram of the different binning options for the muskox distance data. Left: the original bins as collected on the survey. Right: an alternative binning to reduce the effect of heaping. The area of the rectangles is proportional to the number of points within each bin. 34 as stated in previous sections. These analyses are from Marques (2018). The histogram of observed distances with no defined truncation distance (non-truncated data) is like typical Distance Sampling data, perhaps showing some over-dispersion, with notequally-spaced bins (Fig. 10). Given the histogram of binned distances, a strip halfwidth of 𝑤 = 0.75 km was selected (i.e., all observations at distances beyond 750 meters were discarded). This truncation reduced the sample size from 256 to 210 muskox groups for the Distance Sampling analysis. Data truncation is a common procedure because otherwise extra adjustment terms may be needed to fit the long tail of the detection function. Further, little information is lost by truncation, since data observations located more than 0.75 km from each side of the line make a minimal contribution to the abundance estimate. Truncated data With the original binning Option 1 (Fig. 10), there seem to be less than expected observations on the 0.05-0.10 km and 0.30-0.40 km intervals, when compared to the 0.10-0.30 km and the 0.40-0.50 km bins. This suggests heaping, a phenomenon that occurs when observers tend to record some preferred values over others (Buckland et al. 2001) e.g., at round distances that are easily chosen in the absence of a rigorous distance measuring method, as often occurred in this study when rugged terrain forced observers to assign distance bins subjectively. The alternative binning Option 2, with less bins, reduces the influence of potential measurement errors in the observed distances. This alternative binning option includes bin cut points of 0, 0.10, 0.30, 0.50 and 0.75 km (Fig. 10). Both binning options were considered in model fitting, albeit only Option 2 minimizes the effect of measurement error induced by heaping. Since binning Option 1 was not suitable for grouping, only the analyses whose fitted models consider the second binning option are illustrated below. The advantage for choosing the second binning option is that it results in more reliable detection functions. However, owing to fewer degrees of freedom, the small number of bins affects the 𝜒2 Goodness-of-Fit tests following model fitting. The regression analysis suggests that distance is not a statistically significant variable explaining group size (Table 5). Table 5. Summary of the coefficient characteristics of the GLM between observed truncated distance and the muskox group size while considering a Poisson distribution. Parameter Estimate Standard Error z-value p-value Intercept 1.752 0.050 35.18 0.00000 Distance -0.091 0.155 -0.58 0.55926 Note: AIC = 1564.2, Null Deviance = 879.5, Residual Deviance = 879.2. 35 Detection function models fitted with the first binning option did show poor fitting, including the best fit within this group, since these presented several adjustment terms, due to the heaping phenomena. Below, the detection functions are fitted to truncated data and considering the second binning option (Fig. 10), which improves the models. Truncated data – Distance Sampling Models For these models, every combination of key function and adjustment terms was tested. The only additional covariates assessed were observer and group size, considering 𝑤 = 0.75 km. A summary of the information from each model fitted to the data (Table 6) provides a simple overview of several models, and includes the respective key functions, adjustment terms, model formula , 𝜒 2 Goodness-of-Fit test p -value, estimates of the detection probability, respective standard error ( se (𝑃 a) ), and Δ𝐴𝐼 𝐶 comparison between each model and the model with the lowest AIC . The best model fitted to the data possesses the lowest change in AIC value (Δ𝐴𝐼𝐶 = 0). For the 2018 muskox survey data, the Half-Normal key function was selected because it was the most flexible key (AIC = 550.65). The second-best model was Uniform with cosine adjustment terms of order 1 (AIC = 550.692, i.e., Δ𝐴𝐼𝐶 = 0.042). The best fitted detection function parameters superimposed with the observed distances’ histogram, indicates that neither group size nor observer were relevant covariates in detectability for the truncated data. (Table 6). The estimated averaged probability of detection for the North region was 𝑃 a = 0.557 (se = 0.033, Table 6). Figure 11. The detected distances truncated data with the estimated detection function overlaid, considering the binning option 2 that reduces the effect of heaping. 36 Each group size has its separate detection function, corresponding to different estimates for the probability of detection (Fig. 12). For non-truncated dataset, a group size of 2 muskoxen presents an estimated probability of detection of 0.273, a group size of 10 has an estimate of 0.45, a group size of 20 has an estimate of 0.74, while a group size of 30 has an estimate of 0.96 (Fig. 12). With increasing group size, the probability of detection also increases. For the truncated dataset, the probability of detection (0.557) was the same regardless of muskox group size. The Kolmogorov-Smirnov and Cramérvon Mises tests (Appendix 2) cannot be applied since the distances were represented as a discrete variable. Figure 12. Estimated probabilities of detection, given non-truncated and truncated dataset, for each observed group size (muskoxen) obtained with the fitted model. Illustrates group size not important in truncated dataset. Figure 13. Muskox encounter rate (groups per km) per line transect, survey 2018: exploratory analysis of truncated data, illustrating sub-area: Sisimiut (blue), Sisimiut South (orange) and Angujaartorfiup (purple). 37 Encounter rate estimates suggest the Angujaartorfiup sub-area has the most muskoxen, since its estimate is larger than the other sub-areas (Table 7, Fig. 13). Visualization of the detected muskox distribution shows this was somewhat continuous only for line transects in the Angujaartorfiup sub-area, with two ‘hot’ spots, i.e., east, and western sides of the sub-area (Fig. 14). Muskox distribution in Sisimiut was sporadic and rare, and in Sisimiut South nonexistent. The DS design-based estimates for muskox abundance and density also reveal that Angujaartorfiup is the sub-area with the most muskoxen (Tables 8, 9, Fig. 15). Figure 14. ‘Heat’ map illustrating relative distribution of muskoxen observations along the line transects in the North region, West Greenland. For interpretation of underlying map see Fig. 1. In March 2018, the DS design-based population estimate for the combined Angujaartorfiup and Sisimiut sub-areas was approximately 21,746 muskoxen (95% CI: 11,061 – 42,751), with a CV of 28.5% (Table 8). The DS density estimate was 1.1 muskoxen per km2, with 95% CI: 0.56 – 2.16 (Table 9, Fig. 15). Separately, the Maniitsoq muskox harvest management area (Angujaartorfiup subarea) had a DS design-based population estimate of 18,906 muskoxen (95% CI: 8,726 – 40,960), with a CV of 31.5%. Similarly, Sisimiut muskox harvest management area had a population estimate of 2,840 muskoxen (95% CI: 662 – 12,178), with a CV of 56.8%. 38 Table 6. Model comparison across the three Conventional Distance Sampling models and models considering muskox group size and observer as covariates. Key function Formula 𝒙𝟐 p-value 𝑷 a se (𝑷 a) ∆ AIC Half-normal (same result for Cosine/Simple Polynomial/Hermite) 1 0.932 0.557 0.033 0.000 Uniform with cosine adjustment terms of order 1 NA 0.913 0.565 0.025 0.042 Uniform with Hermite polynomial adjustment term of order 4 NA 0.448 0.604 0.024 1.509 Half-normal Group size 0.717 0.557 0.033 1.673 Half-normal Observer 0.709 0.557 0.033 1.953 Hazard-rate (same result for Cosine/Simple Polynomial/Hermite) 1 0.616 0.584 0.054 2.111 Uniform with simple polynomial adjustment terms of order 2,4,6 NA 0.498 0.577 0.047 2.325 Hazard-rate Group size NA 0.596 0.052 3.575 Hazard-rate Observer NA 0.589 0.053 4.004 Note: under Formula explanatory variables are as follows: Group size = group size as variable, 1 = for Uniform key, Observer = observer as variable, NA = no explanatory variables/covariates. There were not enough degrees of freedom for the 𝜒2 Goodness-of-Fit test, thus the ‘NA’ values. For each key function, all three series expansions (Cosine, Simple Polynomial, Hermite (Appendix 2, Table 15) were applied. Table 7. Encounter rate (ER) estimates per sub-area for muskox groups considering truncated data, three sub-areas (strata), four bins, and a Half-Normal detection function fitted without covariates. Surveyed Sub-area* Muskox harvest Management area Encounter rate Standard Error (se) Coefficient of Variance (cv) Sisimiut Sisimiut 0.020 0.007 0.360 Angujaartorfiup Maniitsoq 0.414 0.092 0.222 TOTAL Sisimiut + Maniitsoq 0.173 0.067 0.387 *Third sub-area, Sisimiut-South, does not appear in the table because there were no observations of muskoxen. Table 8. March 2018, the Distance Sampling design-based muskox abundance estimates per stratum (sub-area) in the North region, considering truncated data, three sub-areas (strata), four bins and a Half-normal detection function with no covariates. Sub-area* Muskox harvest Management area Population Size Estimate SE CV 95% Confidence Interval Lower Upper Sisimiut Sisimiut 2,840 1,613 0.568 662 12,178 Angujaartorfiup Maniitsoq 18,906 5,948 0.315 8,726 40,960 TOTAL Sisimiut + Maniitsoq 21,746 6,194 0.285 11,061 42,751 Note: SE = Standard Error, CV = Coefficient of Variance. *Third sub-area, Sisimiut-South, does not appear in the table because there were no observations of muskoxen. Table 9. March 2018, the Distance Sampling design-based muskox density estimates per stratum (sub-area) in the North region, considering truncated data, three sub-areas (strata), four bins and a Half-normal detection function with no covariates. Sub-area* Muskox harvest Management area Density Estimate SE CV 95% Confidence Interval Lower Upper Sisimiut Sisimiut 0.224 0.127 0.568 0.052 0.962 Angujaartorfiup Maniitsoq 2.650 0.834 0.315 1.223 5.742 TOTAL Sisimiut + Maniitsoq 1.099 0.313 0.285 0.559 2.160 Note: SE = Standard Error, CV = Coefficient of Variance. *Third sub-area, Sisimiut-South, does not appear in the table because there were no observations of muskoxen. 39 Figure 15. Distance Sampling design-based muskox density estimates with corresponding confidence intervals for the two sub-areas with muskoxen, Sisimiut and Angujaartorfiup, and finally for the total North region. GAM & DSM After splitting the transects into segments, their centroids were determined and intersected with the shapefiles concerning the spatial variables. These variables were GPS coordinates, vegetation, elevation, slope, and the aspect (geographical compass direction, which relates to sun exposure and temperature of a particular location) (Appendix 6). The vegetation covariate was excluded from the analysis since the current information available lacks the necessary resolution. Furthermore, a prediction data set was generated for the whole study region by converting the region into small 1.5 km sided square cells (i.e., 2.25 km2, note 𝑤 = 0.75 km). Once the distance model was fitted, 𝑃 a could be determined for each group size and thus 𝑛𝑖, the number of muskoxen detected in group i, can be corrected as 𝑛i = 𝑛𝑖 𝑃 a. These predicted values for counts were then modelled using a GAM fitted to the spatial covariates along longitude (lon) and latitude (lat) considered jointly. Within the model, the previously estimate probability of detection was defined as the offset term along with the cell area of 2.25 km2. The distribution considered in model fitting were Tweedie and Negative Binomial. These provide a flexible alternative to the quasiPoisson distribution, which does not capture the response overdispersion. Table 10. GAM model summary table relative to smooth terms for covariates and considering Tweedie distribution. Effective Degrees of Freedom Chi-square p-value s(lon, lat) 19.235 9.021 0.0000 s(elevation) 4.024 4.238 0.0009 Note: 𝑅𝑎𝑑𝑗 2 = 0.378, Deviance explained = 63.6%. 40 Table 11. GAM model summary table relative to smooth terms for the covariates and considering a Negative Binomial distribution. Effective Degrees of Freedom Chi-square p-value s(lon, lat) 19.667 159.429 0.0000 s(aspect) 1.000 4.215 0.0401 s(elevation) 3.692 15.404 0.0123 Note: 𝑅𝑎𝑑𝑗 2 = 0.091, Deviance explained = 67.6%. Fitting GAM/DSM The two distribution models, Tweedie, and Negative Binomial (Tables 10, 11) were fitted and the best model within each was selected for further work. The Tweedie distribution used longitude, latitude, and elevation. The Negative Binomial distribution (NegBin) used longitude, latitude, aspect, and elevation. The third best NegBin model was similar to the first (i.e., Table 11) but lacked elevation and the ∆AIC was low, 0.8799. The second best, which included all covariates, had a ∆AIC difference of 0.61. Thus, these models were almost identical. Figure 16. Fitted DSM plot corresponding to the fitted smooth of the elevation covariate (Y-axis is the fitted smoothing parameter). Green shading illustrates the standard error of the estimates, while the respective observations are represented on the horizontal axis. Despite the similarity among results, the model considering the Tweedie distribution was chosen owing to its lower AIC value (AICTweedie = 1350.1 and AICNegBinom = 1461.5). A larger Effective Degrees of Freedom (EDF) resulted for the bivariate smooth associated with the paired longitude and latitude relative to the other environmental covariates, since more basis functions are required to fit a surface than a line. Each of these environmental covariates does not appear to be linearly related with the response, since the respective EDF is larger than 1. Regarding the smooth functions, the relationship between the response variable and each explanatory variable appears to be Elevation (m) 41 non-linear for the elevations below 1000m (Fig. 16). Specifically, elevations from 300m to 700m seem to be preferred by the muskoxen in the early March period of the survey and density was positively correlated with elevation until around ca. 400 m when it began to fall slightly with elevation (Fig. 16). At high elevation, e.g., above 900m, variability was greater, and density fell. Muskoxen appear to prefer south-facing slopes (90°-270°), since these involved most observations and there was the least variability in the estimates. Furthermore, the results suggested that among south facing slopes, southwest aspects were most preferred by muskoxen. Spatial prediction The number of muskoxen within each cell was then predicted using the GAM model (𝑛i) and spatially represented. The heat map presenting the prediction of muskox abundance in the North region indicates that muskoxen are scarce in most of the North region, excepting the Angujaartorfiup sub-area, which clearly illustrated the muskoxen distributed into two clusters (hot spots) at the time of the survey (Fig. 17). The GAM/DSM model-based approach produced a predicted mean muskox density value for the entire North region of 3.69 muskoxen/km2 (95% CI: 2,87 – 4,74). The predicted range was from 0 to 125 muskoxen/km2 (Fig. 17). The maximum value occurred rarely and only on a remote section of line transect 14, which was on the far west side of the Angujaartorfiup sub-area. This was one of the two hot spots of clumped distribution where most of the darker cells were in the range of 25-30 muskoxen/km2. The clumping coincided with unusually high elevation. The GAM/DSM model-based approach produced a population estimate of 23,256 (95% CI: 18,102 – 29,877) muskoxen for the entire North region. Additionally, a variability map was produced with the CV for each estimate in the survey region (Fig. 18), the CV estimate was 11.36%. Table 12. Comparison of Distance Sampling design-based and GAM/DSM analyses, for entire North region. Analysis Population Size CV for probability of detection Estimate SE CV 95% Confidence Interval Lower Upper Distance Sampling 21,746 6,194 0.285 11,061 42,751 5.98% GAM / DSM 23,256 0.114 18,102 29,877 11.36% Note: SE = Standard Error, CV = Coefficient of Variance. To compare the model-based GAM/DSM to the design-based DS approach we used the CV for probability of detection of the muskoxen (i.e., not for abundance or density). The GAM/DSM had a probability of detection CV of 11.36% (like the estimate’s CV), while the DS had a probability of detection CV of 5.98% (Table 12). For the overall analyses, i.e., DS and GAM/DSM, the probability of detection CV was 12.84%. 48 season and during the March trophy season Maniitsoq muskoxen inhabited unusually high elevations, >700 m, we suggest that current anthropogenic disturbance associated with hunting may play a major role in causing muskox use of high elevations. Consequences of winter foraging in sub-optimal habitat Regardless of the cause(s), in 2018 most Maniitsoq muskoxen in the Angujaartorfiup sub-area foraged at elevations of 700-800 m in winter. This is suboptimal habitat. Given the latitude (ca. 67°N), winter forage available at 700-800 m elevation is inferior and less plentiful than that found in lowland valleys (Körner 2007), and specifically for the Angujaartorfiup sub-area high elevations are often fell field, abrasion plateaus and bare ground (Tamstorf et al. 2005). Foraging in poor habitat at high elevation could negatively affect winter body condition and ultimately survival of Maniitsoq muskoxen (see Appendix 9 for details). Further, it is common knowledge that among large herbivores most fetal growth occurs in late gestation and that this growth substantially increases maternal energy requirements. Muskox calves are usually born mid-April through June, with birthing dates as early as 05 April and as late as 19 June (Lent 1988) as there is no breeding synchrony as seen in caribou (Lent 1966). Depending on the date for parturition, late gestation for muskoxen ranges from February through May. Thus, Maniitsoq cows foraging at high elevation poor pasture could become nutritionally deficient in late gestation and birth weaker calves of low weight (Olesen et al. 1994). Calf percentage Forage quality varies widely seasonally, with winter associated with lowest forage quality (Dermanet et al. 2015). The Angujaartorfiup sub-area straddles the Arctic Circle (66°33’ N). At these latitudes, forage decreases rapidly with rising elevation (Tamstorf et al. 2005, Körner 2007). Muskoxen preferring to forage in lowlands year-round is therefore no surprise. We suspect that observed winter use of vegetation poor high elevations, 700-800 m, by muskoxen in Angujaartorfiup sub-area, could compromise their body condition, with negative consequences for body condition and finally winter survival and calf production. Since muskoxen are capital breeders (Desforges et al. 2019), cows rely on body reserves for parturition and early lactation. Thus, muskox cows entering parturition with low body condition may birth calves but lack sufficient reserves for milk production, which may reduce calf survival. Thus, the winter hunting season coinciding with late gestation is not expected to be compatible with a high percentage of calves. The general decline of calf percentage in the 2000-2010 period and additionally the March 2018 calf (age <1-year) percentage of ca. 18.4%, appears to support this. The March 2018 value is at the low end of the range 17-24% reported for expanding Alaskan and Canadian muskox populations (Jingfors & Klein 1982, Gunn et al. 1984). Further, it is well below the excellent 32% that characterized the Maniitsoq 49 muskox population in the 1970’s and into the late 1980’s (Olesen 1993, Cuyler & Witting 2004, Appendix 9). In addition to the possible role of the winter hunting season in reducing calf percentage, too many animals in an area relative to the available quality and quantity of forage can bring forth density-dependent factors that may reduce calf production in large herbivores (Kie & White 1985, Ouellet et al. 1997). Density-dependent factors are suspected of playing a role in the steady decline of Maniitsoq muskox calf percentage in the 2000–2004 period (Cuyler & Witting 2004). At that time, observed densities on preferred lowland winter forage averaged from 4 to 16 muskoxen/km2 (Cuyler & Witting 2004), with maximums of 21-29 muskoxen/km2 in valleys near the town of Kangerlussuaq (Cuyler et al. 2001). The minimum counts of 2000, 2001, 2002 and 2004 observed similar numbers of muskoxen (Cuyler & Witting 2004). Initial calf percentage was 26%, and dropped with each count, 25%, 21% and 18%, respectively. The latter two were also below those predicted by Bayesian modelling (Cuyler & Witting 2004). Meanwhile, the rough density of muskoxen in lowland elevations appeared to rise, being ca. 2.5/km2 in 2000 and 3.1/km2 in 2004 (Cuyler & Witting 2004). The high densities were observed prior to the winter hunting season and made densitydependent factors a plausible cause behind the declining calf percentages in 2000-2004. The observed decreasing winter cow rump fat depth and pregnancy rates for the same period support this, since not due to temporary factors, e.g., weather events, since similar weather conditions applied across years (Cuyler & Witting 2004). Further, winter harvesting may also have played a role in the falling calf percentages. Already from winter 2000, once snowmobile use commenced for winter hunting then most muskoxen left the valleys and moved into higher elevations (Cuyler unpublished). Given high elevation pasture is poor relative to lowlands and previous muskox densities in lowlands, this movement of muskoxen to higher elevations may have exacerbated the density-dependent factors suspected of already operating in the lowlands. In winter 2018 most of the estimated population of 18,906 Maniitsoq muskoxen were observed at elevations above 700 m, which provide sub-optimal habitat for foraging. Since calf percentage was only approximately 18% calves, population growth is possible (albeit likely slow) but not certain. Under current harvest regimes that calf percentage may be insufficient for population size stability. With almost 19,000 Maniitsoq muskoxen on poor high elevation pasture in winter, density-dependent effects may be expected. Calf percentage might improve if Maniitsoq muskox cows were able to forage undisturbed in lowlands during late gestation. 50 Muskox population size & density Survey coverage of the North region (23,303 km2) was 10.6%, which is usually sufficient to facilitate reliable estimates. Both DS designand GAM/DSM model-based approaches were applied to the dataset to obtain comparable estimates of muskox abundance and density. DS design-based estimates are based on the selected transects, which may or may not adequately represent every feature within the study region, as some features may be over-represented, while others under-represented. Meanwhile, GAM modelling considers the environmental covariates of the whole region, allowing the spatial representation of the estimates obtained and a visualization of patterns in abundance. Most muskoxen were observed in the Angujaartorfiup sub-area (Maniitsoq muskox harvest management area), few in the Sisimiut sub-area (Sisimiut muskox harvest management area), and none in the Sisimiut South sub-area. Considering the combined Angujaartorfiup and Sisimiut sub-areas, the DS analysis estimated muskox abundance at ca. 21,746 muskoxen (Table 12). That estimate, however, had a wide 95% CI (11,061 – 42,751), which indicates the range of possible abundance is broad. Further the CV was large, 28.5%, illustrating the magnitude of imprecision on the estimate. The GAM/DSM analyses produced a slightly larger estimate of 23,256 muskoxen, which had a tighter 95% CI (18,102 – 29,877) and a lower CV of 11.4%. The two abundance estimates may be considered similar because they differ by only 7-8% and 95%CIs overlap. The DS population estimate was not precise as evidenced by the large CV. The main factor contributing to imprecision was the high variability within the dataset. For example, most transects flown had zero muskox observations while a few transects had many. Further, the number of muskox group observations was low (n = 210, truncated dataset) combined with a large range in the number of muskoxen individuals within groups. Specific to the Angujaartorfiup sub-area (Maniitsoq muskox harvest management area) was the unusual, clumped distribution of the muskoxen into two hotspots and few muskoxen elsewhere. Specific to the Sisimiut sub-area (Sisimiut muskox harvest management area) was the enormous area surveyed and the dominance of zero observations. Considering the sub-areas surveyed, the DS estimate for the Angujaartorfiup sub-area was ca. 18,906 Maniitsoq muskoxen (95% CI: 8,726–40,960; CV = 31.5%) (Table 8). Similarly, the Sisimiut sub-area was ca. 2,840 Sisimiut muskoxen (95% CI: 662–12,178; 51 CV = 56.8%). For each sub-area, like with the overall estimate, the DS provided imprecise estimates with a broad range of possible population size. At first glance, the GAM/DSM estimate for the entire North region appears best, however, to permit comparison of the DS and GAM/DSM approaches we used the CV for probability of detection of muskoxen within the entire region. The result was that DS was better than the GAM/DSM, with CVs of 5.98% and 11.36%, respectively (Table 12). This indicates that despite uncertainty the design-based DS approach provided a more reliable population estimate than the model-based GAM/DSM for this dataset. Regarding harvest management decisions, we suggest using the design-based DS estimates for abundance, albeit fully aware of large imprecision. The DS also provides estimates specific for each population, i.e., Maniitsoq and Sisimiut, which can be helpful for applying population specific harvest management. Given the uncertainty in the estimates, erring towards caution is suggested regarding abundance, i.e., consider values in the lower range of the 95% CIs. Regarding muskox density, the design-based DS density estimates are assumed more immediately usable than those from the GAM/DSM. The design-based DS estimate for the combined Angujaartorfiup and Sisimiut sub-areas was 1.1/km2 (95% CI: 0.56 – 2.16). Specifically, the Angujaartorfiup sub-area density was ca. 2.65 muskoxen/km2 and the Sisimiut sub-area ca. 0.22 muskoxen/km2. These muskox densities intuitively match observer experience. However, density estimates from the model-based GAM/DSM approach did not. The GAM/DSM density estimate for the entire North region was 3.69 muskoxen/km2 (95% CI: 2,87 – 4,74), which although not realistic is reasonable given the high concentration of muskoxen at two hot spots. Similarly improbable, estimated GAM/DSM densities had a maximum of 125 muskoxen/km2 (Fig. 17) for rare high elevation locations in the Angujaartorfiup sub-area. We do not suggest a literal interpretation. Instead, the maximum GAM/DSM value reflects the extreme degree muskox distribution was clumped in the Angujaartorfiup sub-area during the survey period. It is inherently harder to estimate abundance and density when clumping is severe. This dataset had enormous variability, which in addition to clumping also included highly variable group size. Additionally, it is important to consider the variability in the density data associated with the covariates (lon, lat, and elevation). The hotspot with the highest densities was associated with elevations that typically ranged from 300 m to 500 m. However, the second hotspot was associated with lower elevations, which made the density estimates lower, albeit the densities were still large (e.g., 70 muskox/km2). Further, the Sisimiut sub-area has extensive terrain with elevations from 300 m to 500 m (Appendix 6 - Fig. 32), but no hotspot 52 occurred. Instead, muskox detections were sparse in Sisimiut, and the absence of an estimated hot spot is based on the muskox observations. Finally, the largest estimated GAM/DSM densities in the hot spots are associated with the highest uncertainty (Figs. 17, 18). The above makes estimating abundance and density difficult. Nevertheless, the GAM/DSM density pattern (Fig. 17) clearly illustrates the clumping of Maniitsoq muskoxen into areas remote from human influence. As there is less precision on the GAM/DSM density and the CIs from the two approaches do not overlap, we suggest giving the GAM/DSM density less ‘weight’, although values in the lower CI range may be relevant at specific locations. We expect the design-based DS densities are closer to actual densities and recommend their use in management decisions. In 2018 the Maniitsoq muskox population in the Angujaartorfiup sub-area had a density, 2.65/km2. Elsewhere in the Arctic (Canada’s Northwest Territories and Nunavut, and Northeast Greenland), muskox densities are typically much lower, having a range from 0.02 to 1.8/km2 and averaging ca. 0.86 muskoxen/km2 (Cuyler et al. 2001). In the period 2000-2004 Maniitsoq muskoxen density in lowland elevations <400 m was ca. 4-5/km2 (Cuyler et al. 2001, Cuyler and Witting 2004). These values, however, applied only to the northern portion of the Angujaartorfiup sub-area. Regardless, both present and past Maniitsoq muskox densities are relatively high, which attests to the habitat suitability of the Angujaartorfiup sub-area for muskoxen. This is further illuminated by the extreme densities encountered in the 2000-2004 period for muskoxen foraging in lowland elevations <400 m. Examples from core winter grazing areas included, 21.6 muskox/km2 in the valley, Ørkendalen (n=589 muskoxen in area of 20.25 km2), and 29.1 muskox/km2 in the lowlands and river valleys surrounding the Ammalortoq Lake (n = 876 muskoxen in area of 40.5 km2) (Cuyler & Witting 2004). The falling calf percentages of the 2000-2004 period, suggest that, at those densities, density-dependent factors may have begun negatively affecting the population. A study to examine the current herbivore carrying-capacity of the Angujaartorfiup subarea’s lowland grass/sedge Kobresia steppe is warranted, as is a similar study for elevations of 700-800 m. Specifically, examining the effect of foraging and trampling on vegetation at current muskox density, while considering that föhn winds create expanses of snow-free winter pastures where the vegetation is exposed to harsh conditions including temperatures well below zero and powerful winds. 53 Muskox population trend From the 2018 survey, the DS and GAM/DSM analyses provided population size estimates for muskoxen in the North region. These are the first population size estimates based on aerial survey data. Alone, the 2018 estimates cannot indicate current population trend, because this would require at least two additional aerial survey estimate points using similar methods. Additional future surveys will be needed. In the past, Maniitsoq muskoxen (Angujaartorfiup sub-area), received ground-based minimum counts (Appendix 1: Table 13, Appendix 9). However, these cannot be used for comparison to the 2018 aerial estimate, primarily because those counts report only the number of animals observed and do not estimate population size. Regardless, the 2018 estimates of muskox population size clearly illustrate an immense population growth since 27 individuals were translocated to the area over 50 years ago. For the Maniitsoq muskoxen, population growth appears to have been steepest in the late 1980’s and early ‘90’s (Appendix 9). In the 2000-2020 period, ground-based winter minimum counts observed that total muskox number and number of groups became fewer, while simultaneously group size, calf percentage and calf recruitment decreased (Appendix 9). Combined, these indicate a declining Maniitsoq population since the mid-2000s, as does the plot of observed muskoxen from minimum counts for the Angujaartorfiup sub-area (Fig 22). Figure 22. Index of Maniitsoq muskox abundance in the Angujaartorfiup sub-area since translocation to the region in 1963-65. Data are minimum counts of only those individuals observed and cannot be interpreted as estimates of population size. All but 2018 are ground counts. All are of varying timing, effort, and area coverage. Beginning 2000, all are pre-calving counts from either winter or late winter. Those from 2000 to 2010 were pre-harvest, while those 2014-2020 were post-harvest. In the latter muskoxen were avoiding preferred lowlands. The actual number of individual muskoxen observed during the March 2018 helicopter survey () has been included as a type of minimum count, which effort and coverage were the most comprehensive relative to any other minimum count presented. (From Cuyler 2020). 0 1000 2000 3000 4000 5000 6000 1950 1960 1970 1980 1990 2000 2010 2020 2030 Number muskoxen observed Year 54 Reported harvests also suggest declining muskox numbers. Harvest increased after 2002, the peak was in 2008 and generally has declined ever since (Fig. 2, Appendix 8: Table 17). The declining harvests following 2008 (Fig. 2) are reflected in the declining index of muskox abundance starting from about the mid-2000s (Fig. 22). Bayesian analyses of the 2000-2004 minimum counts in the Angujaartorfiup sub-area concluded that harvesting more than ca. 1600 muskoxen annually would likely reduce Maniitsoq muskox population size and harvests that large would not be sustainable over the long term (Cuyler & Witting 2004). Nevertheless, reported annual harvest from 2004 to 2016 always exceeded 1600 muskoxen (Appendix 8: Table 17). Seven of those years exceeded 2000 muskoxen annually. Then, the harvests in 2017 and 2018 were the lowest in over a decade, falling below 1600 muskoxen. Since hunter effort for the winter harvest has remained high and even increased while reported harvests diminished, this supports the 2004 prediction that harvests over 1600 animals would not be sustainable, likely owing to population decline. Alternately, with most muskoxen in high inaccessible terrain in winter, as in March 2018, that would not facilitate hunter success either. Meanwhile, sex-biased harvesting may have occurred. For the period 2011-2013 the commercial and recreational harvests as per særmeldingsskemaerne (special reporting forms, which involve only a portion of the Piniarneq reported harvest) were sex-biased towards cows, i.e., 54.5% cows and 44.2% bulls (Cuyler & Raundrup 2014). Whether this cow-biased harvesting among commercial and recreational hunters continued, or changed has not been investigated, nor whether the sex-biased results from the særmeldingsskemaerne can be applied to the entire harvest, i.e., Appendix 8: Table 17. Obviously, if annual harvests repeatedly took more cows than bulls then population decline could be expected. This century has seen a rapid expansion of local and international demand for qiviut and qiviut products. This motivation likely played a role in raising numbers harvested since 2000, specifically because prices paid for raw skins have skyrocketed and are higher than ever before. By 2017, the increasing worldwide demand for qiviut resulted in international buyers willing to pay Greenland hunters up to DKK 9.500,00 ($1,500 US / 1,280 EUR) for one winter skin (Fernando Alvarez pers. comm.). Since 2014, intact skinless carcasses have been discovered despite attempts to hide them (Cuyler unpublished 2014 field report). Illegal skin-only harvest is an additional mortality of unknown magnitude, which might contribute to possible population decline. Regarding Maniitsoq muskox calf percentage, the 2018 value of 18.4% (albeit a rough approximation) means they are less productive than before 2000, when calf percentage typically was above 25% (Roby 1978, Thing et al. 1984, Olesen 1993, Cuyler & Witting 55 2004, Appendix 9). Already in the 2000-2004 period, calf percentage decreased from ca. 26% to ca. 18% (Appendix 1, Table 13). In that period, muskox densities were often exceedingly high in the lowland elevations, which suggests negative densitydependent factors could have been involved. It may be coincidence that over the same period reported harvest more than doubled from 716 to 1614 muskoxen annually (Appendix 8, Table 17) and muskoxen were noted to move into inaccessible high elevations once the winter hunting season began. Thereafter, excepting 2005, calf percentages continued to decline, 15.5% in 2006, 11.0% in 2009 and ca. 12% in 2010. Again, harvests were heavy, but now muskox densities had fallen i.e., densitydependent factors would have had less affect. Calf percentages from 14.6% to 16.5% in North American caribou populations do not permit population size increase because that level of calf recruitment equals adult mortality, while below those decline is inevitable (Bergerud et al. 2008). Meanwhile for muskoxen, calf percentages of 17-24% can facilitate population growth (Jingfors & Klein 1982, Gunn et al. 1984). The 2018 Maniitsoq calf percentage of 18.4% is in the low end of that range, specifically given the absence of large predators. Predation is not pressing calf production down for Maniitsoq muskoxen. Under current conditions, interpreting whether a value of 18.4% calves is sufficient for population stability or growth is difficult. Future population trend is not obvious. We suspect that a factor behind low calf production is poor cow body condition. Muskoxen are generally not as productive as caribou, and the importance of cow body condition begins already prior to the late summer rut. Muskox cows require 22% body fat to have a 50% probability of pregnancy, while caribou need only 7% body fat (Crête et al. 1993, Adamczewski et al. 1998, Pachkowski et al. 2013). This is important because the probability of successful breeding during the rut increases with the body mass of cows (Rowell et al. 1997, White et al. 1997). Thus, muskox cow pregnancy rates are sensitive to nutritional influences, which if poor lead to reproduction declines (Adamczewski & Flood 1997, White et al. 1997, Adamczewski et al. 1998). Everything that makes vegetation/forage inaccessible, ultimately reduces cow body reserves and thus calf production. For example, winter hunting activities in the vegetation rich lowlands, appear to cause muskoxen to forage at high elevation where vegetation quantity and quality are low. Coinciding with late gestation, cow body condition is likely negatively affected and may be reflected in lower calf survival. With factors(s) causing current 18.4% calves unresolved, the chance for that value to facilitate stability or growth for the Maniitsoq muskox population is debatable. Circumstances for winter harvesting in the past decade may have resulted in a hunting pressure that was not compatible with sustainable use of this renewable resource. The 56 combined results from past counts, densities, harvests, and calf percentages present the possibility of some muskox population decline already having occurred. If true, the Maniitsoq muskox population may in the past have been larger than the current 2018 estimates. There is no indication that abundance grew since 2010. We suggest future population decline is possible and may already be in progress. Given the above, a future stable or growing Maniitsoq muskox population is not certain. This report’s population size estimates, specifically for Maniitsoq muskoxen, are larger than any previous estimate for muskox populations anywhere in Greenland, while the Maniitsoq muskox DS density is much higher than elsewhere in the Arctic. If considering only elevations <400 m that density becomes even greater. It is well known that high animal density can increase exposure of individuals to infectious pathogens (diseases, parasites). Thus, population growth is perhaps not to be recommended given possible density-dependent influences at current population size. Whether the 2018 Maniitsoq muskox population size (ca. 19.000) and density (ca. 2.65/km2) for the Maniitsoq muskox harvest management area (Angujaartorfiup sub-area) are within the current herbivore carrying-capacity of the pasture/range remains to be seen. Acknowledgements This project was financed primarily by the Greenland Government and otherwise by the Greenland Institute of Natural Resources, Nuuk Greenland. Grateful thanks go to Air Greenland Charter and their helicopter pilot Kåre Berli for his safe flying. Thanks also to the Greenland Association of Professional Hunters (KNAPK) for providing an experienced observer, excellent at spotting animals despite poor detection conditions. We thank Katrine Raundrup and Rikke Guldborg Hansen for review of the manuscript. Literature cited Anderson M. & Fergusen M. 2016. Movements and habitat use of muskoxen (Ovibos moschatus) on Barthurst, Cornwallis and Devon Islands, 2003.2006. Status Report 2016-08, Nunavut Department of Environment, Wildlife Research Section, Igloolik, NU. 112 pp. Adamczewski J.Z. & Flood P.F. 1997. Seasonal patterns in body composition and reproduction of female muskoxen (Ovibos moschatus). J. Zool. Lond. 242: 245-269. Adamczewski J.Z., Fargey P.J., Laarveld B., Gunn A. & Flood P.F. 1998. The influence of fatness on the likelihood of early-winter pregnancy in muskoxen (Ovibos moschatus). Theriogenology. 50: 605-614. Bergerud A.T., Luttich S.N. & Camps L. 2008. The return of caribou to Ungava. McGill-Queen’s University Press, Montreal & Kingston / London / Ithica. 586 pp. 57 Beumer L.T, van Beest F.M., Stelvig M. & Schmidt N.M. 2019. Spatiotemporal dynamics in habitat suitability of a large Arctic herbivore: Environmental heterogeneity is key to a sedentary lifestyle. Global Ecology and Conservation 18 (2019) e00647. http://doi.org/10.1016/j.gecco.2019.e00647 Brewer M.J., Butler A. & Cooksley S.L. 2016. The relative perperformance AIC, AICc and BIC in the presence of unobserved heterogeneity. Methods in Ecology and Evolution 7(6): 679–692. Buckland S.T. 1992. Fitting density functions with polynomials. Applied Statistics 41: 63-76. Buckland S.T., Anderson D.R., Burnham K.P. & Laake J.L. 1993. Distance Sampling: Estimating Abundance of Biological Populations. Springer. Buckland S.T., Anderson D.R., Burnham K.P., Laake J.L., Borchers D.L. & Thomas L. 2001. Introduction to Distance Sampling. Oxford: Oxford University Press. Buckland S.T., Anderson D.R., Burnham K.P., Laake J.L., Borchers D.L. & Thomas L. 2004. Advanced Distance Sampling: Estimating abundance of biological populations. Oxford University Press. Buckland S.T., Rexstad E.A., Marques T.A. & Oedekoven C.S. 2015. Distance Sampling: Methods and Applications. Springer. Correia I.J.F. 2020. Estimating caribou abundance in West Greenland using distance sampling methods. MSc. Thesis. University of Lisbon, Portugal. 63 pp. Couturier S., Dale A., Wood B. & Snook J. 2018. Results of a Spring 2017 aerial survey of the Torngat Mountains Caribou Herd. Technical report, Torngat Wildlife, Plants and Fisheries Secretariat . Crête M., Huot J., Nault R. & Patenaude R. 1993. Reproduction, growth and body composition of Rivière George caribou in captivity. Arctic 46(3): 189-196. DOI: 10.14430/arctic1343 Cuyler C. 2020. Biologisk rådgivning for rensdyrog moskusoksefangst i 2020. Advisory document prepared for the Directorate for Fishing, Hunting and Agriculture, 15 May 2020. Pinngortitaleriffik – Greenland Institute of Natural Resources, Nuuk, J. No. GN 40-59/40-00-02-42. 4 pp (Supplementary materials 26 pp). Cuyler C., Landa A., Witting, L., Rosing M., Linnell, J. & Loison, A. 2001. Rådgivning for moskusoksebestanden i Kangerlussuaq, region Nord/Avannaa i 2002, 2003 og 2004. Rapport til Direktoratet for Miljø og Natur. 13 November 2001. Greenland Institute of Natural Resources. J. No. 28.21.21./brev nr. 01168. 14 pp. Cuyler C., Marques T.A., Correia I.J.F. Afonso B.C., Jensen A., Hegelund P. & Wagnholt J. 2021. 2018 status Kangerlussuaq-Sisimiut caribou, West Greenland. Pinngortitaleriffik – Greenland Institute of Natural Resources. Technical Report No. 117. 79 pp. Cuyler C., Nymand J., Jensen A. & Mølgaard H.S. 2016. 2012 status of two West Greenland caribou populations, 1) Ameralik, 2) Qeqertarsuatsiaat. Greenland Institute of Natural Resources Technical Report No. 98, 179 pp. Cuyler C. & Raundrup K. 2014. Svar på spørgsmål angående rensdyrog moskusoksefangst 2014/2015. Advisory document prepared for the Directorate for Fishing, Hunting and Agriculture. Greenland Institute of Natural Resources, Nuuk. 02 May, 2014. J. Nr. 40.00.01.42/14, 3 pp. Cuyler L.C., Rosing M., Egede J., Heinrich R. & Mølgaard H. 2005. Status of two West Greenland caribou populations; 1) Akia-Maniitsoq, 2) Kangerlussuaq-Sisimiut. Pinngortitaleriffik – Greenland Institute of Natural Resources. Technical Report No. 61. Part I-II, 64+44 pp. Cuyler C., Rosing M., Linnell J.D.C., Loison A., Ingerslev T. & Landa A. 2002. Status of the Kangerlussuaq-Sisimiut caribou population (Rangifer tarandus groenlandicus) in 2000, West 64 𝑃 described above assumes that every individual of interest is detected (Miller et al. 2016). Frequently, this assumption cannot be met, specifically if among the individuals of interest there are animals impossible to observe owing to low sightability. Several factors cause low sightability, including topographical barriers, weather conditions, ground surface conditions and many others related to observer training and survey design. The proportion of individuals that were not detected can be estimated using the detection function fitted to the observed distances (Thomas et al. 2002). Once this proportion is estimated, it can be considered to obtain more accurate estimates and then, an extrapolation for a wider region can be done similarly as shown in Equation (2). In Distance Sampling, this proportion of detected objects in the area 𝑎 is defined as the probability of detection, 𝑃𝑎. Therefore, a density estimate can be obtained as per Equation (1) by adjusting 𝑛𝑝𝑙𝑜𝑡 by 𝑃𝑎, i.e., by correcting the detections for those that were missed. Since the latter cannot be known, in general, an estimate must be also obtained, thus Equation (3) where 𝑃 󰆹 a is an estimate of 𝑃 𝑎 obtained from the distance data, and 𝑎 is the area of the sampled region. Usually 𝑎 = 2𝑤𝐿, with 𝑤 as the truncation distance, for both sides of the centreline, and the total transect length 𝐿 = ∑𝑙 𝑗 𝑘 𝑗=1 , where l is the length of transect j. Abundance can be determined using a reasoning analogous to that above (Equation 2). The truncation distance is defined as the distance beyond which distances are not recorded. This can be defined in the field or at the analysis stage. The coefficient of variation of 𝐷 , c 𝑣 (𝐷 ), is related with two random components referred above, encounter rate ( 𝑛 𝑝𝑙𝑜𝑡 /𝐿 ), and 𝑃 󰆹 a , plus a third one that is the estimate of the expected size of detected clusters ( 𝐸 󰆹 ( 𝑠) ). Assuming independence between there, the former is given by Equation (4) An approximation of the standard error of 𝐷 , 𝑠𝑒( 𝐷  ), is defined as Equation (5) 65 Once these are obtained, an approximate 100(1 − 𝛼)% confidence interval (CI) can be determined by Equation (6) Where is the quantile of the N(0,1) distribution 1.96 for a 95% confidence interval). However, the distribution of the 𝐷 is positively skewed, thus an interval assuming that 𝐷  is log-normally distributed has better coverage. According with Buckland et al. (2015), a 100(1-alpha)% confidence interval can be given by Equation (7) where Equation (8) and Equation (9) For further details see Buckland et al. (2001) and Buckland et al. (2015). Probability of detection Given the above, the probability of detecting an object, giving that it is within the area co vered b y the transects, 𝑃  a , needs to b e estimated. F or this pro ject, the ob ject of in terest consists in muskox groups. To illustrate the importance of this probability, consider that an observer walks across a large patch of tundra and detects 8 muskoxen (Fig. 24). While discussing with the local biologist, and considering the biologist’s experience, he/she will state that, on average, only one third of all muskoxen present are detected (i.e. , 𝑃  a = 1/3) meaning that probably there were around 24 muskoxen within that patch of tundra and 16 have been missed. That is where Distance Sampling is useful, since it allows a rigorous framework for the estimation of P a and then an estimate of abundance can be obtained as shown in Equation (3). Distance Sampling methods The detection function, 𝑔(𝑦), describes the probability of detecting an object of interest given that it is at a distance 𝑦, from the centreline (also known as 0-line), thus being a non-increasing function of 𝑦 (Buckland et al. 2015). 66 Figure 24. Example of a patch of tundra with the transect in the middle. Blue dots represent eight observed muskoxen, while orange dots represent the 16 undetected ones. The lines perpendicular to the transect represent the recorded distances. For line transects, 𝑦 is the perpendicular distance from the 0-line to the detected object. Within Distance Sampling methods, the probability of detection is explained recurring to these observed distances (Buckland et al. 2001). Sometimes covariates may be added to explain their relationship with the detection probability. In this situation, we are within the Multiple Covariate Distance Sampling (MCDS) framework (Buckland et al. 2001). Conventional Distance Sampling Conventional Distance Sampling (CDS) occurs when no additional covariates are added to the model. Once the detection function is estimated , 𝑃  a can be obtained via the following equation Equation (10) where 𝜋(𝑦)= 1 𝜔 and, therefore, used to estimate density using Equation (3). For 𝑔(𝑦) it is also specified a flexible semi-parametric model, composed by a key function and some additional series expansions, known as adjustment terms, and their parameters are estimated (Marques et al. 2007). 67 To obtain robust estimates of density, flexible models for 𝑔(𝑦) are needed with the form (Buckland et al. 2001) Equation (11) where 𝑘(𝑦) is the parametric key function and 𝑠(𝑦) represents the additional adjustment terms (Table 15). Table 15. Commonly used key functions and series expansions for the detection function. Adapted from Buckland (2001). The uniform key function has no parameters, while the half-normal and the hazard-rate functions include a scale parameter, 𝜎, which determines the rate at which the function decreases with increasing distance (Fig. 25). Furthermore, the hazard-rate function also includes a shape parameter , 𝑏 , that provides greater flexibility to this function comparing to the others (Buckland et al. 2001). It is not always necessary to include adjustment terms, and in such cases, these models are referred to as “key only” models. When the key functions are not enough for fitting 𝑔(𝑦) , some series expansions terms may be added to modify its shape (Fig. 26). These terms can be either cosine, simple polynomial or Hermite polynomial (Table 15). It is important to note that these adjustment terms do not depend directly on 𝑦 but on 𝑦𝑠 which is a scaled value of 𝑦 , where 𝑦 𝑠 = 𝑦 𝜔 with 𝜔 being the truncation distance. This allows independence between the shape of the series expansion and the units used for 𝑦 (Marques et al. 2007). 68 Figure 25. Half-normal (top row) and hazard-rate (bottom row) detection functions without adjustments, varying scale (σ) and, only for hazard-rate, shape (b) parameters. Values tested are presented above the plots. On the top row from left to right, the study species becomes more detectable (higher probability of detection at larger distances). The bottom rows show the hazard-rate model’s more pronounced shoulder. Adapted from Buckland et al. (2001). 69 Figure 26. Possible shapes for the detection function when cosine adjustments are included for half-normal and hazard-rate models. Adapted from Buckland et al. (2001). 70 Right truncation of the data, or the removal of the largest distances, is a common procedure that aids model fitting. Some precision might be lost with truncation; however, it is usually slight. On the other hand, precision is increased since the data is easier to model and, consequently, fewer parameters and adjustment terms are required to model the detection function (Couturier et al. 2018). Multiple Covariate Distance Sampling CDS methods can be extended to MCDS, so that 𝑔(𝑦) is modelled as a function not only of distance, but also of a vector of 𝐽 additional covariates for each of the 𝑛 objects of interest, z i = z 𝑖1 , ..., z 𝑖𝐽 , 𝑖 = 1, ..., 𝑛 . Accordingly, the function that describes the probability of detection at a given distance, is represented by 𝑔(𝑦, z ) . These additional covariates can either be discrete or continuous, such as observer and group size, and are assumed to affect only the scale, 𝜎, of the detection function (Marques et al. 2007; Miller et al. 2016). For line transects, 𝑃 ( z i) , i.e., the probability of detecting the 𝑖 -th object of interest given its respective vector of covariates z i can be estimated using the formula presented in Equation (12). Equation (12) with 𝜋(𝑦)= 1 𝜔. Considering the three key functions previously presented, only the uniform key is excluded from MCDS since it does not have a scale parameter. Halfnormal and hazard-rate functions can have their scale parameter written as a function of the covariate values as Equation (13) Where 𝛽0and all the 𝛽𝑗’s are the J + 1 coefficients to be estimated with J being the total number of covariates. The estimation of the parameters for both CDS and MCDS is typically done via maximum likelihood (Marques et al. 2007). Once the detection function is estimated, according with (Buckland et al. 2004), density can be estimated as Equation (14) where 𝑎 is the total area surveyed , 𝑃 ( z i ) is the estimated probability of detecting the 𝑖 -th object of interest given its respective vector of covariates zi. 71 Finally, Marques et al. (2007) states that MCDS methods potentially offer improved inference in four situations, when comparing to CDS methods: 1. when a subset of data is used to estimate density, e.g., by strata, where this information can be introduced as a factor covariate. In CDS, the strategy is more complex, either to estimate 𝑃 𝑎 for each stratum and thus, stratum-level estimates for density or to use a global estimate for the probability of detection, but this second introduces bias, for example, if one stratum favours the animals when compared to other strata which uses fewer parameters than a fully stratified detection function model; 2. where pooling robustness does not hold for CDS analyses, e.g., when survey intensity varies according with pre-defined strata to increase efficiency, or when the detection probability faces extreme heterogeneity due to different object habitats or behaviors, for example, showy males contrasting with cryptic females in animal surveys; 3. reduces the variance of density estimates by modelling the heterogeneity in the detection function; 4. if there are covariates of interest to be included in the model. Model selection Since the estimator of density is closely linked to the detection function, it is of critical importance to select models for the detection function carefully. Three properties desired for a model for 𝑔(𝑦) are, in order of importance, model robustness, a shape criterion and estimator efficiency (Buckland et al. 2001, 2015; Miller et al. 2016). The most important property of a model for the detection function is model robustness. According with Buckland et al. (2001, 2015), this means that the model is a general, flexible function that can take a variety of plausible shapes for the detection function. The concept of pooling robustness is also included here. Models of 𝑔(𝑦) are pooling robust if the data can be pooled over many factors that affect detection probability and still yield a reliable estimate of density. A model is pooling robust if, for example, a stratified estimation for density, 𝐷  st , and a pooled estimation for density , 𝐷  p, are approximately the same. In the first scenario, the data is stratified by factors, such as observer or habitat type, and an estimate for density in each stratum is made. Then these estimates are combined into 𝐷  av , an average density estimate. In the second scenario, all data could be pooled, regardless of any stratification, and a single estimate computed, 𝐷  p. A model is pooling robust if 𝐷  av ≈ 𝐷 p. 72 According to Buckland et al. (2001), the shape criterion consists in the fact that the detection function should have a ‘shoulder’ near the line (Fig. 27), i.e., detection remains nearly certain at small distances from the sampling unit’s track line (𝑔 ′ (0) = 0) . This allows the reliable estimation of object density (Thomas et al. 2002). Generally, good models for 𝑔(𝑦) will satisfy the shape criterion near the zero-distance track line, which is especially important in the analysis of data where some heaping at zero distance is suspected. Figure 27. A good model for the detection function should have a shoulder, with probability of detection staying at or close to one at short distances from the centreline or point. At larger distances, it should fall away smoothly. The truncation distance 𝜔 corresponds to the strip half-width (for Line Transect Distance Sampling). Adapted from Buckland et al. (2001). Estimator efficiency is the third most important property (Buckland et al. 2001), which means that it is desirable to select a model that provides estimates that are relatively precise, i.e., that have small variance. This property is of benefit only for models that are model robust and have a shoulder near zero distance, otherwise the estimation might be precise but biased. Besides these three criteria, the model should be a monotonic function of distance from the line, that is, the probability of detection at a given distance cannot be greater than the probability of detection at any smaller distance (Fig. 27) (Buckland et al. 2001). There is no fixed standard method to select the best fitting model, i.e., choosing the most appropriate key function and series expansion (Marques et al. 2007). It is usually done by applying the Akaike’s Information Criterion (AIC), Kolmogorov-Smirnov test, Cramér-von Mises test and the 𝜒2 Goodness-of-Fit test (GOF test). The likelihood ratio test can also be used but, since it is only applicable for nested models, AIC is the recommended method (Marques et al. 2007). A proper model should be simple with an adequate fit without overfitting the data. 73 Akaike Information Criterion The relative fit of alternative models may be evaluated recurring to AIC, or AICc, in case of small samples, providing a small sample bias correction (Buckland et al. 2001). These criteria can be determined as follows Equation (15) Equation (16) where ℒ is the likelihood function, 𝑞 is the number of estimated parameters in the model, and 𝑛 is the sample size. This measure provides a trade-off between bias and variance. AIC includes two terms, one related with the fitted model, and the other working as a penalty considering the excess of parameters in the model (Brewer et al. 2016). Kolmogorov-Smirnov test The Kolmogorov–Smirnov test is one of the tests that can be applied to the detection function to assess model fit (Buckland et al. 2004). This test is only applicable for continuous data, being preferable to the 𝜒2 GOF test for MCDS methods. Considering the cumulative distribution function (c.d.f.) 𝐹 (𝑥) = 𝑃 (𝑋 ≤ 𝑥) and the empirical c.d.f. (e.d.f.) 𝑆(𝑥) , the null hypothesis to be tested is 𝐻 0 ∶ 𝐹 (𝑥) = 𝐹 0 (𝑥), ∀𝑥 . The alternative hypothesis states that both functions differ for at least some value of 𝑥 . In practice, 𝐹 (𝑥) is replaced by its estimate, and 𝐻0 states that the assumed model is the true model for the data (Buckland et al. 2004). The largest absolute difference between 𝐹 󰆹 (𝑥) and 𝑆 (𝑥) , denoted 𝐷𝑛, is the test statistic (Gibbons & Chakraborti 2011). The corresponding 𝑝-value can be approximated by Equation (17) Cramér-von Mises test Similar to the Kolmogorov-Smirnov test, the Cramér-von Mises test shares the same null hypothesis and basis on differences between c.d.f. and e.d.f. However, instead of considering only the largest difference between the two functions, this test is based on their entire range (Buckland et al. 2004). The test statistic can be given by Equation (18) 80 Figure 30. Group of 23 muskoxen, which may have included 1-2 calves (age ca. 10-months). Figure 31. Group of 38 muskoxen, which may have included 7-8 calves (age ca. 10-months). Eight adults kept themselves somewhat separated from the main group. 81 Appendix 6 Maps illustrating elevation, aspect, and slope in the North region Elevation Figure 32. Map illustrating elevation in the surveyed area. Legend colour codes the elevation. Aspect Figure 33. Map illustrating aspect in the surveyed area. Legend colour codes the compass headings for aspect. 82 Slope Figure 34. Map illustrating slope in the surveyed area. Legend colour codes the angle of slope. 83 Appendix 7 Angujaartorfiup sub-area, March 2018: Photos of survey conditions and snow cover observed. Photos by C. Cuyler Figure 35. View from surface of lake ice during a pause between survey lines, illustrating lack of snow common in the Angujaartorfiup sub-area. Also illustrating flat-light (no shadows). Figure 36. View as helicopter flew over surface of lake ice, illustrating the typical sparse snow cover on the ground ahead, line transect 18 Angujaartorfiup sub-area. Also illustrating flat-light (no shadows). 84 Figure 37. View while flying line transect 17, illustrating sparse snow-cover common to the Angujaartorfiup sub-area. Also illustrating flat-light (no shadows). Figure 38. Overview of the north side of Tasersuaq Lake (large flat white surface on left: lake elevation 312 m), illustrating flat-light (no shadows) and sparse snow-cover common to the Angujaartorfiup sub-area. 85 Figure 39. Overviews illustrating flat-light (no shadows) and the windswept terrain with sparse snow-cover common to the Angujaartorfiup sub-area. Horizontal white surfaces are lake or fjord ice. 86 Figure 40. Minimal March snow cover typical of lowland valleys in Angujaartorfiup sub-area. View NW into the Maniitsut Atannguisa Kuussuat river valley from line transect 19. Valley bottom elevation <100 m. Figure 41. Observation platform for aerial survey: AS350 helicopter on terminal moraine with Greenland Ice Cap in background, illustrating windswept rocky terrain with sparse snow-cover. 87 Appendix 8 Hunting seasons, quotas and reported harvests, primarily for the Maniitsoq muskox management area. Table 16. Muskox hunting seasons and quotas for the Sisimiut and Maniitsoq muskox management areas in the 2015-2020 period (commercial and recreational combined) (Nuka M. Lund pers comm (Ministry for Fisheries & Hunting: APN). Trophy hunting seasons and quotas not included. Maniitsoq muskox management area corresponds to the Angujaartorfiup sub-area surveyed by helicopter in March 2018, and the Sisimiut management area corresponds to Sisimiut sub-area surveyed. Muskox management area Hunting area1 WINTER HUNT AUTUMN HUNT Season Quota Season Quota Maniitsoq 2015 1 10 January – 10 March Open 1 August – 15 October Open2 2 Closed 0 Closed 0 3 10 January – 10 March Open 1 August – 15 October Open 2016 1 10 January – 10 March Open 1 August – 15 October Open 2 Closed 0 Closed 0 3 10 January – 10 March Open 1 August – 15 October Open 2017 1 10 January – 10 March Open 1 August – 30 September Open 2 22 January – 31 January 400 Closed 0 3 10 January – 10 March 800 1 August – 30 September Open 2018 1 10 January – 15 February Open 1 August – 15 October Open 2 Closed 0 Closed 0 3 10 January – 15 February 800 1 August – 15 October Open 2019 1 25 January – 15 February Open 1 August – 15 October Open 2 Closed 0 Closed 0 3 25 January – 15 February Open 1 August – 15 October Open 4 25 January – 15 February Open 1 August – 15 October Open 2020 1 25 January – 14 February 100 1 August – 15 October Open 2 Closed 0 1 August – 15 October Open 3 25 January – 14 February 1200 1 August – 15 October Open 4 25 January – 14 February 1 August – 15 October Open Sisimiut 2015 --- 10 January – 10 March 520 1 August – 15 October 400 2016 --- 10 January – 10 March 520 1 August – 15 October 400 2017 --- Closed 0 1 August – 15 October 400 2018 --- Closed 0 1 August – 15 October 400 2019 --- Closed 0 1 August – 15 October 400 2020 --- Closed 0 1 August – 15 October 400 1 See Figures 40-41 for the Greenland government’s muskox management hunting areas for Maniitsoq (Angujaartorfiup sub-area surveyed). 2 Open = No limit on the number of muskoxen that may be harvested. 88 Figure 42. Map for the 2010-2019 period illustrating the Government of Greenland’s Maniitsoq muskox harvest management area with specific hunting areas for regulating harvest of the Maniitsoq muskox population. Hunting areas are separated by stippled lines and labelled with “Piniarfik/Jagtområde” and a number. Areas where hunting is prohibited are marked with diagonal lines and/or the label, “Piniarfigeqqusaangitsoq/Jagtfrit område”. The largest encompasses the airport town of Kangerlussuaq and near vicinity including Lake Ferguson, the beginnings of Ørkendalen and the Sandflugtdalen river valley up to and including the Inland Ice Cap. The smaller and triangular area encompasses the inner portion of the valley, Arnangernup Kuua (Paradisdalen) and its side valleys (Naalakkersuisut 2019). 89 Figure 43. Map illustrating the entire Maniitsoq muskox harvest management area with the current hunting sub-areas for regulating harvest of the Maniitsoq muskox population in the Qeqqata municipality. The hunting sub-areas are defined by the Government of Greenland. Hunting areas are labelled with “Piniarfik” and a number. Hunting area 4, established in 2018, was first implemented from summer 2019 (Naalakkersuisut 2019, bilag 4). Areas where hunting is prohibited are either marked with diagonal lines and the label, “Piniarfigeqqusaangitsoq”, or are a solid red colour. Since 1984, hunting has been prohibited in the inner portion of the valley, Arnangernup Kuua (Paradisdalen), which is indicated by the solid red within Piniarfik 4. (Naalakkersuisut 2020). 96 Figure 3. Muskox populations in Greenland. Native populations (blue) are in North and Northeast Greenland. Inglefield Land muskoxen are a mix of native and translocated stock. In western Greenland, specifically from Cape Atholl to Nanortalik, all populations are either translocated (orange) or the result of expansion (green) from a neighbouring translocation. Population borders are approximations and size does not reflect animal abundance, e.g., Greenland’s largest muskox population, Maniitsoq, inhabits a relatively small area. Boundary of the National Park in northeast Greenland is indicated. Currently, regulated harvest is permitted on most populations, however hunting is prohibited in Washington Land, National Park, Nanortalik and Nuuk. 97 Maniitsoq muskox population In 1963 and 1965, 27 juvenile muskoxen were translocated from East Greenland to the area south of the Kangerlussuaq (Søndre Strømfjord) International airport in West Greenland. Previously, this population was commonly called the Kangerlussuaq muskoxen (Cuyler & Witting 2004), but hereafter will be referred to as the Maniitsoq muskoxen as per APNN (2019a). All hunting was prohibited for 25 years, first harvest was 1988, and by 1990 the muskox population was considered well established (Olesen 1993). Already in the 1980’s, the muskoxen clearly selected lowland elevations under 400 m. At that time, this included the lowlands surrounding the Kangerlussuaq (Søndre Strømfjord) international airport and several large valleys, i.e., Ørkendalen (Bioshytte), Ammalortoq Lake, and Paradisedalen (Arnangarnup Qoorua) (Olesen 1993). The total area available to Maniitsoq muskoxen is 7,853 km2, (corrected area from this study is 7133 km2), however only ca. 950 km2 is habitat under 400 m elevation (Olesen 1993). Calf production indicated that lowland habitat within the 7,853 km2 region was well suited to muskoxen. Although neither is normal for muskox populations elsewhere, in the 1970-80’s, calves accounted for 24-28% of the Maniitsoq muskox population, while many cows produced a calf every year (Roby 1978, Thing et al. 1984, 1987, Olesen 1993). To manage the Maniitsoq muskox harvest, APNN subdivided the region into four hunting management areas. Initially (2010) there were three hunting management areas (Piniarfik 1, 2, 3) but a fourth was added in winter 2018 and first implemented in summer 2019 (Figure 4). Figure 4. APNN’s four hunting areas (Piniarfik 1, 2, 3, 4) for harvest management of the Maniitsoq muskox population in Qeqqata municipality (APNN 2019b). Top right and indicated by red diagonal lines is the huntingprohibited area (Piniarfigeqqusaanngitsoq) that surrounds the Kangerlussuaq (Søndre Strømfjord) international airport and includes Sandflugtdalen Valley and Isunngua areas. Hunting has always been prohibited in the inner Paradise Valley, indicated by solid red. 98 Current population trend, Maniitsoq muskoxen A 2018 helicopter survey provided a population size estimate of ca. 20,334 muskoxen (95% CI 9,386 – 44,055; SE 6397; CV 0.31) for the Maniitsoq muskox population (Marques 2018). This number cannot be compared with previous ground-based minimum counts, because it is the first and only estimation of total abundance of muskoxen in the entire area south of the Kangerlussuaq international airport. The high 2018 estimate provides evidence of the explosive increasing trend for abundance that occurred since 27 individuals were translocated to the area in the 1960’s. The steep growth curve in the 1980’s and ‘90’s (Figure 5) supports this. Regardless, alone the 2018 estimate cannot indicate current population trend. Instead, past, recent, and current trends are indicated by the minimum count and calf percentage data gathered over a 40-year period (Figures 5, 6). These suggest the current trend is steady decline in abundance since 2010. The most recent minimum count, in 2020, also supports that population size and production have decreased. The abundance in the 2000-2010 period might have been greater than the 2018 aerial survey estimate of ca. 20,334 muskoxen. Below follows a brief discussion of the minimum counts, calf production, 2018 aerial survey, muskox density, and disturbance. Minimum counts 1963 – 2020, Maniitsoq muskoxen Since the initial release of 27 muskoxen in 1963-65, ground-based minimum counts from the 1980’s to 2020 provide indices of abundance. Although not population size estimates, they illustrate general population trends (Figure 5). Minimum counts in this area are usually done by two observers travelling by skidoo in the January-April period and using binoculars and telescopes recording all muskoxen seen. The 2018 minimum count was the actual number of muskoxen seen during a helicopter survey for caribou. Certainly, the inconsistent timing, effort and coverage of the diverse minimum counts makes comparisons across time less reliable. For example, the 1977-1995 minimum counts were carried out before the hunting season and covered hunting areas 1 and 2 (see Figure 4). Minimum counts in 2000-2010 occurred also prior to winter harvesting, but typically covered most of hunting areas 1, 2, 4 and northern portion of area 3. In 2010, lack of snow and sea ice limited the minimum count to hunting areas 1 and 2 and a small portion of 3, while 4 was impossible to cover. Meanwhile, winter harvest influenced minimum count results post-2014 because counts were during or after the winter hunting season. During these counts, animals had probably moved to less accessible terrain, which prevented their detection. Furthermore, 2014-2020 area coverage was highly variable and often much reduced. The exception being the 2018 minimum count from helicopter, which had the greatest effort and coverage of any count to date. Owing to different timing and effort, muskox abundance trends should be considered with caution. Nevertheless, general trends are apparent. The index of muskox abundance illustrates a general rapid expansion in the 1980’s and early 1990’s, a possible peak from the mid 1990’s until 2010, and a decline thereafter (Figure 5). Neither the magnitude nor slope of the decline are likely as steep as they appear, since 99 minimum counts after 2014 occurred after hunting had altered the distribution and detectability of the animals. Regardless, the suggested declining trend in population size since 2010 is supported by simultaneous declining number of muskox groups observed, maximum group size, as well as a ca. 37% reduction in calf recruitment (Appendix 2). Furthermore, minimum count data from hunting area 2 (Appendix 3), which received consistent effort and coverage on all ground-based counts, although timing varied, similarly supports an abundance decline in the Maniitsoq population since 2010. Figure 5. Index of Maniitsoq muskox abundance since translocation to the region in 1963-65. Data are minimum counts of only those individuals observed and cannot be interpreted as estimates of population size. All but 2018 are ground counts. All are of varying timing, effort, and area coverage. Beginning 2000, all are pre-calving counts from either winter or late winter. Those from 2000 to 2010 were pre-harvest, while those 2014-2020 were post-harvest. In the latter muskoxen are avoiding lowlands owing to recent hunting. The actual number of individual muskoxen observed during the March 2018 helicopter survey () has been included as a type of minimum count, which effort and coverage were the most comprehensive relative to any other minimum count. (Data from Roby 1978, Thing et al. 1984, Olesen 1993 and Cuyler unpublished). Calf production 1977-2020, Maniitsoq muskoxen Until at least 1981 calf percentage for Maniitsoq muskoxen was excellent (Figure 6) resulting in the highest recorded population growth rate, 35%, for any muskox population globally (Thing 1984, Olesen 1993). Population growth for Maniitsoq muskoxen was above the 1724% reported for a similarly translocated population in Alaska (Jingfors & Klein 1982) and the 21% from Queen Maud Gulf, Canada (Gunn et al. 1984). Maniitsoq’s population increment, ca. 30%, continued almost unabated until 1990 (Olesen 1993), and even into 2000-2001 (Cuyler et al. 2001). Beginning in 2002 however, calf percentages declined and have not recovered to pre-2000 values (Figure 6). Calf percentage seemingly improved since the 2010 minimum count, however, the 2014-2020 period involved only post-harvest counts. Resulting calf percentages may just be an artifact of hunters shooting all except calves, thereby reducing the total number while calf number remained the same, which artificially 0 1000 2000 3000 4000 5000 6000 1950 1960 1970 1980 1990 2000 2010 2020 2030 Number muskoxen observed Year 100 raised calf percentage in the counts. Specifically, the 2014 value is likely unreliable and artificially high because the count was post-harvest, no-hunting zones did not yet exist and the small size of surveyed area, which was easily accessible to hunters. In the 2016-2020 period, a no-hunting zone was established, and minimum count effort was concentrated there. Thus, the ca. 15% calves for this period provides the best current expected calf percentage, which is ca. 40% below the pre-2000 values and supported by the ca. 37% known decline in calf recruitment (calves per 100 cows; Appendix 2). Figure 6. Calf percentages in the Maniitsoq muskox population from 1977 to 2020. All results are from ground counts. Beginning 2000, all results are from winter (pre-calving) counts, 2000-2010 pre-harvest and 2014-2020 post-harvest. The April 2014 outlier value () is likely unrealistic and artificially elevated, primarily owing to the timing of the count, which was post-harvest, and count effort involved a limited area easily accessible to hunters. Thus, while total muskox number was reduced calf number remained the same (not shot), which artificially raised calf percentage. (Data from Roby 1978, Thing et al. 1984, Olesen 1993 and Cuyler unpublished). The 1980’s and 1990’s were associated with rapid population growth, when the percentage of calves was usually ca. 25% (Figure 6) of the total number of muskoxen observed. In 1990, near Kangerlussuaq (Søndre Strømfjord) international airport, Ørkendalen Valley (Bioshytte/Bios cabin) and Ammalortoq Lake, density was almost 2 muskoxen/km2 at elevations under 400 m, however, density was greater in the most commonly used valley areas, and yet there were no observations of density-dependent mortality i.e. no increased calf mortality had been observed even in the areas of highest densities (Olesen 1993). This suggests that for Maniitsoq muskoxen, a density of ca. 2 muskoxen/km2 in lowland habitat was no impediment to calf production. By winter 2001, density grew to at least 4-5 muskoxen/km2 in the same area (Cuyler et al. 2001). Again, because winter grazing was often concentrated, some valley locations had much higher densities, e.g., ca. 29 muskoxen/km2 0 5 10 15 20 25 30 1960 1970 1980 1990 2000 2010 2020 2030 Calf percentage (%) Year 101 near Bios cabin Ørkendalen (589 muskoxen in 20.25 km2) and ca. 22 muskoxen/km2 around Ammalortoq Lake (876 muskoxen in 40.5 km2) (Cuyler et al. 2001). Even at those densities, calf percentage was ca. 25% and annual population growth ca. 30% (Cuyler et al. 2001), which indicates that the carrying capacity of that lowland habitat was excellent. However, something changed in the period 2001-2004, because calf percentage decreased each year and this trend continued until 2010 (Figure 6). Typical causes for poor calf production include predation, disease, adverse weather, and maternal body condition. Elsewhere, large predators (Reynolds et al. 2002, MarquardPetersen 1998) and fatal disease outbreaks (Ytrehus et al. 2008) can negatively affect calf percentages. However, for Maniitsoq muskoxen, large predators are absent and fatal disease has not been observed. This leaves poor cow body condition and adverse weather as primary factors. Given the relative stability of the xeric continental climate of the region, weather conditions or events are not expected to explain the observed continuous decline in calf production in the 2000-2010 period, or the relatively stable but low values thereafter. This makes cow body condition the likely primary factor. Causes of poor cow body condition would include muskox densities too high for the forage quantity, quality, and availability in the region, and excessive disturbance preventing cows from regaining body reserves in summer or maintaining those reserves through winter. Simultaneously with the declining trend for calf percentages, in the period 2001-2004 there occurred a significant (p<0.05) reduction in cow rump fat depth, and pregnancy rates (cows bearing a foetus) fell from 74.6% to 61.5% (Cuyler & Witting 2004). Poor cow body condition (fat reserves) suggests the possibility of density-dependent effects and/or excessive disturbance. Regarding the former, the high densities observed in 2000 and 2001 did not appear to change calf percentage. Regarding disturbance, at least for the 2001-2004 period we know harvests doubled (Cuyler & Witting 2004). It may only have been coincidence that calf % dropped just as harvesting increased, specifically winter harvesting. On the other hand, the increased disturbance occasioned by winter hunting activities may have negatively influenced cow body condition and ultimately calf production. Winter hunting began in 1994 and initially harvest disturbance was minimal. For example, the first six years winter quotas were small, averaging 227 muskoxen, or about half the annual quota. The harvests involved a small number of commercial hunters from only the towns of Sisimiut and Maniitsoq, hunting was only in what are now referred to as hunting areas 1 and 2, and Sisimiut hunters used almost exclusively dog sleds for transport. Maniitsoq hunters had a government granted dispensation to use skidoos. In the period 2000-2010 several things occurred. The Maniitsoq winter harvest quotas began at 500 muskoxen and rose quickly thereafter. The harvest continued to involve a growing but still limited number of commercial hunters, still mostly from Sisimiut and Maniitsoq. Many Sisimiut hunters replaced dog-teams with motorized vehicles. Sport hunters were admitted. Finally, by the end of the period a third hunting area was added, which meant winter hunting now occurred over the entire region, excepting Paradise Valley and the near airport area. In sharp contrast to the early days, since ca. 2010 and continuing today, winter harvest involves large numbers of 102 commercial and sport hunters, using many and diverse motorized vehicles, which penetrate deeply into the region and specifically all lowlands. Meanwhile, trophy harvesting has always been permitted into April, ending just as muskox calving begins. In the beginning the trophy harvest involved few agents and hunters. Like the regular harvesting, however, today the trophy harvest industry is greatly expanded involving multiple actors and concession districts covering much of the region. Muskox calving normally begins in mid-April and lasts to past mid-June (Lent 1988). Disturbing parturient cows in late gestation is not expected to be compatible with good calf production. The current January-February harvest likely disturbs normal winter distribution, as well as foraging and rumination time for muskoxen. Likely this also applies to the MarchApril trophy harvest. Each winter since 2000, hunting appears to cause Maniitsoq muskoxen to avoid lowlands and seek vehicle-inaccessible terrain at relatively high elevations (Cuyler unpublished). Avoiding areas accessible to hunters was apparent during the March 2018 aerial survey (Figure 7). Although 2 weeks after the end of the regular 2018 winter hunting season, most muskoxen had not returned to lowlands, and the hunting of trophy bulls was still on-going at the time. Avoiding essential lowland habitat in winter likely has a negative effect on the winter body condition of pregnant cows. Thus, winter harvest may have influenced the steep decline in calf percentage in the 2000-2010 period (Figure 6). Finally, the 2020 minimum count involved only hunting areas 1 and 2, which contain the 950 km2 under 400 m elevation described in Olesen (1993). The 2020 post-harvest count yielded 16.9% calves, and a density of 0.9 muskoxen/km2 in Olesen’s lowlands (Cuyler & Mølgaard unpublished). Both are well below the 2000 values: 26% and 4-5 muskoxen/km2 respectively (Cuyler et al. 2001). Although elsewhere in the Arctic muskox densities are typically lower, for Maniitsoq muskoxen the relatively high 2000-2001 densities of 4-5 muskoxen/km2 in lowland habitat (under 400 m elevation) did not appear to influence calf percentage, which remained high. This suggests that today’s 2020 post-harvest density of 0.9 muskoxen/km2 is unlikely to be a limiting factor to population growth, prerequisite on that cows can utilize lowland habitat as they did prior to and including 2000-2001. Instead, harvest may be the primary limiting factor, specifically its associated disturbance, which causes muskoxen to avoid lowland habitat, although total catch is likely also important. Improving cow body condition would promote calf production. Strategies would include reducing the duration of human disturbance. Shorter hunting season(s) could increase undisturbed foraging in optimal habitat. A shorter summer hunting season, e.g., beginning after mid-August, would likely facilitate rebuilding cow body condition lost during the winter and due to calving and lactation, resulting in more cows achieving the prerequisite 22% fat for ovulation during the rut. Regarding the winter harvest, a shorter winter season ending well before parturition would likely facilitate successful gestation and birthing. Furthermore, the magnitude of disturbance associated with an otherwise shorter seasons might be minimized if fewer persons and motorized vehicles were involved than is currently normal. 103 2018 survey Maniitsoq muskoxen March 2018 was the first ever systematic transect distance sampling aerial survey for Maniitsoq muskox abundance and distribution in the region south of the Kangerlussuaq international airport in West Greenland. The 2018 population size estimate for Maniitsoq muskoxen was ca. 20,334 muskoxen (95% CI: 9,386 – 44,055; SE = 6,397; CV = 0.31). The large Coefficient of Variance (CV) reflects poor accuracy on the population estimate. This estimate cannot provide population trend, because it was the first using this method, and there are no other estimates for comparison. The March 2018 aerial survey occurred two weeks post-harvest (albeit trophy harvest continued) and muskoxen were still avoiding lowlands and other areas easily accessible to hunters. This is likely the primary factor causing the nonuniform distribution observed across the region (Figure 7). Few muskoxen (n=8, 0.5%) were observed in the southernmost portion of the region near the Sukkertoppen Ice Cap, likely because that area is characterized by barren high elevations of around 1000 m or more. Despite the large area of lowlands in the north only a few muskoxen (n=81, 5.5%) utilized that area. This is remarkable, given that lowlands are the preferred habitat for muskoxen. This area is, however, within easy reach for hunters coming from the Kangerlussuaq international airport town. Many muskoxen (n=200, 13.6%) were in an area halfway between the airport and the Sukkertoppen Ice Cap. Meanwhile, the majority concentrated in two subareas on either side of the region. Most (n=632, 43.1%) were in the east near the Greenland Ice Cap and the rest (n= 546, 37.2%) were in the west. This was not surprising because in the winter 2018, muskox hunting was prohibited in hunting areas 2 and 4, while motor vehicles were prohibited in hunting area 1. These restrictions coincide with the two greatest concentrations of muskoxen post-harvest. Regardless of where they were observed, few muskoxen were utilizing lowlands. Figure 7. March 2018 distribution of the 1,467 muskoxen observed two weeks post-harvest (albeit still trophybull season) during the aerial survey of the Maniitsoq muskox population. Boundaries are approximate. (Cuyler unpublished). 104 Supplementary Materials Cuyler (2020: Appendix 2) Maniitsoq muskox population Population trends, 2000 – 2020 Data from ground-based minimum counts, using skidoo and ATV. Data from 2000-2010 period is pre-harvest, while that from 2014-2020 is post-harvest. Figure 9. Number of groups of Maniitsoq muskoxen observed during minimum counts in the period 2000-2020. Figure 10. Number of groups of Maniitsoq muskoxen observed during minimum counts in the period 20002020. y = -27,769x + 56154 R² = 0,5817 0 100 200 300 400 500 600 700 800 900 1000 1995 2000 2005 2010 2015 2020 2025 Number of groups observed Year y = -0,1115x + 232,48 R² = 0,239 0 2 4 6 8 10 12 14 1995 2000 2005 2010 2015 2020 2025 Average group size of muskoxen Year 105 Figure 11. Maximum Maniitsoq muskox group size observed during minimum counts in the period 2000-2020. Figure 12. Winter calf recruitment as represented by calves (age ca. 10 months) per 100 cows, observed during four demographic counts completed in the period 2006-2020 for Maniitsoq muskoxen. The recruitment values in 2017 and 2019 are lower than those obtained in 2006 and 2010 (Figure 11). Given the former are post-harvest values, although relatively low, these are probably artificially high since cows are commonly harvested but calves are not. Actual calf recruitment in 2017 and 2019 was likely lower than shown. y = -1,9725x + 4018,5 R² = 0,465 0 10 20 30 40 50 60 70 80 90 100 1995 2000 2005 2010 2015 2020 2025 Maximum muskox group size observed Year 0 5 10 15 20 25 30 35 40 45 50 2004 2006 2008 2010 2012 2014 2016 2018 2020 Calves per 100 Cows Year 112 Appendix 10 Recommendations for improving future surveys for muskoxen. Aerial survey methods & design: The approximately 11% survey coverage in 2018 promotes accuracy of abundance estimates and should be continued in future to facilitate population trend. Given muskoxen are stationary and the winter hunting season appears to disturb their distribution, in future the best option may be helicopter surveys specific for muskoxen and flown before the winter hunting season begins. This is because hunting seems to disturb normal distribution, which resulted in ‘hot-spots’ of high concentrations of muskoxen that made for high variance among line transects. To obtain greater accuracy and precision for estimates of muskox abundance and density, then the problem of uneven distribution and lack of presence in lowlands must be rectified. Alternatives might include a later survey period than currently used e.g., late March or early April. Albeit, all winter surveys should end before the onset of calving, which for muskoxen can begin by mid-April. Another alternative would be to prohibit hunting and motor vehicle activity before (>2 weeks) and during the survey period. A one-month period, or more, without hunting or other human disturbance may be necessary before muskox vigilance relaxes, as exhibited by foraging in the preferred lowland elevations and moving into known core winter range in hunting areas 1 and 2. Muskox Demographics: Obtaining accurate demographics was not possible during the Distance Sampling survey. Even at the slow helicopter speed flown while flying line transects, when a group of muskoxen was detected, there was usually only sufficient time to obtain total group size, distance from the 0-line and sometimes the group’s behavioral reaction to the fly-by. Accurate identification of the relatively small-bodied calves (age 9-10 months) was difficult, owing to animals often milling about and calves typically remaining hidden among or behind larger members of the group regardless of their distance from the track line. Muskox herd structure data must be collected in a specific effort separate from flying the line transects for Distance Sampling data. Future budgets must permit a minimum 3 to 6 hours helicopter time devoted exclusively to muskox demographic data collection. Alternately, a separate ground effort using aerial drones could obtain accurate demographics if used prior to hunting season(s) and on core range. Helicopter Logistics: Always check whether other helicopter options are available. To date, the smallest helicopter available is the AS350 and only from Air Charter (Air Greenland). The AS350 permits limited vision for rear observers, owing to the small window size containing several bar/struts, and which under cold ambient temperatures always fog with ice-frost. These factors reduce visibility of terrain. 113 [Empty page]