Convergent geographic patterns between grizzly bear population genetic structure and Indigenous language groups in coastal British Columbia, Canada

Landscape genetic analyses of wildlife populations can exclude variation in a broad suite of potential spatiotemporal correlates, including consideration of how such variation might have similarly influenced people over time. Grizzly bear (Ursus arctos) populations in what is now known as coastal British Columbia, Canada, provide an opportunity to examine the possible effects of a complex set of landscape and human influences on genetic structure. In this collaboration among the Nuxalk, Haíɫzaqv, Kitasoo/ Xai’xais, Gitga’at, and Wuikinuxv First Nations and conservation scientists, we characterized patterns of genetic differentiation in the grizzly bear, a species of high cultural value, by genotyping 22 microsatellite loci from noninvasively collected hair samples over a 23,500 km2 area. We identified three well-differentiated groups. Resistance surfaces, which incorporated past and present human use, settlement, and landscape resistant features, could not explain this pattern of genetic variation. Notably, however, we detected spatial alignment between Indigenous language families and grizzly bear genetic groups. Grizzly bears sampled within an area represented by a given language family were significantly similar to those sampled within that language family (P = 0.001) and significantly divergent to those sampled outside the language family (P = 0.001). This spatial co-occurrence suggests that grizzly bear and human groups have been shaped by the landscape in similar ways, creating a convergence of grizzly bear genetic and human linguistic diversity. Additionally, grizzly bear management units designated by the provincial government currently divide an otherwise continuous group and exclude recently colonized island populations that are genetically continuous with adjacent mainland groups. This work provides not only insight into how ecological and geographic conditions can similarly shape the distribution of people and wildlife but also new genetic evidence to support renewed, locally led management of grizzly bears into the future.


