Legacies of Past Exploitation and Climate affect Mammalian Sexes Differently on the Roof of the World - The Case of Wild Yaks

In polar environments, a lack of empirical knowledge about biodiversity prompts reliance on species distribution models to predict future change, yet these ignore the role of biotic interactions including the role of long past human exploitation. To explore how mammals of extreme elevation respond to glacial recession and past harvest, we combined our fieldwork with remote sensing and used analyses of ~60 expeditions from 1850–1925 to represent baseline conditions for wildlife before heavy exploitation on the Tibetan Plateau. Focusing on endangered wild yaks (Bos mutus), we document female changes in habitat use across time whereupon they increasingly relied on steeper post-glacial terrain, and currently have a 20x greater dependence on winter snow patches than males. Our twin findings—that the sexes of a cold-adapted species respond differently to modern climate forcing and long-past exploitation—indicate that effective conservation planning will require knowledge of the interplay between past and future if we will assure persistence of the region's biodiversity.

Legacy Effects as Inferred from Historical Data. As indicated, legacy effects were operationally defined as those behavioural attributes of a population shaped by biotic interactions over time. While many cases will therefore be evolutionary in nature 1 , in our case we rely on decadal or relatively immediate (i.e. proximate) responses reflected by behavioural plasticity especially those associated with extreme exploitation. The detection of legacy effects will be enhanced best by analyses across multiple generations and comparative study. References from Table 1 are as follows: migration [2][3][4][5] , distributional shifts [6][7] , activity [8][9] , personality [10][11] , and diet selection [12][13] . Other taxa with documented or suspected legacy effects of intense exploitation include cetaceans 14 , primates [15][16] , and marsupials 17 .
To check the possibility of legacy effects in wild yaks, we examined 59 historic expeditions across the Tibetan-Himalayan region to establish baseline values  of habitat use prior to periods of widespread shooting (which began after the 1930s). The first historical account, written by Bogle in 1774 and made available in 1876 20 , was followed by Moorcroft and Trebeck in 1841 21 and included their 1819-1825 journeys to part of Ladakh, Kashmir, and associated regions. The last was Kingdon Ward's (1937) Across Southern Tibet in 1935, which includes the crossing of 3 passes >5,500 m and reference to domestic but not wild yak 76 . In most cases explicit mention of yak presence or absence is recorded: "The wild yak is not found anywhere in this region" (Wollaston 1922:7) 77 . Expeditions were taken primarily during winter because summer travel was inefficient due to rivers and sodden areas 18,19,22 .
Surveyors often reported habitat, groups, or sex. For instance, Prejevalsky (1876), "the southern slope of this range that we saw herds of them" or "I was clambering over the mountains when I suddenly caught sight of three yaks lying down" 33, 50 , Welby (1898) "I spotted in a valley… a single yak grazing" 23 ; and Rawlins 1905 "whilst immense herds of wild yak were to be seen grazing in all the valleys" 22 .
Twenty-eight expeditions yielded 217 descriptions in which yak habitat, group, or sex were inferable. The word nullah was used in 9 of these, 2 for putative female groups and 7 of males in small groups or singly. Because it was unclear where nullahs were situated and because the sample size for nullah usage itself was insufficient for classification as a habitat, these 9 cases were censored. Likewise, of the total 208 remaining accounts, 60 did not include simultaneous mention of both habitat and sex (or its inference), so these were excluded, rendering 148 usable accounts (Fig. S1).
To assess whether groups noted in the historical record were of males or females, we relied on three methods. First, we used direct surveyors' notation and excluded cases in which we could not reasonably make a determination. Second, we based inferences on currently known proclivities of male and female sexually dimorphic ungulates to associate in different-sized groups [78][79][80] . However, this pattern was noted for wild yaks nearly 125 years ago: "the cows are generally to be found in herds varying in numbers from ten to one hundred, while the old bulls are for the most part solitary or in small parties of three or four" (Kinloch 1892:118) 27 , and well before the heavy human exploitation period. The groups we observed fit well within these earlynoted size categories for the two sexes. For instance, in 2012 the largest female group was ~210.
Median group size for females was 23. For males it was 1, and 2/3rds of those consisted of 1-2  80 .
Finally, we assumed groups to be of males when 7 or less, an assumption that we previously assessed by contrasting group sizes between those having been documented with known males and those with presumptive males. If differences in these categories of males existed, the assumption would be invalid. However, means and 95% CI of known and putative male groups To evaluate the probability of resource selection by sex as a function of group size for glacial and primary productivity-related spatial covariates, we used RSF 85 . We applied mixed-effects models to account for potential differences in yak distribution using a random intercept for survey years.

Analytical Approaches and Modeling of Resource Selection.
To test for broad differences in habitat use and to relate it to historical surveys, we first compared frequencies of use between sexes and as a function of group size. We also conducted direct comparison of resources used by males and females, using logistic regression to check the possibility of sex differences. We next evaluated resource selection with a used-available design comparing the covariate values at locations of yaks compared to random locations within a 3-km buffer of the survey route.
This approach was conservative because it effectively assumed equal sightability within our survey buffer although yaksdue to their larger size and flat sparsely vegetated survey areasmay be detectable at greater distances 86 . We generated 1 random location for each 1-km of each survey in each year for comparison to actual yak locations.
For all models, we used scatterplot matrices (i. e. Fig. S2) and a threshold correlation coefficient of r=0.5 to screen for collinearity among spatial covariates; model selection was conducted using Akaike Information Criteria [87][88][89] . In screening for co-linearity amongst spatial covariates with scatterplot matrixes and correlation coefficient of r=0.5, we noted our two different measures of distance to glaciers, their edge and centroid, were too correlated to include in the same model (r>0.7). Distance to glacier edge explained yak resource selection more than distance to center of glacier, so we retained distance to glacier in subsequent modeling. No other covariates were correlated (r>0.5), though distance to glacier edge and glacier centroid were close (r=0.43). Model selection was conducted using Akaike Information Criteria (AIC) using the number of independently observed yak groups as the sample unit. Kfolds cross validation of the top RSF model was performed 87 with all analyses conducted in R 90 .
Latitude:longitude plots for the sexes suggested no skew in a given direction (Fig.   S2). And, while as the effects of sex on resource selection for NDVI, snow and slope were strong ( Fig. 3; Table 2), AIC evidence ratio of averaged group vs sex effects was 9.1 (Table S 1). Thus, we retained group-size effects in the best RSF model, not sex, although both are interchangeable for wild yaks given our sex-group size distributions.
Model fit was not improved by inclusion of a random intercept for survey year (i. e. AIC between the top model in Table S1 with and without the random intercept revealed the fixed-effects model in Table S1 had better model fit). The overall model fit was highly significant with a likelihood ratio test X 2 =6.43, p<0.00001, and a Hosmer-    shown. Parameter definitions are as follows: NDVI = normalized difference vegetation index for the previous summer's growing season (July-August); Snow = percentage snow cover occurring during the survey; Slope = degree's slope of the terrain; DG = distance to nearest glacier edge in km; A_glacier = area of the nearest glacier edge in km 2 ; Grp = group size interactions from Eq. 1; Sex = sexspecific interactions with the same.