INTRODUCTION
This work is dedicated to the late Nuaqawa (Evelyn Windsor), who taught us that people learned their language and way of life from the bears.
Landscape genetic studies have traditionally focused on a limited suite of geographic and anthropogenic landscape features within narrow temporal scales to explain patterns of genetic structure. Recent land cover and land use are the most common variables examined, appearing in 83% of surveyed studies that investigated landscape effects on functional connectivity, followed by roads and similar linear features in 24% of studies (Zeller et al. 2012). Although roads can create significant resistance to gene flow (Epps et al. 2005, Riley et al. 2006, water bodies can also exert similar influence (Frantz et al. 2010, Hartmann et al. 2013), yet they are rarely considered. Further, although the influence of modern landscape resistant features on gene flow can be significant, some genetic markers can also be informative for archaeological timescales (Brown et al. 2009, Harris andTaylor 2010). Indeed, anthropogenically shaped landscapes of the past can leave signatures in population structure observed today (Holzhauer et al. 2006), especially in species with long generation times (Spear andStorfer 2008, Zellmer andKnowles 2009).
Although anthropogenic settlement and land use are commonly included as resistance factors, human activities may not be categorically resistant, especially in historical timelines and considering less conflict-prone cultural approaches to coexistence (Clark andSlocombe 2009, Bhattacharyya andSlocombe 2017). For instance, whereas roads often pose barriers to gene flow, they are not equally and universally resistant (Frantz et al. 2012); indeed, some species use roads as movement corridors (Oleksa et al. 2015). Furthermore, historic human landscape use and settlement, even where dense, may not have the same resistance effects as modern anthropogenic use. For example, dense historic settlements had no effect on gene flow patterns in African elephant populations, whereas modern habitat selection patterns are best explained by contemporary settlement (Epps et al. 2013). The limited focus of many landscape genetic studies on modern land cover and linear features provides an opportunity for the inclusion of a broader suite of spatiotemporal geographic and anthropogenic variables. https://www.ecologyandsociety.org/vol26/iss3/art7/ Humans and human-human interactions not only shape landscapes (Power et al. 2018) but also can be shaped by the environment in similar ways to non-human animals. Broadly, this parallel response is so prevalent that it has resulted in global patterns of co-occurrence of biological diversity in wildlife and of linguistic diversity in human societies (Maffi 2005, Gorenflo et al. 2012). This co-location of biological (i.e., species) diversity and socio-cultural (i.e., cultural or linguistic) diversity occurs at local (Stepp et al. 2005), regional (Manne 2002), and global scales (Gorenflo et al. 2012). For example, in North America north of Mexico, high human linguistic, ethnolinguistic, and cultural diversity were all significantly correlated over space with high species diversity of trees (Smith 2001). As a component of cultural diversity, spatial patterns in linguistic diversity have been found to originate from long-term social interactions that relate to variation in a suite of biophysical factors, including geographic distance (Honkola et al. 2018), topography (Axelsen and Manrubia 2014), resource productivity (Codding and Jones 2013), and resource diversity (Moore et al. 2002). Although biophysical factors can influence and shape linguistic diversity over time, spatial patterns of linguistic diversity also reflect an enduring history of socio-cultural interaction and integration of communities of speakers. Our objective is to contribute to a small but growing body of work (e.g., Polfus et al. 2016, Schulz et al. 2019, Gagnon et al. 2020 investigating if and how spatial patterns in fine-scale components of the biocultural relationship, i.e., human linguistic diversity and wildlife genetic diversity, might align. Grizzly bears (Ursus arctos) in the Central Coast region of what is now known as British Columbia (BC), Canada, have shared habitat and resources with people on the landscape for millennia. Indigenous oral histories and songs from the region depict grizzly bears as kin and cultural or customary enforcers for their human neighbors (Housty 2014). Today, the Nuxalk, Haíɫzaqv, Kitasoo/ Xai'xais, Gitga'at, and Wuikinuxv Nations are increasingly engaging in Nation-led scientific research and management regarding land-use planning, fisheries management, and coexistence strategies with grizzly bears and black bears (U. americanus) in the region. Given the scholarly, management, and cultural context, the coastal grizzly bear system presents unique opportunities: (1) to examine the potential for genetic structure in a highly mobile species with a diversity of potential drivers over long time scales, (2) to investigate possible associations between genetic structure and Indigenous cultural diversity in the area, and (3) to assess the appropriateness of current management unit designations with regard to genetic diversity in a species of conservation concern.
A variety of factors has likely shaped biocultural variation on the Central Coast. In this region, industrial disturbance and road networks are limited, with a complex coastline comprising many channels, inlets, and estuaries that may supplant roads Stenhouse 2014, Proctor et al. 2020) as key linear features potentially affecting wildlife genetic structure. Additionally, the spatiotemporal influence of human activity on wildlife has been dynamic and complex in this region. In recent history, however, pronounced changes have occurred in human presence, activity, and interaction with wildlife. Within the most recent two centuries, violence-and disease-mediated genocide following European colonization reduced Indigenous populations by orders of magnitude and interrupted their role as prominent ecological agents within these coastal temperate rainforests (Harris 1994, McLaren et al. 2015. A conservative estimate of the population of the Northwest Coast of North America numbered > 180,000 people during the late 1700s, whereas less than a century later, only 35,000 people remained (Boyd 1999).
Despite these recent changes, such previously high-density societies, particularly when they had strong resource use overlap with bears over millennia, might have influenced bear movement and gene flow. The density of resources and high-relief topography of the Central Coast of BC in combination with the long-term establishment of Indigenous territories and autonomy, despite extensive cultural exchange, enabled a rich diversity of distinct Indigenous language families, languages, and dialects to develop in parallel over numerous generations (Nichols 1992, Beck 2000. Given broadly similar mobility and foraging ecology between bears and Indigenous peoples of the area over millennia, as well as evidence of co-locations of biological and cultural diversity elsewhere (e.g., Maffi 2005), we might expect similar spatial patterns defining bear genetic and Indigenous linguistic groups. Such an association could be mediated through two processes, which are not mutually exclusive: (1) combinations of biophysical factors resistant to people and bears created boundaries among both human linguistic and bear genetic groupings (Nichols 1992, Axelsen andManrubia 2014), and (2) broader spatial patterns in ecological heterogeneity (i.e., in the distribution of food and other resources) that may have, in part, structured humans and their language cultures also structured grizzly populations (Maffi 2005, Muñoz-Fuentes et al. 2009, Housty et al. 2014. Relationships between grizzly bears and humans in a contemporary management context also motivated our research collaboration and this specific work. We aimed to examine grizzly population genetic structure data to evaluate the spatial designation of grizzly bear population units (GBPUs) identified by the provincial government along the Central Coast of BC. Although the designation of GBPUs in southern BC considers the consequences of landscape and anthropogenic factors on genetic isolation, neither genetic diversity nor genetic structure are currently considered in the designation and management of the remote GBPUs of the Central Coast; boundaries are instead defined primarily by heights of land separating watersheds (Province of British Columbia 2012). Accordingly, evidence from genetic data could refine these management units to inform longterm management and conservation and ensure the viability of grizzly bear populations in this region. Such management is increasingly led by Indigenous Stewardship offices . https://www.ecologyandsociety.org/vol26/iss3/art7/ Against this background, we formed a number of analytical hypotheses to guide our investigation. Given the mobility of grizzly bears, the modest spatial scale of the study area, and the relatively homogenous resource distribution across the study area, we hypothesized that genetic structure analysis would reveal a semicontinuous population influenced by both landscape and anthropogenic variables (past and present). Specifically, we expected large waterways (Paetkau et al. 1998a), rugged terrain (Proctor et al. 2012, Lewis et al. 2015, and the presence of snow and ice (Lewis et al. 2015) to impose potential resistance to gene flow. Predicting responses to human activity in recent and more distant periods was more difficult. We hypothesized that contemporary human settlement (Proctor et al. 2012) and industrial forestry (Apps et al. 2004) could impose resistance to gene flow, but acknowledged that overall human density and industrial activity on the landscape today is limited. Similarly, we predicted that the dense and widespread presence of Indigenous peoples and communities over the past several millennia may have affected gene flow, but we recognize that cultural norms, which prioritize coexistence (Housty et al. 2014), might have limited such effects. Considering this complex background, our models included the possible resistance provided by archaeologically recorded, pre-industrial, Indigenous settlement sites and present human settlement and activity. Finally, given the aforementioned rich diversity of coastal Indigenous languages, shared ecological niche, potentially similar responses to biophysical factors, and kincentric relationship between bears and Indigenous peoples in coastal BC, we also hypothesized that there could be a spatial association between bear genetic groups and Indigenous language families. Hair snag stations were part of a long-term bear monitoring project in collaboration with the Nuxalk, Haíɫzaqv, Kitasoo/ Xai'xais, Gitga'at, and Wuikinuxv Nations, covering 23,500 km² of coastal BC. This multi-Nation collaboration that extends to partnering academia and of which this genetics research project is part, is a multifaceted, applied, and scholarly collection of research and management activities. Extensive background about how the voices of the communities, their history, and their relationships with grizzly bears and the environment influence this program is covered in depth elsewhere (Adams et al. 2014, Housty et al. 2014Artelle et al., unpublished manuscript). Broadly, this partnership, now a decade old, aligns with recent literature on biocultural research approaches that support equitable decision-making and connecting knowledge to actions (Groffman et al. 2010, Pretty 2011, Smith et al. 2013. Specifically and briefly, applied research on bears in the area was initiated with a project by a Haíɫzaqv nonprofit (QQS Projects Society) and led by a Haíɫzaqv scientist and manager (Housty et al. 2014 The genetics research described here was initiated because of interest expressed by several participating Nations in understanding the ancestry of grizzly bear individuals colonizing islands off the Central Coast (Service et al. 2014). Additionally, this project was intended to address the long-term monitoring goals of the Bear Working Group by identifying empirically informed spatial population units for management, given that existing Provincial GBPUs were created in the absence of local data. Accordingly, from project initiation to interpretation of results, the voices of participating community members are included in this scholarly form of dissemination, which complements important "inreach" and policy development activities of participating Nations.

DNA extraction and microsatellite genotyping
DNA extraction, amplification, and sequencing were conducted by Wildlife Genetics International in Nelson, BC. We extracted DNA from selected hairs with a Qiagen DNeasy Blood and Tissue Kit (Qiagen, Maryland, USA) using the spin column protocol with a double elution (for detailed extraction methods, see Appendix 1).

Population genetic analyses
We estimated observed and expected heterozygosity per population using the R package "adegenet" ; see Appendix 1 for tests for Hardy-Weinberg proportions). To identify the spatial pattern of grizzly bear population structure, we calculated the mean center of all detections per individual and employed one frequentist and two Bayesian methods. We first conducted a spatial principal components analysis (sPCA), which is a frequentist, non-model based approach that does not rely on Hardy-Weinberg assumptions conducted in R version 3.2.4 (R Core Team 2018) using the "adegenet" package with a Delauney triangulation, with eigenvalues chosen interactively and data not scaled to unit variance. To complement the frequentist, spatially explicit PCA, we also implemented a model-based Bayesian clustering algorithm without a priori populations, because of no prior knowledge of genetic groupings, using the software STRUCTURE  with the parameters USEPOPINFO = 0, MAXPOPS = 6, NOADMIX = 0, ALPHA = 1, with 10 iterations of 50,000 burn-ins and 450,000 Markov chain Monte Carlo repetitions. These parameters were selected following the guidance of Porras-Hurtado et al. (2013) in accordance with our sample numbers, loci set, possible estimated maximum number of populations given the relatively small geographic area, and realistic probability of admixture among genetic groups. To account for potential unbalanced sampling, we conducted an additional analysis by adjusting alpha to 0.17 (Wang 2017). We determined the appropriate number of populations (K) using Delta K, a second-order change rate in likelihood that peaks at the most probable K as implemented in Structure Harvester (Evanno et al. 2005, Earl andvonHoldt 2012). To represent STRUCTURE results visually, we interpolated population Q-values (individual cluster membership proportions) using the inverse distance weighting (IDW) method, which uses the assumption of spatial autocorrelation to interpolate values across unmeasured areas, in ArcGIS 10.3. We combined IDW surfaces for all populations using the Composite Bands function (Chiocchini et al. 2016). We used Geneland version 4.0.6  as an additional Bayesian spatially explicit clustering program with default parameters to assess up to 10 genetic partitions with no uncertainty on coordinates, 100,000 iterations, 100 thinning, and an uncorrelated allele frequency model. We did not calculate F ST (or analogs) because of the potential that our microsatellite markers and observed populations do not follow the assumptions of a stepwise mutation model (Meirmans and Hedrick 2011) and have a higher migration than mutation rate (Wright 1943, Balloux andLugon-Moulin 2002). Because of their geographic isolation, we removed 22 CID samples from Kitimat from further analyses, resulting in 128 samples (N = 125 noninvasive hair samples, N = 3 CID samples) for connectivity and resistance estimation. We imposed a map of GBPUs onto the map of the IDW-interpolated STRUCTURE Q-values to assess visually the extent of overlap between genetic group and GBPU boundaries.

Landscape variables
Based on the literature (Nichols 1992, Paetkau et al. 1998a, Apps et al. 2004, Proctor et al. 2012, Axelsen and Manrubia 2014, Lewis et al. 2015, we created resistance layers for the variables we expected to affect grizzly bear gene flow. All layers are in the EPSG:3005 coordinate system with the NAD83 datum and BC Ecology and Society 26(3): 7 https://www.ecologyandsociety.org/vol26/iss3/art7/ Albers Equal Area projection. The layers included terrain ruggedness, waterways (Geospatial Data Downloads 2018), ice and snow (Government of Canada 2015), archaeological (preindustrial) anthropogenic settlement as indicated by the presence of recorded shell midden deposits and fish traps documented during archaeological survey efforts (primary sources in Table  A1.2 in Appendix 1), and modern anthropogenic disturbance in the form of a simple variable indicating either presence or absence of forestry activity (Hermosilla et al. 2018), as well as current settlement defined by pockets of human-dominated areas (Province of British Columbia 2016). We acknowledge that recorded archaeological sites are subject to the uneven distribution of survey effort, do not encompass all archaeological research in the area, and likely represent only a small proportion of actual habitation and use.
Although harvesting of grizzly bears occurred on the Central Coast in historic and modern times, we infer from data described below that its influence on genetic structure was minimal, and thus, we did not consider its potential effects in our analyses. An analysis of animal remains recovered in archeological sites suggests that black bears were killed during the Pleistocene-Holocene transition in Haida Gwaii (McLaren et al. 2005), and the only excavated archaeological site in the study area with identified bear bones contained a small number of remains from black bears only (Cannon 1991). In the 1800s, the Hudson Bay Company had trading posts on the coast. Although grizzly bears were clearly killed by commercial hunters and trappers, they were not a preferred target (Meilleur 2002). During much of the past century, when hunting was managed by the Province, only a relatively small part of our study area (the Bella Coola valley) showed evidence of recent overexploitation (Artelle et al. 2013). This limited contemporary harvest on the Central Coast is likely associated with restricted access because of the near absence of roads. Given this extended history of intermittent and generally low levels of killing, we reason that hunting has had minimal influence on genetic structure compared to continuous factors such as landscape resistance and habitat fragmentation (McLellan et al. 2019).
Recognizing a visual overlap between Indigenous language family boundaries and the identified bear genetic groups, we also considered the approximate physical location (i.e., linear extent) of the boundaries of Indigenous language families on the Central Coast as representative of a potentially unique combination of biophysical factors that may be resistant to both people and bears (Nichols 1992, Axelsen andManrubia 2014). These boundaries did not spatially align completely with any other single considered resistance factor, supporting their inclusion as a unique variable. We considered multiple sources, which generally align, to define extents and boundaries among language families. The boundaries we used (digitized from Lepofsky and Turner 2013) are largely consistent with the First Peoples' language map of BC (First Peoples' Culture Council of British Columbia 2020), and are based on older linguistic renderings of relationships among language families with updated spellings and place names (Thompson andKinkade 1990, Beck 2000). As background to their estimation and in the absence of written records, modern and historic relationships among individual languages and language families have been gleaned through analyses of typological and grammatical similarities (Thomason and Kaufman 1988). The Salishan language family (here referred to as Salishan Nuxalk due to the inclusion of It7Nuxalkmc or Nuxalk in English, the language of the Nuxalk Nation) diverges strongly from the two other language families in the study area, i.e., the Tsimshian and Wakashan language families, which were identified as fundamentally distinct linguistic groups from the late 1800s (Tolmie andDawson 1884, Boas 1887); relationships among other language families on the coast were refined and debated through later research (Swadesh 1951, Thomason andKaufman 1988, Thompson andKinkade 1990). Although classified only relatively recently by linguists, language family boundaries on the Central Coast have a deep history (Embleton 1985), which is supported by archaeological and oral history evidence of long-term human occupation (Cannon 2003, McLaren et al. 2015. The categorization of this area as part of a larger residual zone (an area defined by diverse coexisting languages maintained by an established intergroup equilibrium) provides some evidence that the pattern of linguistic distribution may be older because of the long time periods typically needed to form an equilibrium and foster grammatical autonomy in language groups (Nichols 1992). These language families have considerable meaning and place-specific relevance to Nations participating in this work, with Gitga'at and Kitasoo members speaking languages in the Tsimshian family; Haíɫzaqv, Wuikinuxv, and Xai'xais Nations using languages encompassed by the Northern Wakashan family; and the Nuxalk Nation represented by the Salish-isolate Salishan Nuxalk family. Although spellings can differ within and among languages, grizzly bears are known as Nan in the Haíɫzaqvla, It7Nuxalkmc, 'Wuiḱala, and Xai'xais languages spoken by the Haíɫzaqv, Nuxalk, Wuikinuxv, and Xai'xais Nations, respectively. In the Sgüüxs language of the Kitasoo Nation and the Sm'algyax language of the Gitga'at Nation, grizzly bears are known as Medi'ik.
The delineation of language family boundaries in this study was made with the understanding that these boundaries are fuzzy, imperfectly resolved, and potentially reflect a relatively limited temporal depth. Additionally, human travel and exchange among neighboring linguistic groups was common, and many individuals likely had fluency in multiple languages (Ames 2002). Given this information, we used analysis methods focused on either a rough approximation of these language family borders or the entire spatial area encompassed by each of the three language families. We additionally emphasize that the focus of this work is on the identification of bear genetic patterns and the factors potentially underlying them. Although we identify spatial alignment between language families and bear genetic groups, the focus on bear genetic data limits our ability to comment extensively on the derivation of language family boundaries.

Connectivity and resistance estimation
We assigned categorical resistance values to variables on a scale of 1 (very low) to 5 (maximum; Table 1), which assisted in our ability to compare resistance surfaces. Because of the potential importance of waterways on the coast, we tested an additional parameterization scheme with all waterways assigned a maximum resistance value of 5 and all other landscape cells assigned a value of 1. We used the circuit theory-based software program Circuitscape   Lewis et al. (2015). ‡ Paetkau et al. (1998a). To validate the resistance surfaces, we constructed a Euclidean distance matrix of Q-values assigned by STRUCTURE as the dependent variable of genetic distance (Balkenhol et al. 2014). The resulting genetic distance matrix was used in the statistical validation of matrices of effective distances created by the input of resistance surfaces into Circuitscape. Because of high multicollinearity between significant variables identified by Mantel tests (Appendix 1), we used commonality analysis (CA), which employs variance partitioning to address collinearity issues (Prunier et al. 2015). Although CA parsed out the unique and common contribution of each variable, the high collinearity of the top two landscape variables (one being geographic distance) made interpretation of their respective contributions difficult. This result prompted us to remove the apparent effect of isolation by distance in our data set to avoid erroneously inflating the explanatory power of the other landscape variables.
Finally, to partial out the effect of geographic distance, we calculated the residuals from a linear regression of each landscape variable using the lm and residuals functions in the R package "stats". This procedure allowed us to identify the unique contribution of variables to the pattern of genetic structure while accounting for the contribution of geographic distance (i.e., isolation-by-distance). The variable residuals were combined using CA with logistic regression on distance matrices (LRDM) with the CAlrdm function, with 1000 bootstrap iterations and 10,000 permutations for tests of significance. Within the CAlrdm framework, we tested two separate models: an "archaeological" resistance surface that included landscape (waterways, snow and ice, terrain ruggedness) and anthropogenic (middens, fish traps, Indigenous language family boundaries) variables, and a "modern" surface containing the same landscape variables plus modern anthropogenic variables (human dominated areas, forestry).

Statistical analysis of language family overlap with bear genetic groupings
In addition to testing the hypothesis that the linear extents of human language family boundaries represent potentially resistant areas to grizzly bear gene flow, we tested for the overall overlap and similarities in spatial structuring among Indigenous language families and bear genetic groups using analysis of similarities (ANOSIM) and multivariate analysis of variance (MANOVA). We chose to use broader Indigenous language families rather than individual languages and territories because of their deeper temporal scale and more stable nature (Beck 2000, Lepofsky andTurner 2013). Using this map, bears were assigned to language family using their mean center of detection. We tested if bears within each language group had more similar Q-values to other bears within that language group or to members of another language group. We analyzed three language family categories and the Euclidean distance matrix of STRUCTURE Q-values for ANOSIM within the R package "vegan". To complement ANOSIM, we used MANOVA, which uses Q-values representing each bear genetic group assignment and language family categories to test for dissimilarities in Q-values between each language family. We did not remove the effect of geographic distance from these analyses because their purpose was to identify similarities in structure between bear and human groups, not to determine the contribution of these language groups to the pattern of bear genetic structure or vice versa. Additionally, given the innate importance of geographic distance in structuring Indigenous languages (Creanza et al. 2015), the removal of geographic distance patterns would likely remove important variation in language grouping. Furthermore, significant spatial overlap of bear genetic groups and human language groups may be due to a convergent response to the effects of geographic distance.

Central Coast grizzly bear genetic groups
Grizzly bears on the Central Coast form three distinct genetic groups. STRUCTURE identified K = 3 (G1, G2, G3) as the most likely number of genetic partitions, with high average cluster assignments (Q G1 = 0.92, Q G2 = 0.91, Q G3 = 0.87; Fig. 2; Fig. A1.1 in Appendix 1). We also identified K = 3 as the best fit partition when we implemented STRUCTURE with alpha = 1/k. Interpolation of STRUCTURE Q-values aided in visually identifying the geographic locations of boundaries among genetic groups (Fig. 2). We corroborated this genetic and geographic structural pattern by the additional Bayesian spatially explicit clustering algorithm implemented in Geneland (Fig. A1.2 in Appendix 1). Geneland confirmed the presence of three genetic groups identified by STRUCTURE and the locations of breaks among them. The sPCA also identified boundaries in population structure concordant to those identified by Bayesian approaches (Fig. A1.3 in Appendix 1). We did not identify any difference in sPCA structural patterns when samples were separated by sex (Figs. A1.4 and A1.5 in Appendix 1). Collectively, we identified three distinct genetic groups with clearly defined boundaries appearing primarily along waterways (Fig. 2). When we investigated the provenance of island individuals, we found that most bears detected on islands (N = 18 of 24; 75%) belonged to the nearest mainland genetic group, which is consistent with recent field observations and Indigenous knowledge, as well as inferences that grizzly bears have been actively colonizing islands in the recent past (Service et al. 2014).

Fig. 2.
Plots of grizzly bear genetic and Indigenous language family overlap. Plots are STRUCTURE bar plot at K = 3 (A) and combined (B) and separate (C) inverse distance weightinginterpolated STRUCTURE population Q-value maps (blue = G1, green = G2, red = G3). Dark grey lines indicate borders among Indigenous language families (approximate overlap: blue = Tsimshian, green = Wakashan, red = Salishan Nuxalk); light grey points indicate mean detection locations of individual grizzly bears.
We also examined the overlap between the most likely number of genetic partitions and the geography of BC provincial GBPU designations (Province of British Columbia 2012). An overlay of the genetic partitions onto the IDW-interpolated map of population genetic structure revealed some inconsistencies (Fig.  3). We identified an overlap with the boundary between G2 and G3 and the GBPU boundary, but a mismatch with the boundary between G1 and G2. This mismatch results in the split of the otherwise continuous G2 group, with 16 of 58 (28%) individuals assigned to G2 isolated within a separate GBPU.  (Table A1.3 in Appendix 1).

Explanatory power of resistance surfaces
Resistance surfaces did not explain a significant proportion of genetic variation. Within a CA-LRDM framework with geographic distance included as a variable, the archaeological surface explained similar variation in the genetic distance matrix (23.24%) to the modern surface (23.18%), with water (P = 0.001), language boundaries (P = 0.05), and geographic distance (P = 0.001) identified as significant variables in the archaeological surface. Forestry (P = 0.046), waterways (P = 0.001), and geographic distance (P = 0.001) were significant in the modern surface. Although significant, language boundaries had a low unique contribution (u = 0.002) to the variation in the genetic https://www.ecologyandsociety.org/vol26/iss3/art7/ Fig. 4. Plots of explained variance (%) for combinations of all tested variables. (A) Commonality coefficients for the archaeological surface, with 95% bootstrap confidence intervals (F1 = language family residuals, F2 = fish trap residuals, F3 = midden residuals, F4 = ice and snow residuals, F5 = waterway residuals, F6 = ruggedness residuals). (B) Commonality coefficients for the modern surface, with 95% bootstrap confidence intervals (F1 = modern settlement residuals, F2 = forestry residuals, F3 = ice and snow residuals, F4 = waterway residuals, F5 = ruggedness residuals). distance matrix. The independent effects of water and geographic distance were difficult to disentangle because of high multicollinearity between these variables (R = 0.901). We observed similar multicollinearity between the water and geographic distance variables in the modern surface. Although significant, forestry had both a negative beta weight (ß = −0.053) and a low unique contribution (u = 0.003). Given the apparent contribution of isolation by distance in this data set, we removed geographic distance from all variables using linear regressions and used Mantel tests to detect correlations among matrices. Finding no substantial correlations (R > 0.7), we included all variables in a modern and archaeological surface framework within a CA-LRDM analysis. With geographic distance removed, the archaeological surface explained only 2.44% of the variation in STRUCTURE Q-values, and the modern surface explained a similarly small percentage (2.40%). Language boundaries, however, remained significant (P = 0.001) in the archaeological surface in addition to forestry (P = 0.001) in the modern surface (Fig. 4). Waterways also remained significant in both modern and archaeological models (P = 0.001), but the issues of low unique contributions for language (u = 0.005) and forestry (u = 0.005), and the negative beta weight of forestry (ß = −0.039; Fig. A1.6 in Appendix 1) called into question their contribution to genetic structure. Although water continued to be important, it explained a negligible proportion of variation in genetic variation (1.93% in archaeological and modern surfaces). To determine if this low explanatory power was due to the parameterization scheme, we tested a model in which all waterways, regardless of width, were assigned maximum resistance. However, waterways continued to have only marginal explanatory power for genetic distances (0.19% in archaeological and modern surfaces).

Spatial overlap of grizzly bear genetic groups and Indigenous language families
Given the finding of significant resistance but very low explanatory power provided by language group boundaries, we found little support for the hypothesis that the relationship between language and bear genetic groupings was mediated by a combination of potential biophysical resistance features present in these border regions. However, the lack of high collinearity (R > 0.7) between language group boundaries and other landscape variables supported the inclusion of language group borders as a unique potential resistance variable. To test the alternate hypothesis that spatial patterns in ecological heterogeneity (i.e., in the distribution of food and other resources) may have resulted in convergent population structure for bears and people, we used ANOSIM and MANOVA to investigate similarities in Q-values. We found that bears within the spatial boundaries of distinct language families were both significantly similar to each other Ecology and Society 26(3): 7 https://www.ecologyandsociety.org/vol26/iss3/art7/  [Lepofsky and Turner 2013]). STRUCTURE Q-values belonging to G1 occurred within the spatial boundary of the Tsimshian language family. STRUCTURE Q-values belonging to G2 occurred within the spatial boundary of the Wakashan language family. STRUCTURE Q-values belonging to G3 occurred within the spatial boundary of the Salishan Nuxalk language family. Boxes represent the interquartile range (middle 50% of range) with a line corresponding to the median, with whiskers displaying 1.5 box height (limited by maximum and minimum data points) above and below the box.
(with some overlap; ANOSIM R = 0.403, P = 0.001) and significantly dissimilar to bears located within the geographic areas of other Indigenous language families (MANOVA P = 0.001; Fig. 5). Thus, Indigenous language families and bear genetic groups showed significant spatial overlap.

DISCUSSION
Although it is increasingly apparent that ecological systems reflect long-term human use and cultural influence (e.g., Levis et al. 2017), the legacies of social-ecological systems and biocultural diversity are rarely considered in wildlife landscape genetic analyses. Our finding of spatial convergence between bear genetic groups and Indigenous language families not only offers a new case study to a body of literature describing the relationship between cultural and biological diversity (Maffi 2005, Bridgewater andRotherham 2019), but also contributes to emerging work on a finer scale of biocultural relationships (i.e., the co-localization of linguistic diversity in humans and genetic diversity in co-occurring wildlife; e.g., Polfus et al. 2016, Schulz et al. 2019, Gagnon et al. 2020. We also showcase here Indigenous-led applied research that provided the first landscape genetics data to inform spatially explicit grizzly bear management on the Central Coast of BC, an area in which First Nations are reasserting their agency to continue managing wildlife and other resources. We found three clearly distinct grizzly bear genetic groups corroborated by three analysis methods. Analyses using resistance surfaces for both modern and archaeological landscape and anthropogenic variables, however, revealed no significant or biologically meaningful contribution of these variables to this landscape genetic pattern. Of the variables tested, geographic distance played a substantial role, followed by waterways. Although significant, the low explanatory power of both water and Indigenous language family boundaries indicates that their role as resistance factors is slight and likely not biologically relevant for the overall process of bear gene flow in this landscape. Though border areas among language families do not appear to be highly resistant, the finding of spatial genetic structure cooccurring with the spatial distribution of language families suggests that there is not extensive gene flow across these boundary areas. Given this finding, we suggest that the spatial areas encompassed by language families may have experienced a higher degree of human-human interaction within rather than between these regions, as well as long-term coexistence with grizzly bear populations that are similarly structured. These observations appear to contrast with historic examples from less geographically complex environments, where persistent human conflicts between adjacent groups is hypothesized to facilitate associated wildlife refugia along conflict zones (Martin andSzuter 1999, Laliberte andRipple 2004). Additionally, the lack of substantial collinearity between language family boundaries and the biophysical variables of ruggedness and snow and ice suggests that topography is possibly not the primary driver of linguistic differentiation in this area. Collectively, these patterns and inferences suggest the role of an underlying influence of isolation by distance in grizzly bear genetic structure and the potential contribution of other unexamined landscape, humanrelated, and demographic variables. Specifically, female natalbiased dispersal (Shirane et al. 2018) may contribute to the observed pattern of population structure, although we did not see differences between structural patterns in males and females. Additionally, the pattern might be driven by a response to spatial variation in resource distribution at a smaller spatial scale than we could assess.
The primary limitations of this study include the selection and parameterization of only a handful of potential variables, and the inherent limitations associated with the assumptions of https://www.ecologyandsociety.org/vol26/iss3/art7/ parameterization (Zeller et al. 2012). However, we selected variables from relevant literature. Because we used expert opinion guided by relevant literature to parameterize resistance surfaces, we might have under-or overestimated the influence of considered resistance factors facing grizzly bears on the coast (Zeller et al. 2012). Our parameterization scheme, also informed by available literature, is limited. We did conduct a coarse sensitivity analyses with the waterways variable, but the patterns remained unchanged. Additionally, although we used the coarser scale of language family rather than individual Indigenous language dialects, we acknowledge that there remains uncertainty surrounding the precise boundaries and temporal stability of language family boundaries ( Thomason and Kaufman 1988). Although we expect some influence of landscape variables on the distribution of Indigenous language families, we acknowledge that socio-cultural and historical factors were likely as or more important in shaping the geographic distribution of language families. Finally, difficulty predicting the temporal lag between the creation and observation of these genetic groups limits the understanding of whether archaeological or modern anthropogenic features do not have an effect on genetic structures or are simply asynchronous. If this uncertain temporal lag is longer than typically expected with microsatellite markers in a highly mobile species (maximum 15 generations; Landguth et al. 2010), postglacial sea-level change and the presence of snow and ice could provide different levels of resistance (e.g., Mackie et al. 2011Mackie et al. , 2018. Additionally, this bear genetic structuring pattern, if older than anticipated, may also in part reflect a history of glacial refugia, with genetic groups created through a series of isolation and recolonization events (Waits et al. 1998).
The slightly higher explanatory power of the archaeological surface, the co-localization of bear genetic groups with longstanding Indigenous language families, and the observations of current long-distance bear movements collectively suggest that the grizzly bear genetic structure pattern is at least, in part, historical. Given this inference, the lack of association between human settlement and activities (on both preindustrial and current time frames) and bear genetic structure, despite high archaeological human densities and with year-round intensive landscape use , Maxwell et al. 1997, Boyd 1999, Cannon 2003, is somewhat surprising. Other work that considers contemporary human settlement and activities has found that these forces can have pronounced disruptive effects on wildlife, including bears (Proctor et al. 2012, Støen et al. 2015. This absence of a relationship could relate to the previously mentioned asynchronous nature of the grizzly bear genetic structure pattern and archaeological human settlement or the relatively low and variable archaeological survey effort conducted in this region thus far. Alternatively, this result could reflect an Indigenous placebased management model of respect for and reciprocity with bears (Housty et al. 2014) that prioritizes coexistence over dominion, despite considerable resource niche overlap and potentially mediated through avoidance and domestic dog guardians (McKechnie et al. 2020). Such possibilities emphasize the utility of considering a broad suite of landscape parameters, spatial and temporal scales, and human cultural norms toward wildlife when investigating landscape genetic patterns. The assumption of resistance to human activity might be less likely to hold when considering (even dense) human populations of the recent past and varied approaches to coexistence.
The geographic locations of the three observed genetic groups with regard to GBPU designations, the genetic groups of island individuals, and the boundaries of Indigenous language families carry management implications. The finding that the boundary between G1 and G2 was not reflected in GBPU boundaries provides an opportunity for consideration of genetic data in refinements to the designation of GBPUs. Considering that GBPU distinctions are used to set population targets and for landuse planning and priorities (Province of British Columbia 2012), erroneous subdivision of the continuous genetic group with the lowest heterozygosity (G2) could lead to incorrect inference regarding the long-term viability of G1 and G2 and the movement and mating patterns of grizzly bears on the landscape. Additionally, our finding that most island individuals (75%) belong to the adjacent mainland genetic group supports the inclusion of these recently colonized island populations (Service et al. 2014) within the nearest GBPU. The identification of three genetic groups on the Central Coast suggests that low heterozygosity and population substructuring in grizzly bears are not unique to southern BC. Moreover, owing to the inability of ruggedness to explain the pattern of differentiation, the heights of land between watersheds may not be appropriate for capturing genetic structure on the Central Coast. Future directions to further reassess the delineation of Central Coast GBPUs include the use of whole nuclear and mitochondrial genome markers to assess levels of adaptive differentiation among these groups. Additionally, the co-localization of Indigenous language families and grizzly bear genetic groups emphasizes the need for continued locally led management of this culturally important species, which will likely include population unit delineation by the sovereign Indigenous Nations with the authority and agency to do so.
Co-localization of the spatial extent of Indigenous language families and bear genetic groups implies a potentially similar relationship between each and the landscape, with both producing comparable structural patterns when confronted with the constraints of geographic distance. The strong genetic structure present in grizzly bear populations within the relatively small scale of Central Coast area is consistent with the linguistic diversity present in the area created through the long-term enduring coexistence of distinct languages and communities over time (Beck 2000). The parallel nature of these structures may indicate more than solely similar responses to geographic distance. Specifically, each language family and grizzly bear genetic group may be occupying its own broad resource space with allowance for fine-scale ecological heterogeneity and use by many autonomous Indigenous communities within a single language family (Beck 2000). The degree of either Indigenous language convergence or grizzly bear genetic admixture may be mediated by the intensity of contact among neighboring groups and also the opportunities afforded by the abundant resources on the landscape to support local grizzly bear genetic and human linguistic variation. Although this co-localization does not explicitly show any relationship between grizzly bear and human groupings, it does imply that the same landscape pressures that shaped Indigenous language families, should they be resource or geographically mediated, also could have shaped grizzly bear genetic groups. The spatial convergence of language groupings and grizzly bear genetic groups comes as no surprise to Elders and community members in the territories in which our research Ecology and Society 26(3): 7 https://www.ecologyandsociety.org/vol26/iss3/art7/ collaboration occurs, where people and grizzly bears have shared space and resources for millennia. More broadly, this finding contributes to a body of biocultural research that describes the spatial co-occurrence of cultural and linguistic diversity with biological diversity across many different environments and ecosystems (Maffi 2005, Gorenflo et al. 2012. Specifically, our finding of spatial overlap between grizzly bear genetic and human linguistic groups, and the implications of a similar response to biophysical factors, situates this work within scholarship investigating how global patterns of biocultural relationships can manifest at finer scales of genetic and linguistic diversity.
Responses to this article can be read online at: https://www.ecologyandsociety.org/issues/responses. php/12443

Acknowledgments:
We appreciate and recognize the contribution and collaboration of > 100 people to the long-term collection of grizzly bear hair samples, with special gratitude extended to our Gitga'at, Haíɫzaqv, Kitasoo/ Xai'xais, Wuikinuxv, and Nuxalk leadership and Guardian Watchmen programs, as well  We prepared hair samples for extraction by selecting 10 guard hairs when available and 6 supplementing with five underfur hairs per missing guard where needed. When no guard 7 hairs were present in a sample, 30 underfur hairs were chosen. 8 9 For microsatellite genotyping, we labeled primers with FAM, HEX, HEX, TET, or NED dye 10 groups. We amplified DNA on a MJ Research PTC-100 thermocycler with PCR reagent 11 concentrations optimized for each primer pair (Table A1.1). We utilized a quality control 12 protocol that involved subsampling each sample and removing samples that were poor quality or 13 had three or more alleles at a locus (Paetkau 2003). This protocol has previously resulted in error 14 rates of 0.002-0.005 per locus per sample (Kendall et al. 2009).

16
The amelogenin locus for sex determination was amplified using 10 pM of each primer 17 (Forward:CAGCCAAACCTCCCTCTGC Reverse:CCCGCTTGGTCTTGTCTGTTGC), 200uM 18 dNTPs, and 0.9 units of Taq polymerase on a MJ Research PTC-100 thermocycler. We 19 distinguished between male and female individuals using gel electrophoresis with female sample 20 producing a single 280bp fragment and male samples producing a 280 and a 217 bp fragment. To 21 avoid Y allele dropout, we only sexed samples that produced high confidence scores for all other 22 microsatellite loci (Paetkau 2003). This method has previously produced error rates of 0.0007 23 per locus per sample (Kendall et al. 2009).

25
Hardy Weinberg proportions 26 27 Deviations from Hardy Weinberg equilibrium (HWE) were investigated with the web-based 28 version of Genepop (Rousset 2008). Identifying deviations is important as they can provide 29 information on population size, gene flow, and the presence of selection (Allendorf et al. 2012).

30
Additionally, HWE is an essential assumption underlying the model-based clustering algorithms 31 performed by STRUCTURE and Geneland  Global deviation from HWE found in one genetic group 42 43 We found significant (p < 0.05) heterozygote deficiency at three loci (G10C, p = 0.017; G10L, p 44 = 0.026; and MSUT2, p = 0.030 in population G3, G2, and G1 respectively). We observed 45 significant (p < 0.05) heterozygote excess at four loci: G10B (p = 0.028) and G10U (p = 0.040) 46 in population G1, and MU59 (p = 0.023) and X145P07 (p = 0.030) in population G3. Using 47 global HWE tests, we did not find any populations with a significant heterozygote deficit, 48 whereas G1 was the only population that showed significant (p < 0.050) global deviation from 49 HWE in the form of heterozygote excess (p = 0.015). This heterozygote excess suggests that G1 50 may be receiving gene flow from other areas or represents the unification of previously separated 51 populations (Allendorf et al. 2012). Though these deviations from HWE can be problematic for 52 the use of STRUCTURE and Geneland, deviation is only globally present in one genetic group 53 and the breaks between groups are confirmed with sPCA, a method that does not require the 54 assumptions of HWE to be met . 55 56